# Lamellae alignment by shear flow in a model of a diblock copolymer

## I Introduction

We study a mesoscopic model of a block copolymer to describe the re-orientation of a lamellar structure by an imposed uniform shear flow that is either constant or periodic in time. This is a first step towards understanding known phenomenology pertaining to the response of the block copolymer microstructure to shear flows near the isotropic to lamellar transition . The model that we use is based on the free energy of the diblock copolymer obtained by Leibler , and later by Ohta and Kawasaki , to which an advection term is added to incorporate the effect of the externally applied shear flow. We identify spatially periodic solutions that correspond to a lamellar structure, and determine their stability against a number of long wavelength perturbations.

Modulated phases are ubiquitous in physical and chemical systems . They generally result from the competition between short and long range forces. Additional symmetries of the system (e.g., translational or rotational invariance) often lead in practice to rich textures, especially in systems of extent that is large compared with the characteristic wavelength of the modulation. Modulated phases often have interesting macroscopic behavior, and exhibit a complex response to externally applied forces. While it is possible to devise approximate constitutive laws to describe the macroscopic response of such phases, it is often necessary to explicitly address their evolution at the mesoscopic scale, and to determine how microstructure evolution influences the macroscopic response.

We focus here on the lamellar phase observed in diblock copolymers below the order-disorder transition . Diblock copolymers are formed by two distinct sequences of monomers, A and B, that are mutually incompatible but chemically linked. At sufficiently low temperatures, species A and B would segregate to form macroscopic domains, but the chemical bonding between the two leads to a modulated phase instead. The detailed equilibrium microstructure depends on the relative molecular weight of the chains  and has been studied in detail within a mean field approximation .

We follow in this paper the approach of Leibler who introduced an order parameter field $`\psi (𝐫)`$ that describes the local number density difference of monomers A and B. The order parameter is defined to be zero above the order-disorder transition, and is finite and nonuniform below. Leibler’s analysis was restricted to the weak-segregation limit (close to the order-disorder transition) within which the thickness of the interface separating the A-rich from the A-poor regions is of the order of the wavelength of the microstructure. Later, Ohta and Kawasaki extended Leibler’s free energy to the strong segregation range, and showed the importance of long ranged effective interactions that arise from the connectivity of the polymer chains. We use this latter free energy as the driving force for the re-orientation dynamics, allowing also for passive advection of the order parameter by an imposed shear flow. The model studied is similar to that considered by Fredrickson , except that we neglect thermal fluctuations and assume that both phases have the same viscosity.

The stability of a lamellar structure to secondary instabilities has already been addressed in the literature, although in the absence of shear flow . In fact, the similarity between the equations governing the motion of the lamellae and the Swift-Hohenberg model of Rayleigh-Bénard convection  gives rise to a common phenomenology . The lamellar structure is found to be stable only within a range of wavenumbers. At higher wavenumbers it undergoes an Eckhaus instability which generally results in a decrease of wavenumber, whereas for wavenumbers below that range the structure undergoes a zig-zag instability. In this paper, we extend these stability results to explicitly include fluid advection by the imposed shear. We find that the stability boundaries are modified with respect to the zero velocity case in a way that depends not only on the amplitude of the shear, but also on the orientation of the lamellae relative to the flow. Of course, the latter dependence is absent in earlier treatments that neglected advection.

Our results are a first step towards understanding the complex re-orientation phenomenology that has been observed experimentally . In this initial analysis, we introduce a number of restrictive assumptions that we plan to relax in future work. First, our calculations are primarily two dimensional and thus can only address the so-called parallel and transverse orientations. Second, and more importantly, we neglect thermal fluctuations and any viscosity contrast between the two phases, elements that have been argued to be important in determining the main qualitative features of the re-orientation process. We also neglect flow induced by the lamellae themselves in response to the applied shear. These secondary flows could become important for the late stage coarsening of the lamellar structure. Finally, we have confined our study to locating the boundaries of several secondary instabilities of the lamellar structure, but have not addressed the evolution following the instabilities, nor the coarsening of the resulting textured pattern . However, our results concerning the periodic base solution under flow, and its stability against long wavelength perturbations are the prerequisite building blocks of a more general theory.

## II Mesoscopic model equations

Following Leibler , we introduce an order parameter field, $`\psi (𝐫)`$, function of the local density difference of monomers A and B. For a block copolymer with equal length sub-chains, the order parameter is $`\psi (𝐫)=\frac{\rho _A(𝐫)\rho _B(𝐫)}{2\rho _0}`$, where $`\rho _X,X=A,B`$ is the density of monomer $`X`$, and $`\rho _0`$ is the total density, assumed constant (incompressibility condition). A mean field free energy $`\left[\psi (𝐫)\right]`$ was derived by Leibler for a monodisperse diblock copolymer melt , and later by Ohta and Kawasaki . In units of $`k_BT`$, where $`k_B`$ is Boltzmann’s constant and $`T`$ the temperature, the free energy is comprised of two terms, $`/(\rho _0k_BT)=_s+_l`$. The term $`_s`$ incorporates local monomer interactions,

$$_s=\text{}\text{r}\left[\frac{\kappa }{2}|\psi |^2\frac{\tau }{2}\psi ^2+\frac{u}{4}\psi ^4\right],$$

and is formally identical to the Ginzburg-Landau free energy commonly used to describe phase separation in a binary fluid mixture . The parameters $`\kappa ,\tau `$ and $`B`$ can be approximately related to the polymerization index $`N`$, Kuhn’s statistical length $`b`$ and the Flory-Huggins parameter $`\chi `$ through the relations $`\kappa =\frac{b^2}{3}`$, $`\tau =\frac{2\chi N7.2}{N}`$ and $`B=\frac{144}{N^2b^2}`$ .

Long range interactions arising from the covalent bond connecting the two sub-chains are contained in $`_l`$,

$$_l=\frac{B}{2}d𝐫d𝐫^{}G(𝐫r^{})\psi (𝐫)\psi (𝐫^{})$$

where the kernel $`G(𝐫r^{})`$ is the infinite space Green’s function of the Laplacian operator $`^2G(𝐫𝐫^{})=\delta (𝐫𝐫^{})`$. The nonlocal interactions arising from the connectivity of the chains lead to a thermodynamic equilibrium state with a nonuniform density. In our case of equal length sub-chains, the equilibrium configuration is a periodic lamellar structure, with a characteristic wavelength of the order of 100 Å for a typical system.

Given this mean field free energy, a phenomenological set of equations that govern the temporal relaxation of equilibrium thermal fluctuations of $`\psi (𝐫)`$ and of fluid velocity $`𝐯`$ has been derived close to the order-disorder transition . A similar phenomenological description can be used below the order-disorder transition under the assumption that the local relaxation of the order parameter field at the mesoscopic scale is still driven by minimization of the same free energy . Under this assumption, $`\psi `$ obeys the time-dependent Ginzburg-Landau equation,

$$\frac{\psi }{t}+𝐯\psi =M^2\frac{\delta }{\delta \psi },$$
(1)   
where $`M`$ is a phenomenological mobility coefficient, and $`\delta /\delta \psi `$ stands for functional differentiation with respect to $`\psi `$. Equation (1) includes the effect of advection by a local velocity field v, which satisfies an extended Navier-Stokes equation

$$\frac{𝐯}{t}+(𝐯)𝐯=\nu ^2𝐯\frac{p}{\rho }+\frac{\delta }{\delta \psi }\frac{\psi }{\rho },$$
(2)   
where $`\nu `$ is the kinematic viscosity of the fluid, assumed constant and independent of $`\psi `$, $`p`$ is the fluid pressure, and appropriate boundary conditions for both $`\psi `$ and v must be introduced. The last term on the right-hand side of Eq. (2) is required to ensure that there cannot be free energy reduction by pure advection of $`\psi `$ . This term is sometimes referred to as osmotic stress, and it leads to the creation of rotational flow by curved lamellae that is directed towards their local center of curvature.

We focus on a layer of block copolymer, unbounded in the $`x`$ and $`y`$ directions, and being uniformly sheared along the $`z`$ direction (Fig. 1). The layer is confined between the stationary $`z=0`$ plane, and the plane $`z=d`$ which is uniformly displaced parallel to itself with a velocity $`v_{\mathrm{plane}}=sd`$ in the case of a steady shear, and $`v_{\mathrm{plane}}=\gamma d\omega \mathrm{cos}(\omega t)`$ in the case of an oscillatory shear. $`s`$ is the dimensional shear rate in the steady case, and $`\gamma `$ is the dimensionless strain amplitude in the case of an oscillatory shear of angular frequency $`\omega `$.

The general problem defined by Eqs. (1) and (2) can be considerably simplified by noting that under typical experimental conditions inertia is negligible ($`\omega d^2/\nu 1`$). Furthermore, we will neglect in this paper the term $`(\delta /\delta \psi )\psi /\rho `$ in Eq. (2). Under these conditions, Eq. (2) admits a simple solution that satisfies the specified conditions at the moving plane: $`𝐯=sz\widehat{i}`$ for a steady shear, and $`𝐯=\gamma d\omega \mathrm{cos}(\omega t)z\widehat{i}`$ for an oscillatory shear, where $`\widehat{i}`$ is the unit vector in the $`x`$ direction. Therefore, the problem reduces to a single governing equation for the order parameter field $`\psi `$ (Eq. (1)) under a prescribed advection velocity $`𝐯`$. As discussed in the introduction, previous theoretical work on the formation and stability of lamellar structures further neglected advection of $`\psi `$ in Eq. (1). The results presented in this paper are free of this restriction.

Since the base state to be considered is comprised of spatially uniform lamellae advected by the shear flow, it is convenient to introduce a new frame of reference in which the velocity vanishes. Define a new system of non mutually orthogonal coordinates $`(x_1,x_2,x_3)`$ by $`x_1=xa(t)z,x_2=y`$ and $`x_3=z`$. The dimensionless quantity $`a(t)=st`$ for a steady shear, and $`a(t)=\gamma \mathrm{sin}(\omega t)`$ for an oscillatory shear. All the calculations reported in this paper, both analytical and numerical, have been performed in this new frame of reference. Analytical calculations consider an unbounded geometry in the $`x_1`$ and $`x_2`$ directions and periodic boundary conditions along the $`x_3`$ direction, whereas the numerical computations have been conducted in a two dimensional, square domain on the $`(x_1,x_3)`$ plane and consider periodic boundary conditions along both $`x_1`$ and $`x_3`$. Note that both frame of references coincide at $`t=0`$, and at equal successive intervals of one half the period of the shear in the case of oscillatory shear.

Dimensionless variables are introduced by defining a scale of length by $`\sqrt{\kappa /\tau }`$, a scale of time by $`\kappa /M\tau ^2`$, and an order parameter scale by $`\sqrt{\tau /u}`$. In the transformed frame of reference and in dimensionless variables, Eq. (1) reads,

$$\frac{\psi }{t}=^2(\psi +\psi ^3^2\psi )\frac{B\kappa }{\tau ^2}\psi ,$$
(3)   
with

$$^2=\left[1+a^2(t)\right]\frac{^2}{x_1^2}2a(t)\frac{^2}{x_1x_3}+\frac{^2}{x_3^2}+\frac{^2}{x_2^2}.$$

There is only one dimensionless group remaining $`B\kappa /\tau ^2`$, which will be simply denoted by $`B`$ in what follows.

We will first show in Sec. III that below (but close to) the order-disorder transition point (in the weak-segregation limit), Eq. (3) admits periodic solutions. Their stability against infinitesimal long wavelength perturbations is the subject of Sec. IV.

## III Lamellar solution in the weak-segregation limit

In the absence of shear ($`a(t)=0`$) the uniform solution of Eq. (3), $`\psi =0`$, loses stability at the order-disorder transition. In a mean field approximation the transition occurs at $`B_c=1/4`$. This is a supercritical bifurcation with a critical wavenumber $`q_c=\sqrt{1/2}`$. Near threshold $`(0ϵ=(B_cB)/2B_c1)`$ there exist periodic stationary solutions of the form

$$\psi (𝐫)=2A\mathrm{cos}(𝐪r)+A_1\mathrm{cos}(3𝐪r)+\mathrm{},$$
(4)   
with $`A^2=\frac{q^2q^4B}{3q^2}O(ϵ)`$, and $`A_1`$ of higher order in $`ϵ`$. This solution only exists for a range of wavenumbers $`q`$ such that $`\sigma (q^2)=q^2q^4B0`$.

For nonzero shear, we seek solutions of Eq. (3) of the form of Eq. (4), with $`𝐫=(x_1,x_2,x_3)`$ expressed in the sheared frame basis set $`\{𝐞_1=\widehat{i},𝐞_2=\widehat{j},𝐞_3=a(t)\widehat{i}+\widehat{k}\}`$. Wavevectors are expressed in the reciprocal basis set $`\{𝐠_1=\widehat{i}a(t)\widehat{k},𝐠_2=\widehat{j},𝐠_3=\widehat{k}\}`$. Therefore we keep the same functional form as for nonzero shear, but allow a time-dependent amplitude $`A(t)`$. Note that the components of the wavevector $`𝐪`$ are assumed to be independent of time and given by $`q_1=q_x(t=0),q_2=q_y(t=0)`$ and $`q_3=q_z(t=0)`$ respectively. The wavevector itself depends on time through the time dependence of the reciprocal basis set. Such a solution corresponds to a spatially uniform lamellar structure with a time-dependent wavevector that adiabatically follows the imposed shear in the laboratory frame (see Fig. 1). Inserting Eq. (4) into Eq. (3), we find to order $`ϵ^{3/2}`$ ($`\sigma `$ is itself of order $`ϵ`$),

$$\frac{dA}{dt}=\sigma \left[q^2(t)\right]A3q^2(t)A^3,$$
(5)   
with $`q^2(t)=q_1^2+\left(a(t)q_1q_3\right)^2+q_2^2`$, and $`\sigma (q^2)=q^2q^4B`$. This nonlinear equation with time-dependent coefficients can be solved exactly in the two cases of steady and oscillatory shear flow.

In the case of a steady shear $`a(t)=st`$. We find,

$$A(t)=\left\{\frac{e^{2H(t)}}{A(0)^2}+6e^{2H(t)}_0^t𝑑t^{}e^{2H(t^{})}q^2(t^{})\right\}^{1/2},$$
(6)   
with

$$H(t)=(q_0^4+Bq_0^2)t+(12q_0^2)sq_1q_3t^2+\frac{2q_1^2s^2(2q_3^2+q_0^2)}{3}t^3q_1^3q_3s^3t^4+\frac{q_1^4s^4}{5}t^5.$$
(7)   
The constant quantity $`q_0=\sqrt{q_1^2+q_2^2+q_3^2}`$ is the initial wavenumber, and $`A(0)`$ is the initial amplitude. For the special case $`q_1=0`$, $`A(t)`$ simply relaxes to its equilibrium value in the absence of shear $`A^2=\frac{q_0^2q_0^4B}{3q_0^2}`$. This corresponds to an initial orientation of the structure which has no component transverse to the flow. For any other initial orientation, the shear induces changes in the lamellar spacing in the laboratory frame of reference (Fig. 1). As a result, the amplitude $`A(t)`$ decreases and approaches zero at long times. Hence, the structure melts and reforms with a different orientation which we cannot predict on the basis of our single mode analysis. The emerging structure presumably results from the amplification of thermal fluctuations near the point at which the amplitude $`A(t)`$ vanishes, and they have been neglected in our treatment. Thermal fluctuation effects have been accounted for by others .

For an oscillatory shear $`a(t)=\gamma \mathrm{sin}(\omega t)`$. We first examine the stability of the uniform solution $`\psi =0`$ against small perturbations. Linearization of Eq. (5) leads to,

$$\frac{dA(t)}{dt}=\sigma \left[q^2(t)\right]A(t),$$
(8)   
with $`\sigma (t+T)=\sigma (t)`$ and $`T=2\pi /\omega `$. Equation (8) constitutes a one-dimensional Floquet problem. The solution $`A=0`$ is unstable when

$$\overline{\sigma }=_0^T\sigma (t)𝑑t>0.$$
(9)   
The resulting neutral stability curve is given by,

$$B=q_0^2q_0^4\frac{3q_1^4\gamma ^4}{8}\frac{(2q_0^2+4q_3^21)\gamma ^2q_1^2}{2}.$$
(10)   
Instability modes can be conveniently classified by considering the relative orientation of the lamellae at $`t=0`$ and the shear direction. We define a parallel orientation, $`q_30,q_1=q_2=0`$; a perpendicular orientation, $`q_20,q_1=q_3=0`$; and a transverse orientation, $`q_10,q_2=q_3=0`$. The following instability points are identified depending on the orientation of the critical wavevector: a transverse mode with

$$B_c=\frac{1}{2}\frac{(2+\gamma ^2)^2}{8+8\gamma ^2+3\gamma ^4},q_{1c}=\sqrt{\frac{4+2\gamma ^2}{8+8\gamma ^2+3\gamma ^4}},$$
(11)   
a mixed parallel-perpendicular mode with

$$B_c=\frac{1}{4},q_{1c}=0,2q_{2c}^2+2q_{3c}^2=1,$$
(12)   
and a mixed parallel-transverse mode defined by

$$B_c=\frac{1}{4}\frac{7\gamma ^2+16}{15\gamma ^2+16},q_{2c}=0,q_{1c}=2\sqrt{\frac{1}{15\gamma ^2+16}},q_{3c}=\sqrt{\frac{3\gamma ^2+8}{30\gamma ^2+32}}.$$
(13)   
Note that the threshold corresponding to perturbations of wavevectors that do not have a projection along the transverse direction are not affected by the shear. Furthermore, neither the stability boundaries nor the values of the critical wavenumbers depend on the angular frequency $`\omega `$.

In what follows, we consider mainly two dimensional solutions in the plane $`q_2=0`$ (transverse and parallel orientations) to make contact with two dimensional numerical calculations. As an example, Fig. 2 shows the neutral stability curve in the $`(q_1,q_3)`$ plane for mixed parallel-transverse modes at $`ϵ=0.04`$, and for several values of the dimensionless strain amplitude $`\gamma `$. Recall that $`q_1=q_x(t=0)`$ and $`q_3=q_z(t=0)`$ define the initial orientation of the lamellae. The figure shows that the shear does not modify the neutral stability curve in the vicinity of $`q_1=0`$ (parallel orientation), whereas the curve is shifted near $`q_3=0`$ (transverse orientation). Large changes are observed for oblique wavevectors, including the complete suppression of the instability at sufficiently large values of the strain amplitude.

Above threshold, Eq. (5) can be solved to yield the time-dependent amplitude of the lamellar structure under oscillatory shear. We find,

$$A(t)=\left\{\frac{e^{2(I(t)c_4c_5)}}{A(0)^2}+6e^{2I(t)}_0^t𝑑t^{}e^{2I(t^{})}q^2(t^{})\right\}^{1/2}.$$
(14)   
The function $`I(t)`$ is given by

$$I(t)=c_1t+c_2\mathrm{sin}(2\omega t)+c_3\mathrm{sin}(4\omega t)+c_4\mathrm{cos}(\omega t)+c_5\mathrm{cos}^3(\omega t),$$
(15)   
with

* $`c_1=\left[\frac{3q_1^4\gamma ^4}{8}+\frac{(2q_0^2+4q_3^21)\gamma ^2q_1^2}{2}+q_0^4+Bq_0^2\right]`$
* $`c_2=\left[\frac{q_1^4\gamma ^4+(2q_0^2+4q_3^21)\gamma ^2q_1^2}{4\omega }\right]`$, $`c_3=\frac{q_1^4\gamma ^4}{32\omega }`$
* $`c_4=\left[\frac{(4q_0^22)\gamma q_1q_3+4\gamma ^3q_1^3q_3}{\omega }\right]`$, and $`c_5=\frac{4\gamma ^3q_1^3q_3}{3\omega }`$.

We note that the stability condition Eq. (10) is equivalent to $`c_1=0`$. Hence, the asymptotic behavior of $`A(t)`$ at long times changes qualitatively depending on the sign of $`c_1`$. For $`c_1>0`$, $`lim_t\mathrm{}e^{2I}=0`$, so that the integral in Eq. (14) tends to a finite constant. Since the prefactor $`e^{2I}`$ diverges exponentially, $`A(t)`$ decays to zero. If, on the other hand, $`c_1<0`$, $`A(t)`$ becomes periodic at long times. To prove this statement, we first rewrite the second term inside the curly brackets in Eq. (14) as

$$(t)=6e^{2f(t)}_0^t𝑑t^{}e^{2c_1(t^{}t)2f(t^{})}q^2(t^{}),$$
(16)   
with $`f(t^{})=c_2\mathrm{sin}(2\omega t^{})+c_3\mathrm{sin}(4\omega t^{})+c_4\mathrm{cos}(\omega t^{})+c_5\mathrm{cos}^3(\omega t^{})`$. Since both $`f(t^{})`$ and $`q^2(t^{})`$ are periodic with period $`T=2\pi /\omega `$, we can decompose Eq. (16) into

$$(t)=6e^{2f(t)}\left[\underset{j=1}{\overset{n}{}}e^{2jc_1T}_0^T𝑑t^{}e^{2c_1t^{}2f(t^{}+t)}q^2(t^{}+t)+e^{2c_1t}_0^{tnT}𝑑t^{}e^{2c_1t^{}2f(t^{})}q^2(t^{})\right],$$
(17)   
where $`n`$ is an integer such that $`0<tnT<T`$. In the limit of large $`t`$ and with $`c_1`$ negative, the last term on the right-hand side of Eq. (17) vanishes while the sum $`_{j=1}^ne^{2jc_1T}`$ converges to $`1/(e^{2c_1T}1)`$. Combining Eqs. (14) and (17) yields an asymptotically periodic solution for $`A(t)`$,

$$A(t)=\left[\frac{6e^{2f(t)}}{e^{2c_1T}1}_0^T𝑑t^{}e^{2c_1t^{}2f(t^{}+t)}q^2(t^{}+t)\right]^{1/2}.$$
(18)   
The condition $`c_1=0`$ can also be understood in terms of a critical strain amplitude $`\gamma _c`$ above which an existing lamellar structure of a given orientation at $`t=0`$ will melt (i.e., $`A(t)`$ will decay to zero at long times) The value of $`\gamma _c`$ that corresponds to $`c_1=0`$ is given by,

$$\gamma _c=[(b+\sqrt{b^24dc})/2d]^{1/2},$$
(19)   
with $`b=(2q_0^2+4q_3^21)q_1^2/2,c=q_0^4+Bq_0^2`$, and $`d=3q_1^4/8`$. Note again that the critical strain amplitude is independent of the angular frequency $`\omega `$.

In order to test the approximations involved in Eq. (4), namely that the wavevector $`𝐪`$ adiabatically follows the flow, and the single mode truncation for small $`ϵ`$, we have undertaken a numerical solution of the model equation in a two dimensional, square geometry (see the Appendix for the details of the numerical method). As a first example, we consider an oscillatory shear of angular frequency $`\omega =0.02`$ imposed on a lamellar structure of initial wavevector $`(q_1,q_3)=(0.687,0.098)`$. The critical strain amplitude for this initial orientation is $`\gamma _c=0.695`$. Figure 3 shows the temporal evolution of $`A(t)`$ for two values of $`\gamma `$, one larger and one smaller than $`\gamma _c`$. The solid lines are the predictions of Eq. (14), and the symbols are the results of the numerical calculation. The agreement in both cases is excellent.

## IV Secondary instabilities of the lamellar pattern

In order to address the stability of the lamellar pattern, we next consider long wavelength perturbations of the base state with wavevector $`𝐐=(Q_1,Q_2,Q_3)`$ such that its components are also constant in the sheared frame of reference. Close to threshold, perturbations evolve in a slow time scale compared to the inverse frequency of the shear. We therefore assume that the wavenumber of any long wave perturbation would adiabatically follow the imposed flow. Specifically, we consider a solution of the form,

$$\psi (𝐫,t)=[A(t)+\delta A_+e^{i𝐐r}+\delta A_{}e^{i𝐐r}]e^{i𝐪r}+\mathrm{c}.\mathrm{c}.$$
(20)   
where $`A(t)`$ is the nonlinear solution obtained in Section III. Substituting Eq. (20) into Eq. (3), and linearizing with respect to the amplitudes $`\delta A_+`$ and $`\delta A_{}`$, we find,

$$\frac{}{t}\left[\begin{array}{c}\delta A_+\\ \delta A_{}\end{array}\right]=L(t)\left[\begin{array}{c}\delta A_+\\ \delta A_{}\end{array}\right],$$
(21)   
with

$$L(t)=\left[\begin{array}{cc}l_+l_+^2B+6A(t)^2l_+& 3A(t)^2l_+\\ 3A(t)^2l_{}& l_{}l_{}^2B+6A(t)^2l_{}\end{array}\right],$$

and $`l_\pm =(𝐪\pm 𝐐)^2=(1+a(t)^2)(q_1\pm Q_1)^2+2a(t)(q_1\pm Q_1)(q_3\pm Q_3)(q_2\pm Q_2)^2(q_3\pm Q_3)^2`$. In general, the matrix elements $`L_{ij}`$ are complicated functions of time, and we have not attempted to solve Eq. (21) analytically. For $`\gamma <\gamma _c`$ the operator $`L`$ contains terms that are both periodic in time and decaying transients. At long enough times, $`A(t)`$ is given by Eq. (18), and the linear system Eq. (21) has periodic coefficients. Hence, it reduces to a two dimensional Floquet problem for the amplitudes $`\delta A_+`$ and $`\delta A_{}`$ .

In order to gain some insight into the stability problem, we first briefly review the known results for zero shear . In this case $`A(t)`$ is a constant, and the matrix elements of $`L(t)`$ are independent of time. An eigenvalue problem results by considering solutions of Eq. (21) of the form $`\delta A_\pm e^{\sigma _+t}+e^{\sigma _{}t}`$, and instability follows when either eigenvalue is positive. Two modes of instability are obtained: a zig-zag (ZZ) mode that leads to a transverse modulation of the lamellae ($`𝐐q=0`$), and an Eckhaus (E) mode that is purely longitudinal in nature $`𝐐q=Qq`$. In the zig-zag case, $`\sigma _+(𝐐)`$ has a maximum at

$$Q_{max,ZZ}^2=\frac{12q^23A^2}{2}.$$
(22)   
The eigenvalue $`\sigma _+(𝐐)`$ changes sign on the line $`q=q_c`$, which therefore defines the zigzag stability boundary. In the Eckhaus case, we find after some straightforward algebra that the perturbation with the largest growth rate is

$$Q_{max,E}^2=\frac{64\delta q^4(ϵ4\delta q^2)^2}{64\delta q^2},$$
(23)   
with $`\delta q=qq_c`$. Therefore the Eckhaus stability boundary is given by $`ϵ=12\delta q^2`$. These results are schematically summarized in Fig. 4. The hatched area is the region of stability of a lamellar solution in the absence of shear flow. It is worth pointing out that this stability diagram is identical to that of the Swift-Hohenberg model of Rayleigh-Bénard convection . Shiwa  has recently shown that in the weak-segregation limit ($`ϵ1`$), and in the absence of shear flow, the amplitude equation describing slow modulations of a lamellar solution is the same as the amplitude equation of the Swift-Hohenberg model near onset of convection. The same stability diagram has been derived by Kodama and Doi  by examining free energy changes upon distortion of a lamellar pattern.

We now return to the Floquet problem of Eq. (21) when $`A(t)`$ is a periodic function of time (Eq. (18)). Since $`A(t+T)=A(t)`$ ($`T=2\pi /\omega `$), the solution of (21) is given by,

$$\left[\begin{array}{c}\delta A_+\\ \delta A_{}\end{array}\right]=e^{\sigma t}\left[\begin{array}{c}\varphi _+(t)\\ \varphi _{}(t)\end{array}\right],$$
(24)   
with $`\varphi _\pm (t+T)=\varphi _\pm (t)`$. Equation (21) is then transformed to an eigenvalue problem within $`(0,T)`$,

$$\frac{}{t}\left[\begin{array}{c}\varphi _+(t)\\ \varphi _{}(t)\end{array}\right]=\sigma \left[\begin{array}{c}\varphi _+(t)\\ \varphi _{}(t)\end{array}\right]+L(t)\left[\begin{array}{c}\varphi _+(t)\\ \varphi _{}(t)\end{array}\right].$$
(25)   
Given that the function $`A(t)`$ is quite complicated, we have solved this eigenvalue problem numerically. The eigenvalue $`\sigma `$ can depend in principle on the wavevector of the base state $`𝐪`$, on the wavevector of the perturbation $`𝐐`$, and on the amplitude $`\gamma `$ and frequency $`\omega `$ of the shear. For ease of presentation, we have focused on the case $`ϵ=0.04`$ although extension to other values of $`ϵ`$ is straightforward.

Figures 5 and 6 summarize our results for the cases $`\gamma =0.2`$ and $`\gamma =0.4`$ respectively, and show the stability boundaries in the plane $`(q_1,q_3)`$, as well as the neutral stability curve already shown in Fig. 2. As before, $`(q_1,q_3)`$ is the wavevector of the lamellar structure at $`t=0`$. At fixed $`ϵ,\gamma `$ and $`\omega `$, these curves have been obtained by determining the loci of $`𝐪`$ at which the function $`\sigma (𝐐)`$ changes from a maximum to a saddle point at $`𝐐=0`$. First we note that any orientation of the lamellar pattern that is not initially close to either parallel or transverse is unstable to moderate shears. Second, and up to our numerical accuracy, these curves appear to be independent of angular frequency. Finally, and contrary to the case of no shear, the reciprocal basis vectors are not time independent. Since the components of both $`𝐪`$ and $`𝐐`$ are independent of time in the sheared frame, their mutual angle is not (except for the case in which they are parallel). However, the following statements can be made about the type of secondary instability. We have found that the secondary instability is of the longitudinal type only when either $`q_1=0`$ or $`q_3=0`$ (intersections between the lines marked with circles and the axes in Figs. 5 and 6). Otherwise, the angle between $`𝐪`$ and $`𝐐`$ is time-dependent. The lines on both figures marked with squares have the property that even though both $`𝐪`$ and $`𝐐`$ are functions of time, their angle oscillates periodically around $`90^o`$.

The cases discussed up to now concern long wavelength instabilities of the base periodic pattern that are associated with the broken translational symmetry of the original system by the appearance of a periodic pattern. We now show that it is possible to obtain analytical expressions for the stability boundaries against finite wavelength perturbations that may have some experimental relevance as well. In some experimental protocols, the lamellar pattern is first obtained in the absence of shear. The resulting configuration comprises regions or domains of locally parallel lamellae but with a continuous distribution of orientations. A shear flow is then initiated and the reorientation of the pattern studied as a function of time. The pattern obtained in the absence of shear may be now unstable to several finite wavenumber perturbations that would not have been observable in the case in which flow is present throughout the ordering process. In the latter case the unstable orientations would have decayed away during the process of formation of the lamellae. In addition, the approximation that we derive below is generally valid when $`Q_3`$ cannot approach zero, as is the case in a system of finite extent in the direction of the velocity gradient.

We first define the following linear transformation,

$$\left[\begin{array}{c}\delta _+\\ \delta _{}\end{array}\right]=\left[\begin{array}{cc}\frac{3A(t)^2l_{}}{\sigma _+\sigma _{}}& \frac{\sigma _+\sigma _{}+L_{22}L_{11}}{2(\sigma _+\sigma _{})}\\ \frac{3A(t)^2l_{}}{\sigma _+\sigma _{}}& \frac{\sigma _+\sigma _{}+L_{11}L_{22}}{2(\sigma _+\sigma _{})}\end{array}\right]\left[\begin{array}{c}\delta A_+\\ \delta A_{}\end{array}\right]=T(t)\left[\begin{array}{c}\delta A_+\\ \delta A_{}\end{array}\right]$$
(26)   
which diagonalizes matrix $`L(t)`$, and where

$$\sigma _\pm (t)=\frac{L_{11}+L_{22}\pm \sqrt{(L_{11}L_{22})^2+36l_+l_{}A(t)^4}}{2}.$$
(27)   
Combining Eqs. (21) and (26), we find

$$\frac{}{t}\left[\begin{array}{c}\delta _+(t)\\ \delta _{}(t)\end{array}\right]=\left[\begin{array}{cc}\sigma _+& 0\\ 0& \sigma _{}\end{array}\right]\left[\begin{array}{c}\delta _+\\ \delta _{}\end{array}\right]\frac{A(t)^2l_{}}{\sigma _+\sigma _{}}\left[\begin{array}{cc}\dot{M_1}& \dot{M_2}\\ \dot{M_1}& \dot{M_2}\end{array}\right]\left[\begin{array}{c}\delta _+\\ \delta _{}\end{array}\right],$$
(28)   
where $`\dot{M_1}=\frac{}{t}\left[\frac{\sigma _+\sigma _{}+L_{11}L_{22}}{2A^2l_{}}\right]`$ and $`\dot{M_2}=\frac{}{t}\left[\frac{\sigma _{}\sigma _++L_{11}L_{22}}{2A^2l_{}}\right]`$. For finite $`Q`$, $`\frac{36A(t)^4l_+l_{}}{(L_{11}L_{22})^2}1`$. Assuming $`L_{11}L_{22}>0`$ (the other case leads to no extra complications), $`\sigma _+=L_{11}`$ and $`\sigma _{}=L_{22}`$. Also $`M_1=\frac{L_{11}L_{22}}{A^2l_{}}`$ and $`M_2=0`$ so that the equation for $`\delta _+`$ decouples from the equation for $`\delta _{}`$. The solution for $`\delta _+`$ is

$$\delta _+(t)=\delta _+(0)e^{_0^t(\sigma _+\frac{}{t^{}}\mathrm{ln}\left[\frac{L_{22}L_{11}}{A^2l_{}}\right])𝑑t^{}}.$$
(29)   
The stability boundary is defined by $`\overline{\sigma }=_0^T(\sigma _+\frac{}{t^{}}\mathrm{ln}\left[\frac{L_{22}L_{11}}{A^2l_{}}\right])𝑑t^{}=_0^T\sigma _+𝑑t^{}=0`$. We have checked that this stability condition agrees with the numerical stability analysis based on Eq. (25) for finite $`Q`$.

We finish by illustrating the re-orientation dynamics of the lamellar structure following a long wavelength instability by direct numerical solution of the governing equation. We focus on the region in which the uniform lamellar structure is linearly unstable. The first example discussed concerns a lamellar structure of initial wavenumber $`(q_1,q_3)=(0.4908,0.4908)`$ being sheared periodically with an amplitude $`\gamma =1`$. Figures 7A,B,C show the sequence of configurations obtained when $`\omega =5\times 10^6`$. A long wavelength transverse modulation of the lamellae is observed (Figure 7A). Subsequent growth leads to the formation of a forward kink band similar to that recently observed experimentally (Figure 7B) . As the strain grows larger, the kink band disappears leaving behind a lamellar structure without any defects and oriented differently relative to the shear (Figure 7C).

In the second example (Figs. 7D,E,F), a structure initially transverse to the flow is being sheared periodically with an amplitude $`\gamma =1`$ and a frequency $`\omega =5\times 10^7`$. A longitudinal perturbation is clearly visible that manifests itself by local compression and dilation of the structure, leading to the disappearance of a pair of lamellae. The overall result is an increase in the lamellar spacing without any change in the orientation.

In summary, we have obtained a nonlinear solution of the model equations that govern the formation of a lamellar structure in the weak-segregation limit. The solution is a periodic lamellar structure with a time-dependent wavevector that adiabatically follows the imposed shear flow, and a time-dependent amplitude which we have computed for the cases of steady and oscillatory shears. In the case of an oscillatory shear, the periodic solution only exists for a range of orientations of the lamellae relative to the shear direction. The width of the region depends on the shear amplitude but not on its frequency. Long wavelength secondary instabilities further reduce the range of existence of stable lamellar solutions. The corresponding stability boundaries depend again on the shear amplitude, but are independent (up to our numerical accuracy) of frequency. We next plan to examine the stability of the nonlinear solution presented in this paper when neither osmotic stresses nor viscosity contrast are neglected.

## Acknowledgments

This research has been supported by the U.S. Department of Energy, contract No. DE-FG05-95ER14566, and also in part by the Supercomputer Computations Research Institute, which is partially funded by the U.S. Department of Energy, contract No. DE-FC05-85ER25000. F.D. is supported by the Microgravity Science and Applications Division of the NASA under contract No. NAG3-1885.

## A Numerical algorithm

We use a pseudo-spectral technique to solve Eq. (3) in two spatial dimensions and in the sheared frame with periodic boundary conditions along both directions. Equation (3) can be written as,

$$\frac{\stackrel{~}{\psi }}{t}=\sigma (t)\stackrel{~}{\psi }q^2(t)\stackrel{~}{\psi ^3},$$
(A1)   
where $`\sigma (t)=q^2(t)q^4(t)B`$ and $`q^2(t)=q_1^2+[a(t)q_1q_3]^2`$. The algorithm we use, due to Cross et al. , is obtained by first multiplying both sides of Eq. (A1) by $`\mathrm{exp}(\sigma (t^{})t^{})`$ and integrating over $`t^{}`$. This gives

$$\mathrm{exp}(\sigma (t)t^{})\stackrel{~}{\psi }_t^{t+\mathrm{\Delta }t}=q^2(t)_t^{t+\mathrm{\Delta }t}𝑑t^{}\stackrel{~}{\psi ^3}(t^{})\mathrm{exp}(\sigma (t)t^{}),$$
(A2)   
where we have assumed $`\sigma (t^{})\sigma (t)`$ and $`q^2(t^{})q^2(t)`$. Next, we write the non-linear term $`\stackrel{~}{\psi ^3}(t^{})`$ as a linear function of $`t^{}`$ in the interval $`tt^{}t+\mathrm{\Delta }t`$, i.e.,

$$\stackrel{~}{\psi ^3}(t^{})\stackrel{~}{\psi ^3}(t)+\frac{\stackrel{~}{\psi ^3}(t+\mathrm{\Delta }t)\stackrel{~}{\psi ^3}(t)}{\mathrm{\Delta }t}(tt^{}).$$
(A3)   
Combining the last two equations finally yields

$`\stackrel{~}{\psi }(t+\mathrm{\Delta }t)`$ $`=`$ $`\mathrm{exp}(\sigma (t)\mathrm{\Delta }t)\stackrel{~}{\psi }(t)q^2(t)\stackrel{~}{\psi ^3}(t)\left[{\displaystyle \frac{\mathrm{exp}(\sigma (t)\mathrm{\Delta }t)1}{\sigma (t)}}\right]`$ (A5)   
$`q^2(t)\left[{\displaystyle \frac{\stackrel{~}{\psi ^3}(t+\mathrm{\Delta }t)\stackrel{~}{\psi ^3}(t)}{\mathrm{\Delta }t}}\right]\left[{\displaystyle \frac{\mathrm{exp}(\sigma (t)\mathrm{\Delta }t)(1+\sigma (t)\mathrm{\Delta }t)}{\sigma ^2(t)}}\right].`$   
Eq. (A5) is first evaluated with the last term on its right-hand side set to zero. The resulting value for $`\stackrel{~}{\psi }(t+\mathrm{\Delta }t)`$ is then used to estimate $`\stackrel{~}{\psi ^3}(t+\mathrm{\Delta }t)`$. Eq. (A5) is finally applied a second time with all three terms on its right-hand side now included in the calculation. The fact that the nonlinear terms are integrated using an explicit procedure in time limits the size of the time step $`\mathrm{\Delta }t`$ that can be used in simulations of the model.

All the numerical results presented in this paper were obtained in the sheared frame of reference with $`128\times 128`$ spectral modes. We have chosen $`B=0.23`$ (which corresponds to $`ϵ=0.04`$) and a time step of maximum size $`\mathrm{\Delta }t=0.2`$, for which no numerical instability was observed. The initial condition $`\psi (𝐫,t=0)`$ is, unless otherwise noted, a lamellar structure obtained by numerical integration of Eq. (A5) with $`a(t)=0`$ (no shear) starting from a random initial condition (a gaussian distribution for $`\psi `$ if zero mean and small variance), for approximately 300,000 iterations until a stationary lamellar structure is reached.