Substrate-mediated pattern formation in monolayer transfer: a reduced model

The formation of regular stripe patterns during the transfer of surfactant monolayers from water surfaces onto moving solid substrates can be understood as a phase decomposition process under the influence of the effective molecular interaction between the substrate and the monolayer, also called substrate-mediated condensation (SMC). To describe this phenomenon, we propose a reduced model based on an amended Cahn–Hilliard equation. A combination of numerical simulations and continuation methods is employed to investigate stationary and time-periodic solutions of the model and to determine the resulting bifurcation diagram. The onset of spatiotemporal pattern formation is found to result from a homoclinic and a Hopf bifurcation at small and large substrate speeds, respectively. The critical velocity corresponding to the Hopf bifurcation can be calculated by means of the marginal stability criterion for pattern formation behind propagating fronts. In the regime of low transfer velocities, the stationary solutions exhibit snaking behavior.


Introduction
Coating solid substrates with regularly patterned surfactant monolayers is exemplary for the purposeful utilization of self-organized patterns formation. Using the Langmuir-Blodgett technique, i.e. withdrawing a substrate from a water-filled trough covered by a monolayer of the phospholipid dipalmitoylphosphatdylcholine (DPPC), transfer of highly regular patterns with periods down to a few hundreds of nanometers can be achieved [1][2][3][4][5][6][7]. The formation of these structures is directly connected to the thermodynamics of the phospholipid monolayers and in particular to the first-order phase transition between the liquid-expanded (LE) and liquidcondensed (LC) phases. In fact, the observed patterns consist of alternating LE and LC domains. The domains form stripes that are either parallel or perpendicular to the water-substrate contact line [3]. The phase transition distinguishes these patterns from similar patterns formed in the Langmuir-Blodgett transfer of fatty acids that does not involve any phase transition [8][9][10][11]. Remarkably, before the transfer the floating monolayer is in the pure LE phase. Thus, the partial condensation necessary for the formation of the observed patterns has to occur during the transfer process itself. In fact, it results from an effective molecular interaction between the substrate and the monolayer and is commonly referred to as 'substrate-mediated condensation' (SMC) [12][13][14][15][16]. In the present paper, we argue that the observed pattern formation can be understood as a phase decomposition process under the influence of the SMC interaction field. This implies that the formation of stripe patterns observed in the Langmuir-Blodgett transfer of monolayers of binary phospholipid mixtures, which is ascribed to phase separation [17,18], should be amenable by a similar theoretical approach.
A recently developed model describes the monolayer transfer process using a set of two coupled nonlinear partial differential equations for the time evolution of the height of the liquid sub-phase and the density of the covering surfactant layer [19,20]. The hydrodynamics of the thin liquid layer is treated within the framework of lubrication theory [21] and is coupled to the monolayer thermodynamics via a surfactant equation of state. This equation quantifies the dependence of the surface tension of the liquid-air interface on the density of the surfactant monolayer and results in variations of the Laplace pressure and a Marangoni force. It can be derived from a free energy functional of the type proposed by Cahn and Hilliard [22] in the context of phase decomposition when applied to describe monolayers in the vicinity of a firstorder phase transition [23]. The results obtained using this model indicate that in this case the hydrodynamics in the meniscus is not essential for the formation of stripe patterns. Instead, the patterns form by phase decomposition under the influence of an external field that describes the interaction between the substrate and the monolayer [19,20]. Therefore, this effect may be captured by a simplified model that allows for a more profound mathematical analysis and, as a consequence, a deeper understanding of the physical mechanisms. Here, we propose to describe the process using a variant of the prototypical model equation for phase decomposition, namely the Cahn-Hilliard equation [24]. Within this reduced model the SMC effect is incorporated via a spatially varying free energy and the moving substrate is modeled by an advective transport term. Similar problems have been investigated by Krekhov [25] as well as Foard and Wagner [26]. Both consider the pattern formation behind a propagating quench front which results, within certain parameter ranges, in similar structures to those observed here. This implies that the analysis that we develop below may be adapted to their systems.
Following a brief outline of the previously developed full thin-film model of monolayer transfer in section 2, the reduced model is introduced in section 3 together with a presentation of the resulting basic two-dimensional (2D) deposition patterns. Subsequently, in section 4 we focus on stripe patterns parallel to the contact line and analyze various steady and time-periodic solutions of the reduced model equations in one dimension and summarize the results in a bifurcation diagram. A detailed analysis of the bifurcations leading to the onset of the formation of stripe patterns at low and high velocities is given in sections 4.2 and 4.3 respectively. 4 of fresh water dragged from the reservoir defines a meniscus and a contact line region that are steady in the laboratory system but move relative to the substrate. As in the corresponding Langmuir-Blodgett experiments [1][2][3][4][5][6][7], the water reservoir is assumed to be covered with a surfactant monolayer in the pure LE phase. The flow carries the surfactant towards the contact line where it is eventually subjected to the interaction with the substrate so that SMC sets in. The employed geometry of the liquid meniscus is similar to that used in the deposition experiments with polymer solutions of Yabu and Shimomura [27].
Here, we describe the liquid film by its height profile h(x, t), which indicates the local film thickness, whereas the surfactant density at the surface above the point x = (x, y) is described by the function γ (x, t). By the choice of a proper set of scales h 0 , l 0 , t 0 , γ 0 and σ 0 for height, length, time, surfactant density and surface tension, the time evolution equations for the height profile H (X, T ) := h(Xl 0 , T t 0 )/ h 0 and the surfactant density (X, T ) := γ (Xl 0 , T t 0 )/γ 0 can be written in dimensionless form as [20] with the generalized pressureP = −(1 − 2 P hom ( ))∇ 2 H + (H ) comprising the Laplace pressure ∼ ∇ 2 H , which reflects the relation of the surface curvature to the pressure difference between the liquid and vapor, as well as the disjoining pressure (H ) that describes the wetting properties of the substrate, i.e. the substrate-liquid interaction [28,29]. Here, P hom ( ) is the surface pressure of a homogeneous surfactant layer of density . The spatial variations of the surfactant density lead to surface tension gradients that induce Marangoni forces. They are accounted for by the terms involving ∇ˆ = −∇ P hom ( ) + −2 ∇ 3 . Here, the dimensionless number = σ 0 /(γ 2 0 κ)(a 6 /a 3 ) 1/3 relates the forces acting at the monolayer covered surface, characterized by σ 0 , γ 0 , and the surfactant domain line tension κ to the forces acting between the substrate and the liquid, characterized by the corresponding Hamaker constants a 3 , a 6 [20]. Both equations (1) and (2) contain an advective term proportional to the transfer velocity V that describes how the liquid and surfactants are dragged by the moving substrate.
The non-mass-conserving contribution E ev δµ to equation (1) arises from the evaporation of the water film. Here δµ = µ w − µ v denotes the difference of the chemical potentials of the water film and the ambient vapor phase, whereas E ev = ηl 2 0 Q e / h 3 0 is the evaporation number with effective rate constant Q e . The pressure in the vapor above the film is assumed to be close to the saturation pressure, allowing us to identify the chemical potential of the water film with the generalized pressure [30,31]; that is, µ w =P and µ v = const.
In the literature, different expressions for the disjoining pressure (H ) have been considered (see [21,32] and references therein for a discussion of possible choices). Here, we employ the dimensionless expression [33][34][35] (H, With this choice the substrate is always covered with a microscopic precursor film of height H p defined by (H p ) = 0 even if it is macroscopically dry. In the chosen scales H p = 1, so that the precursor height is the unit of film thickness.
Further, the value = 1 corresponds to the critical density of the monolayer, i.e. the density at the center of the coexistence region of the first-order LE-LC phase transition. The dependence of the surface tension on the surfactant density is determined by the equation of state where˜ := − 1 denotes the deviation of from the critical value = 1. Two coefficients, M 1 and M 2 , characterize the free energy of the surfactant layer that is modeled as a symmetric double-well centered around the critical density (see [20] for details). By choosing S(H ) = B (H ), where (H ) = dH (H ) with integration constant zero is the potential of the substrate-liquid interaction and B is a positive coupling constant, SMC acts on length scales comparable to the substrate-liquid interaction [19,20]. Note that the model given by equations (1) and (2) assumes that the monolayer covered surface is perfectly mobile; that is, the rigidity of the surfactant layer is not taken into account. It can, however, be estimated that the monolayer rigidity does not play a significant role for phospholipids such as DPPC in the vicinity of the LE-LC phase transition [20]. For other materials, the deformation of the surface layer may involve a significant amount of work. In that case, our model has to be extended following the approach outlined in [36].
The model given by equations (1) and (2) is solved using the boundary conditions where the constant boundary values H 0 and 0 correspond to the film height and the monolayer density of the reservoir. The latter is kept at a value corresponding to the LE phase. In the Y -direction, periodic boundary conditions are applied. The properties of the transferred monolayer depend on the transfer velocity V . For small velocities, the substrate is coated with a homogeneous LC monolayer. Increasing the velocity beyond a critical value V l , the homogeneous transfer becomes unstable and a periodic spatiotemporal pattern is created. This pattern consists of alternating domains of the LE and the LC phases (see figure 2). Beyond an upper velocity bound V u , the transfer becomes homogeneous again and the substrate is coated with a homogeneous LE monolayer. In the remainder of this paper, the range (V l , V u ) is referred to as the patterning range. Within the patterning range, two main transfer modes have to be distinguished: stripes perpendicular to the contact line are observed for lower velocities, while stripes parallel to the contact line are generated for higher V . When decreasing V from large values where the parallel stripes are deposited, they become unstable to transversal instability modes. Therefore, the perpendicular stripes may be seen as resulting from a secondary instability. In a time simulation, parallel stripes always emerge first, i.e. they occur at least as transient states that then transform into perpendicular stripes through an orientation transition [19,20]. Note that the perpendicular stripe pattern is stationary in the laboratory frame, while the parallel stripe pattern is time periodic and consists of domains that are perpetually generated and grow at the contact line, detach from it and finally are advected with the moving substrate.
Once the meniscus has taken a shape that depends on the transfer velocity, it remains almost static and exhibits only small oscillations with the formation of every new domain in the case of stripes parallel to the contact line [19,20]. It turns out that the meniscus oscillations are not a vital part of the mechanism of pattern formation but a mere consequence. The patterns are still observed in numerical simulations with an artificially frozen static meniscus. Therefore, the essential mechanism of the formation of the LE-LC stripes is the decomposition of the initially homogeneous layer due to SMC.
The described main findings of the full thin-film model motivate us to develop in the next section a reduced model that (i) still captures the main features of the patterning process and (ii) also allows for a profound analysis of its solution structure and therefore for a better understanding of the onset of patterning.

The reduced model
The fact that the meniscus is nearly static means that the interaction field S(H (X, T )) may be approximated by a static field S(X). From this perspective, the SMC phenomenon arises because the monolayer is pulled by the moving substrate through a spatially varying free energy landscape. Far to the left of the contact line, the film height is large; therefore the free energy is not affected by SMC and the LE phase is stable. As the monolayer approaches the contact line, the free energy is tilted towards higher densities: the minimum corresponding to the LC phase is lowered, so that condensation is energetically favored. The influence of the hydrodynamics is therefore restricted to the determination of the details of the spatial variation of S. As a result, one is able to capture the essence of the LE-LC domain formation under the influence of SMC by a variant of the prototypical model equation for phase decomposition, namely the Cahn-Hilliard equation, which is augmented by the addition of a space-dependent external field and a dragging term. Let a scalar field c(X, T ) evolve according to where V is a constant vector describing the velocity of substrate withdrawal. Due to the spacedependent SMC field, the free energy F(c, X) is a function of c and the spatial coordinates X.

It is given by
By inserting this expression into equation (5), the evolution equation for c finally becomes Note that in the present problem the advective term Vc cannot be removed from this equation by simply performing a Galilei transform because there exists one privileged frame of reference: the frame in which the contact line is at rest. For simplicity, we focus here on the case of a straight contact line that is chosen, without loss of generality, to be perpendicular to the X -axis, so that ζ (X) = ζ (X ) and V = V e X . The function ζ smoothly interpolates between the region of the symmetric free energy to the left of the contact line, where LE and LC phases are energetically equivalent, and the region of the tilted free energy to the right of the contact line, where the LC phase is energetically favored. That is, ζ models the field S that is determined in the full model (equations (1) and (2)) by the meniscus shape. It turns out that the results are not very sensitive to the particular form of ζ . Here, the smooth transition is modeled by a hyperbolic tangent centered at X = X s , where X s denotes the position of the meniscus.
and the width of the transition is determined by the constant l s . The boundary conditions are similar to those in the full model, i.e.
and periodic in the Y -direction. Although closely related to the problem considered here, the models investigated by Krekhov [25] as well as Foard and Wagner [26] differ from our approach as they use symmetric free energies that switch the sign of the c 2 term at the location of the moving front; that is, the free energies change from a single-well to a double-well potential. Furthermore, in the here considered monolayer transfer with SMC, the inflow of the surfactant is controlled by the transfer velocity V and is not an independent parameter. Equations (7)-(9) are solved numerically in one (1D) and two (2D) spatial dimensions on a domain of size [0, L] and [0, L] × [0, L], respectively. Finite differences are used for the spatial discretization, whereas an adaptive Runge-Kutta-Fehlberg scheme of orders 4 (5) is used for time stepping [37]. The 2D simulations are performed on an NVIDIA GeForce480 using the CUDA framework. For parallelized computation, the grid consisting of 384 × 384 points is decomposed into blocks of 16 × 16 points [38]. The numerical parameter values used in the calculations are: the smooth step function ζ is centered around X s = 10, the step width is l s = 2 and the strength of the free energy tilt is µ = 0.5. The boundary value is chosen to be c 0 = −0.9, i.e. slightly above c = −1, the equilibrium value in the absence of external fields, but still far from the spinodal value at c = −1/ √ 3 ≈ −0.5774. The numerical solutions show that the reduced system qualitatively behaves as the full model given by equations (1) and (2). If the transfer velocity V is below a certain critical value V l ≈ 0.0242, only homogeneous transfer of a layer with c ≈ 0.9 is observed, corresponding to the transfer of a homogeneous LE monolayer. Within a patterning range from V l to an upper critical velocity V u ≈ 0.0572, patterns consisting of alternating domains in the two phases  As for the full model (equations (1) and (2)) [19,20], the formation of the perpendicular stripes is always preceded in time by the formation of a number of parallel stripes which develop a transversal instability, eventually break up and transform into perpendicular stripes. This can be seen in figure 4 that gives snapshots at four different times from the transient stage of the transfer process at V = 0.027. This observation indicates that the onset of the formation of parallel stripes coincides with the onset of patterning in general.
The formation of the parallel stripes can be conveniently investigated by solving the model equations (7)-(9) in one spatial dimension. A plot of the 'wavenumber' k (obtained from the period P by k = 2π/P) of the transferred structures against the transfer velocity V is given in figure 5. The diagram shows that the branch of periodic solutions emerges very steeply from the V -axis at V l . As explained in section 4.2, the corresponding bifurcation is an infinite period bifurcation, so that the branch actually extends down to k = 0. At the upper critical velocity V = V u , the branch ends with a finite wavenumber. For comparison, see the analogous diagram for the full model given in figure 3 of [19].

Stationary solutions and bifurcation diagram
To complement the results of the direct numerical simulations, stationary solutions are calculated using the method of pseudo-arclength continuation [39]. The continuation algorithm is implemented in MATLAB and based on the same spatial discretization as the direct numerical simulations. A solution has to be supplied as a starting point of the continuation. For this purpose, a direct numerical simulation is performed for V = 0.02, i.e. at a velocity for which c relaxes to a stationary state. This state is then used as the first solution of the continuation process. The solution branch that passes through this solution is then traced in both directions. In a second step, a linear stability analysis of all stationary states obtained along the branch is performed. Writing equation (7) in the form with a nonlinear operator L, infinitesimal perturbations η = c −ĉ, whereĉ is a stationary solution, are governed by the linearized equation The stationary statesĉ(X ) are not homogeneous and therefore L c is space dependent. Moreover, L c is not self-adjoint so that its eigenfunctions and eigenvalues are in general complex. Here, they are calculated numerically using MATLAB. By this procedure, the spectrum of each + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + solution on the branch can be determined. Note that although the stability problem itself does not depend explicitly on the space-dependent external field ζ (X ), the function ζ (X ) determines the stationary statesĉ(X ) and therefore implicitly controls the linear stability.
We are now in a position to draw a bifurcation diagram combining the information gathered from the direct numerical simulations and the continuation. Because of the high dimensionality of the system the norm is used as a solution measure. In figure 6, it is plotted against the control parameter V . One finds that the branch of stationary solutions exhibits a large number of folds. However, only in the saddle-node bifurcations S 1 (V = 0.0541, L 2 = 1.133), S 2 (V = 0.0161, L 2 = 1.137) and S 3 (V = 0.0256, L 2 = 1.006) do stable and unstable states meet. At the other folds, more and more eigenvalues cross the imaginary axis, and on all involved branches there is at least one eigenvalue with a positive real part. Twelve profiles from different parts of the branch of stationary solutions are shown in figure 7. The uppermost stable part of the branch (left of S 1 ) consists of stable solutions of the type shown in figure 7(a). The profile c directly goes from its boundary value c = −0.9 to a higher value c ≈ 1.23. The variation of the external field strength around X s = 10 yields merely a small undulation of the profile. Along the unstable part of the branch between S 1 and S 2 (for a typical profile see figure 7(b)), the position of the step from low to high c is continuously shifted to X s = 10. On the next part (between S 2 and S 3 ), stable solutions of the type shown in figure 7(c) are found. In the context of the LB transfer, these solutions can be associated with the homogeneous transfer of a condensed layer. Following the stable part further, the step is shifted away from X s to the right until the solutions become unstable again in the saddlenode bifurcation S 3 . Note that there is exactly one eigenvalue with a positive real part along the unstable part of the branch between the two saddle-node bifurcations S 3 and S 4 (V = 0.00964 and L 2 = 0.924), so that the corresponding solutions on the branch are saddle points in the function space. As can be seen from the solutions presented in figures 7(d)-(i), beyond S 4 each pair of subsequent saddle-node bifurcations (folds) of the branch results in the addition of another bump to the structure that develops in the transition region between low and high concentration. The new bump is always formed at X s and then 'moves' to the right as the next bump forms. A related solution behavior is known from an investigation of other model equations of pattern formation, such as the Swift-Hohenberg equation [40]. There, localized structures with different numbers of bumps or holes may be observed for the same value of some control parameter, just as we observe different numbers of bumps for the same transfer velocity. Usually such localized structures are investigated for large domains with the same conditions at all boundaries or for large period periodic systems, so that all localized solutions are homoclines. Due to this homoclinic character and the snake-like appearance of the solution branches of localized states in the bifurcation diagrams the phenomenon has been termed homoclinic snaking [41][42][43]. Therefore, here it is tempting to use the term heteroclinic snaking to describe the solution structure found, as the localized trains of bumps connect regions of a low and a high value of c. However, even if the states on the boundaries may formally be seen as fixed points of the Cahn-Hilliard equations with different respective homogeneous external potentials, due to the spatial change of the external potential the localized solutions found do not represent heteroclinic orbits of any autonomous system of ordinary differential equations. They do, however, still represent fronts that mediate the transition between different states. As a consequence, the phenomenon observed here could be called front snaking. The snaking of the branch in the bifurcation diagram ends when the whole domain right of the change in the external field is filled with bump structures (cf figure 7(i)). The amplitude of the bumps becomes smaller and smaller as one follows the branch further until finally the profile becomes monotonic again. Increasing V further, the slope in c and its value at the right boundary decrease (see figures 7(j) and (k)). Finally, the branch stabilizes again when at (V = V u = 0.0572, L 2 = 0.607) a pair of complex conjugate eigenvalues crosses the imaginary axis. The branch remains linearly stable for higher velocities. Figure 7(l) shows a typical solution found for V > V u . In the context of LB transfer, these solutions can be associated with the homogeneous transfer of a surfactant layer in the LE phase.
An averaging of the L 2 -norms of the time-periodic solutions obtained from direct numerical simulations with respect to time allows us to include them into the bifurcation diagram in figure 6. The corresponding branch emerges at (V = 0.0242 and 0.980) from the unstable part between S 3 and S 4 of the branch of steady solutions. The bifurcation is not related to a change of linear stability of the steady branch. The branch of time-periodic solutions ends at the bifurcation point at (V = V u = 0.0572 and L 2 = 0.607). In the following two subsections the two bifurcations will be investigated more closely.  The type of bifurcation can be identified by determining the scaling behavior of the temporal period of the oscillations in the vicinity of the bifurcation point [44]. To this end, a large number of direct numerical simulations with V close to the critical velocity is performed. In the velocity range V l ≈ 0.0242 V 0.0256, the branch of the periodic solutions coexists with two other stable branches (see figure 6) and it depends on the initial conditions onto which branch a solution settles. Therefore, it is necessary to follow the branch of the periodic solutions very carefully. One has to start with a profile obtained in the region V > 0.0256 and then decrement the velocity in small steps in a series of simulation runs, whereupon each run uses the final state of the previous run as its initial condition. Close to the bifurcation point, the velocity decrement has to be chosen smaller and smaller, because even a tiny V results in a significant change of the temporal period. In this way, the critical velocity has been determined up to 12 significant digits: V l = 0.024 186 808 7724. Figure 9 shows a logarithmic plot of the temporal period T P of the solutions against the distance from the bifurcation point V − V l . Clearly, T P diverges logarithmically as V l is approached. This scaling behavior indicates a homoclinic bifurcation [45,46]. In such a bifurcation, a limit cycle (the time-periodic solution) collides with a saddle point (the stationary solution). This leads to the following scaling law of the temporal period:

The onset of pattern formation at low velocities
where λ u is the eigenvalue of the unstable direction of the saddle point. bifurcation leading to the onset of the formation of parallel stripes at low velocities is indeed identified as a homoclinic bifurcation.
Here it seems pertinent to point out that the observed transition at V l from stationary solutions to the deposition of stripes on the moving substrate can as well be seen as a depinning transition. However, the details of our transition are unlike the depinning via a homoclinic bifurcation described in [46] where a drop depins from a substrate heterogeneity under the influence of an external driving force, and the entire drop slides along the substrate. Here, after depinning the concentration profile at the static meniscus does not move in its entirety with the dragging substrate. Instead, only a part of it (a depression) starts to move, resulting in the deposition of a line. This process then repeats periodically. Note that in the context of drops, for other parameter values the depinning transition may also be due to a Hopf or a sniper bifurcation [47,48]. The former is related to what we find here at large V (see the next subsection).

Amplitude equation.
The branch of the periodic solutions related to the deposition of stripe patterns ends at (V = V u = 0.0572 and L 2 = 0.607). There, according to linear stability analysis, a pair of complex conjugate eigenvalues crosses the imaginary axis. Actually, this bifurcation is the last one of a series of consecutive Hopf bifurcations as illustrated in figure 10 where selected eigenvalue spectra for V V u are given. The occurrence of several of these bifurcations in the vicinity of V u explains the very long transients that are observed in time simulations close to V u : several modes decay very slowly around this point. Furthermore, the presence of so many eigenvalues close to the imaginary axis limits the applicability of an asymptotic theory to a very narrow range of velocities. Nevertheless, the bifurcation can be investigated perturbatively by means of standard tools: we first determine the amplitude equation for the critical mode. At the bifurcation point V = V u the corresponding eigenvalues are purely imaginary; that is, λ c = ± i ω c . Thus, close to the final bifurcation at V = V u , i.e. for with the perturbation η(X, T ) given by the expansion of the eigenfunctions whereĉ(X ) denotes the stationary solution at the bifurcation point and ξ u (T ) denotes the complex amplitude of the critical mode that becomes unstable at the bifurcation. Note that this approach is based on a separation of timescales; that is, T is a 'slow time', so that ∂ T ξ u ∼ ξ u (see chapter 2.2 of [49]). The complex eigenfunction corresponding to ξ u is denoted by φ(X ) and the imaginary part of the corresponding eigenvalue at the bifurcation point is the frequency ω c . The functions 0 (X, T ) and 2 (X, T ) describe the contributions of the stable modes of the spectrum due to the harmonics e 0iω c T and e 2iω c T . The evolution equation for the perturbation η can be expanded as where the operators L c , L c and L c denote the Frechet derivatives of the nonlinear operator defined in equation (10) evaluated atĉ. The operator L c depends on the control parameter V (see equation (11)). Close to the bifurcation point, it can be approximated by Now, we insert expression (12) into equation (13) and discuss the different orders in e iω c T separately. By using the eigenvalue equation L c φ = iω c φ at the bifurcation point, the resulting equation for O(e iω c T ) iṡ 16 where the asterisk denotes complex conjugation. For O(e 2iω c T ), the equation ∂ 2 ∂ T + 2iω c 2 = L c 2 + 1 2 ξ 2 u L c φφ is obtained. The time derivative of 2 can be eliminated adiabatically [50] and the amplitude 2 can be expressed in terms of ξ u as 2 = X 2 ξ 2 u , where X 2 denotes a function living in the same space as the eigenfunction φ and ξ 2 u is the corresponding amplitude. Then, one obtains the solvability condition for the function X 2 Again, the time derivative is eliminated and 0 is rewritten as 0 = X 0 |ξ u | 2 . Therefore, the term ∼ 2 0 is of order |ξ u | 4 and is negligible. This yields L c X 0 = −L c φφ * .
To obtain an evolution equation for ξ u , the state vectors have to be projected. As the operator L c is not self-adjoint, its eigenfunctions do not obey an orthogonality relation. Instead, the components of the modes have to be calculated by projection on the eigenfunctions of the adjoint eigenvalue problem L † c φ † = λ * φ † . Here, L † c is the adjoint operator of L c and its eigenvectors are denoted by φ † By the use of the scalar product φ|ψ = L 0 dX φ * (X )ψ(X ), equation (14) is then projected onto φ † , the mode of the adjoint operator L † belonging to the eigenvalue λ * = −iω c . This yieldṡ ξ u = a 1 ξ u + a 2 ξ u |ξ u | 2 (15) with the complex coefficients and Equation (15) is the lowest-order normal form of a Hopf bifurcation [49]. If Re a 1 < 0, the solution ξ u = 0 is linearly stable and perturbations are damped. For Re a 1 > 0, the solution ξ u = 0 is linearly unstable. Then, for Re a 2 < 0 the Hopf bifurcation is supercritical and a stable periodic solution exists (for Re a 1 > 0) that has the form ξ u = Re iω c T with the radius However, if Re a 2 > 0 the Hopf bifurcation is subcritical and to lowest order no stable periodic solution exists for Re a 1 > 0. Instead, an unstable periodic solution exists for Re a 1 < 0 [49,51].
In our case, the derivatives of the operator L acting on arbitrary states φ, ψ and χ are given by equation (11) and L c φψ = 6 ĉψ φ + ψφ ĉ +ĉφ ψ + 2ĉ (∇ψ) · (∇φ) + 2ψ ∇ĉ · (∇φ) + 2φ ∇ĉ · (∇ψ) , L c φψχ = 6 (χ ψ φ+χ φ ψ+ψφ χ +2χ (∇ψ) · (∇φ)+2ψ (∇χ ) · (∇φ) + 2φ (∇χ) · (∇ψ)) , The vectors φ, φ † , X 0 and X 2 , as well as the coefficients a 1 and a 2 , are calculated from the discretized versions of L c , L c and L c . This discretization is performed in the same way as in the direct numerical simulation of the problem, i.e. by finite differences. Inserting the results into equations (16) and (17), one obtains respectively. As Re a 1 < 0, the stationary solutions are unstable for = V − V u < 0, i.e. for velocities below V u . This is in agreement with our time simulations and the linear stability results. The finding that Re a 2 > 0 indicates that the final Hopf bifurcation is subcritical. Higher order asymptotics would be needed to obtain further information on the system behavior. However, as such an approach needed to take into account the influence of the other Hopf bifurcations occurring nearby, it would be rather cumbersome and is not pursued here (also see the discussion below).
In fact, the behavior of the solutions slightly below the bifurcation shows no trace of a stable limit cycle. Starting a simulation at V = 0.057 191 33, a velocity that lies only 2 × 10 −8 below V u , a long phase of slow exponential growth of the unstable mode is observed, corresponding to the very small real part of the eigenvalue close to the bifurcation. The unstable mode φ has the shape of an almost sinusoidal wave traveling to the right, whereas the amplitude grows slowly from left to right (see figure 11). The regime of linear growth ends abruptly and further modes are excited. In parallel, one observes the formation of larger bumps or even complete domains with c ≈ 1 close to the right boundary of the integration domain. These bumps leave the integration domain quickly at the right and the dynamics settles back to the linear behavior for some time. Sooner or later, however, the formation of a larger bump occurs again. By projection of the solution onto the unstable mode φ † , the complex amplitude ξ u can be calculated and a phase plot (Re ξ u , Im ξ u ) can be drawn, as shown in figure 12. The linear phases of the dynamics correspond to a slowly growing spiral. The formation of bumps can be seen from the excursion from these planar spirals that follow an obviously higher-dimensional trajectory that eventually comes back to the plane and begins to follow a spiral again. The absolute value of the amplitude, |ξ u |, is shown in figure 13. From these results the following picture emerges: for V < V u the stationary solutions become unstable due to a Hopf bifurcation. The related unstable linear mode results in a time evolution corresponding to an unstable spiral in the phase plane. In the nonlinear regime, the spiral motion is caught by a higher-dimensional structure due to the influence of the stable modes that are next to the imaginary axis. A possible scenario involving two pairs of complex conjugate eigenvalues λ 1,2 = −γ ± iω and λ 3,4 = λ ± iα with λ, γ , α, ω > 0 and λ = γ has been investigated by Ovsyannikov and Shilnikov (see case 3 in [52]). However, explicit calculation shows that it does not suffice to extend the above calculation by the addition of one pair of stable modes ξ s,1 , as the resulting dynamical system for the two amplitudes ξ u , ξ s,1 does not exhibit the characteristic dynamics shown in figure 12. As can be seen from the absolute values of the ξ s,i shown in figure 13, even the third stable pair of modes seems to influence significantly the dynamics during the large excursions that correspond to the domain formation. Thus, it can be concluded that even though the description of the behavior close to the bifurcation in terms of amplitude equations is possible, it will probably be too complicated to be of practical use.

Marginal stability analysis.
However, a change of perspective allows us to derive additional analytic results regarding the considered transition. By switching from the fixed laboratory frame to a frame comoving with the substrate, one can consider the location of the domain formation to be a propagating front which closely follows the contact line. This idea is particularly appealing, since it allows us to calculate the front velocity analytically. It also provides an intuitive explanation of the cessation of the patterning: when the pull velocity exceeds the front velocity, the meniscus escapes the front. However, the proposed picture also has some limitations. First of all, normally, analytical results for front propagation are based on the assumption that the front moves into a homogeneous unstable state [53]. This is not given in our case, since the front moves into the meniscus region which is determined by the boundary conditions at X = 0 and which is clearly non-homogeneous. Nevertheless, the stationary states that develop for velocities V ≈ V u are almost constant at the rear end, i.e. for X ≈ L. If this almost constant value corresponds to an unstable state, one can-at least approximately-calculate the velocity of a front that propagates into the integration domain from the right.
This calculation is based on the marginal-stability criterion [53]. Consider the leading edge of the front, where it deviates only slightly from the homogeneous state, so that it is well described by the linearized version of the evolution equations. The linear equations can be solved by a superposition of plane waves with wavenumber k and amplitudec(k) in the form where ω(k) is the dispersion relation of the linearized system. Assuming the existence of a finite front velocity V f , it is possible to switch to a comoving reference frame yielding whereω(k) = ω(k) − kV f is the dispersion relation in the comoving frame. If the exponent is sufficiently large, integrals like equation (18) can be evaluated by performing a contour integration in the complex plane, thereby using the saddle point approximation [53]. This is convenient in our case, as we are interested in the asymptotic behavior for T → ∞. If the contour is deformed so that it passes through the saddle point and does so in the direction of steepest descent, then the largest contributions of the integrand are located around the saddle point itself. Thus, the integral can be approximated by its value at k 0 .
Here, it is not necessary to carry out the integration. Instead it suffices to know that the integral will be dominated by the contribution around the saddle point. We therefore look for a wavenumber k 0 with ∂ω Since the leading edge of the front is stationary in the comoving frame, its envelope will be stationary, so that there is neither growth nor decay; that is, The conditions (19) and (20) allow us to calculate k 0 and V f . Linearization of equation (7) in a comoving frame yields where we use the shorthand notationF cc := (∂ 2 F/∂c 2 )| c=ĉ . Note that the inhomogeneity described by ζ does not appear in this equation, because we are considering the rear end of the integration domain where ζ can be assumed constant. For propagation into an unstable state, i.e. forF cc < 0, one obtains Here, only the absolute value of V f is important to answer the question of whether or not the front can keep up with the moving contact line. Now we consider the situation close to the final Hopf bifurcation at V u discussed in section 4.3.1. Away from the contact line region, the profiles obtained by continuation are almost completely flat (see figure 7(k)) and can thus be approximately regarded as homogeneous with c = c L = const, where c L denotes the value of c at the right boundary of the integration domain. Figure 14 shows the values c L of the stationary solutions obtained by continuation as a function of the pull velocity V f . With this information, it is possible to establish a direct relation between the pull velocity V and the front velocity V f in the form of the function V f (c L (V )). Of course, the front can only keep up with the contact line as long as V f > V . Figure 15 shows the plot of V f (c L (V )) together with the bisecting line V f = V . Both lines intersect at V = 0.0575, a value that differs by less than 0.6% from the exact location of the bifurcation as determined from the linear stability analysis of the stationary solutions. Given the approximations involved in the calculation outlined above, the achieved level of agreement is considerable and leads to the conclusion that it is indeed justified to conceive the pattern formation process at the moving contact line as pattern formation in the wake of a propagating front.

Conclusions and outlook
We have shown that the essential dynamics of the transfer of surfactant monolayers from a water sub-phase onto a solid substrate can be captured by a reduced model based on an amended Cahn-Hilliard equation. It incorporates the SMC effect-which is responsible for the patterning-via a spatially varying free energy and the moving substrate by an advective transport term. The analysis of the reduced model has shown that it qualitatively reproduces all the main features of the patterning process as obtained with the full thin-film model in [19,20]. Furthermore, the reduced model has allowed for a deeper analysis of its solution structure and therefore for a better understanding of the onset of patterning. In particular, the results of direct numerical simulations and continuation methods have been combined into a bifurcation diagram that has allowed us to discuss how the time-periodic solutions related to the formation of stripes parallel to the contact line emerge from the stationary solutions. The onset at low velocities has been identified as a homoclinic bifurcation based on an analysis of the characteristic scaling of the temporal period. This implies the existence of a hysteresis; that is, there exists a range of substrate velocities where one may find transferred monolayers that are either homogeneous or patterned, depending on initial conditions. Furthermore, it has been found that when at high velocities the stripe patterns cease to exist several Hopf bifurcations are involved. As the final one is a subcritical bifurcation a lowest-order perturbation approach is not sufficient to describe the system behavior close to the transition. As higher order asymptotics would need to take into account the influence of the other Hopf bifurcations occurring nearby, it has not been pursued here. However, using an alternative approach the upper critical transfer velocity has been calculated independently. Switching into the frame comoving with the substrate, we have employed the marginal stability criterion for fronts propagating into unstable states to calculate the maximal substrate velocity that allows for the formation of stripes in very good agreement with the linear stability analysis of the stationary states. Therefore, the formation of stripe patterns in the transfer of phospholipids close to a first order phase transition can be understood as pattern formation in the wake of a propagating front.
Although hardly accessible in experiments, the snaking of localized states near a front that is exhibited by the unstable steady solutions obtained for low velocities is highly interesting from the mathematical point of view and will be the subject of future investigations.
The clear theoretical picture of the onset of the stripe pattern formation at low velocities needs to be substantiated by experimental data. Careful measurement should allow us to resolve the hysteresis behavior resulting from the coexistence of stationary and time-periodic solutions at the same velocities (see figure 6). Such experiments are particularly challenging in the small velocity regime as there the parallel stripes are unstable with respect to transversal instability modes, resulting in a transition of the orientation from parallel to perpendicular stripes. Our 2D calculations show, for both the full and the reduced model, that within the low velocity range only a limited number of parallel stripes are obtained as a transient state. This raises the question of whether the transversal instability might be suppressed deliberately as would be desirable for applications. It has been shown that prepatterned substrates can be used in order to modify the orientation transition [20], but this approach alters the properties of the primary instability as well.
From these considerations it is clear that the orientation transition needs to be further investigated theoretically and experimentally in order to gain a complete understanding of the pattern formation during monolayer transfer. At low velocities fully time-periodic 2D solutions corresponding to the deposition of parallel stripe patterns could be constructed from their 1D counterparts that are available as numerical results. By analyzing the linear stability of these solutions with respect to transversal perturbations the stability of the parallel stripes may be determined. However, this is beyond the scope of this paper and will be the subject of future research.
A further promising task would be the application of the techniques we have employed here to a similar problem of pattern formation in the wake of a quench front as studied by Krekhov [25] and Foard and Wagner [26]. This would allow one to investigate the structure of the stationary solutions as well as the bifurcations leading to the formation of the various patterns observed by these authors.
The formation of stripe patterns at the contact line during the Langmuir-Blodgett transfer is prototypical for a variety of other physical systems that exhibit pattern formation at a receding contact line due to mechanisms that are at the moment not yet well understood. Thus our mathematical approach could be extended to models that describe the experimentally observed formation of regular microscopic stripe patterns of metallic monolayers, e.g. of gold and silver [54], the deposition of stripe-like patterns of polymers or micro-and nanoparticles from solution [27,[55][56][57] and microstripes of organic semiconductors [58]. In particular, we expect a similar bifurcational structure in a recent model of line deposition from evaporating films of colloidal suspensions and polymer solutions [59].
As already pointed out at the end of section 4.2, the system shows similarities to different classes of systems. One example is pinned drops on heterogeneous substrates that show depinning transitions via sniper, homoclinic and Hopf bifurcations [46][47][48]. We have argued that the onset of the deposition of stripes can also be regarded as a depinning transition. Another example is systems related to waterjet cutting processes using abrasive particles. There, in a certain parameter range, the surface of a solid substrate that is moving through the cutting jet is patterned periodically [60] and the pattern is formed behind a front that is localized at the position of the waterjet.
The ubiquity of forced patterning caused by imposed moving fronts and its importance for various technical applications indicates that our reduced model and the analysis we propose might constitute a useful stepping stone toward a more general understanding of pattern formation in such systems.