Stripe patterns orientation resulting from nonuniform forcings and other competitive effects in the Swift-Hohenberg dynamics
Daniel L. Coelho, Eduardo Vitral, José Pontes, Norberto Mangiavacchi
SStripe patterns orientation resulting from nonuniform forcingsand other competitive effects in the Swift-Hohenberg dynamics
Daniel L. Coelho, José Pontes, and Norberto Mangiavacchi
GESAR Group/UERJ - State University of Rio de Janeiro,20940-903 Rua Fonseca Teles 121, Rio de Janeiro, RJ, Brazil
Eduardo Vitral
Aerospace Engineering and Mechanics Department, University of Minnesota,107 Akerman Hall, Minneapolis, MN, USA 55455-0153 (Dated: August 4, 2020)Spatio-temporal pattern formation in natural systems presents rich nonlinear dynamics which lead to theemergence of periodic nonequilibrium structures. One of the most successful equations currently availablefor theoretically investigating the behavior of these structures is the Swift-Hohenberg (SH), which contains abifurcation parameter (forcing) that controls the dynamics by changing the energy landscape of the system.Though a large part of the literature on pattern formation addresses uniformly forced systems, nonuniformforcings are also observed in several natural systems, for instance, in developmental biology and in soft matterapplications. In these cases, an orientation effect due to a forcing gradient is a new factor playing a rolein the development of patterns, particularly in the class of stripe patterns, which we investigate through thenonuniformly forced SH dynamics. The present work addresses the stability of stripes, and the competitionbetween the orientation effect of the gradient and other bulk, boundary, and geometric effects taking part inthe selection of the emerging patterns. A weakly nonlinear analysis shows that stripes tend to align with thegradient, and become unstable when perpendicular to the preferred direction. This analysis is complemented bya numerical work that accounts for other effects, confirming that stripes align, or even reorient from preexistingconditions. However, we observe that the orientation effect does not always prevail in face of competing effects.Other than the cubic SH equation (SH3), the quadratic-cubic (SH23) and cubic-quintic (SH35) equations arealso studied. In the SH23 case, not only do forcing gradients lead to stripe orientation, but also interfere in thetransition from hexagonal patterns to stripes.
I. INTRODUCTION
The Swift-Hohenberg equation (SH)[1] is a widely acceptedmathematical model for describing pattern formation in manyphysical systems presenting symmetry breaking instabilities.A classical example where this symmmetry breaking occurs isfound in the emergence of convection stripes in a thin layer offluid heated from below, for which the distance to the onset ofinstability is given by a control parameter or forcing ε . TheSH equation with a cubic nonlinearity (SH3) first appeared inthe framework of Bénard thermal convection between two “in-finite” horizontal surfaces with distinct temperatures, where J.Swift and P. Hohenberg [1] performed a reduction of the fulldynamics, led by the Boussinesq equations, to the slow modesdynamics. Such reduction was also accomplished for reaction-diffusion systems where a similar model could be derived withan additional quadratic nonlinearity (SH23) [2]. In many cases,it is possible to derive this dynamical equation from a potential,a Lyapunov functional, which leads to a relaxational gradientdynamics that has been explored in materials science and softmatter. This potential is often associated with the system’senergy, allowing for non-local diffusive dynamics, pattern se-lection and emergence of dissipative structures. In this context,the SH model is part of the phase-field theory, which originatesfrom statistical mechanics principles, and whose goal is to ob-tain governing equations for an order parameter evolution (e.g.composition, some microstructural feature); it connects thusthermodynamic and kinetic properties with microstructure viaa mathematical formalism [3–5]. An extension of the original Swift-Hohenberg equation(SH3) consists in adopting a destabilizing cubic term and inadding a quintic one. The resulting quintic equation (SH35)admits the coexistence of stable uniform and structured solu-tions, so that localized patterns may exist under a uniform con-trol parameter [6, 7]. This becomes a desirable physical featurein systems that allow for the coexistence of phases of distinctsymmetry [8]. While localized states have been extensivelystudied through SH35, such states can also appear as solutionsfor SH3 and SH23. The issue for SH3 is that at the Turing point, ε =
0, we have a supercritical bifurcation representing the tran-sition from trivial to modulated solutions. This means that theamplitude of the stripes emerging at ε > ε ,and coexistence between the stripes and the trivial solution isnot possible under uniform forcing. However, by allowing aspatial dependency on the control parameter ε , localized statesbecome possible by varying ε between values above and belowthe Turing point. In turn, this poses the question of what are theconsequences of control parameter gradients to pattern selec-tion. While such gradients have been known to induce patternorientation, we here propose a comprehensive numerical studyof these orientation effects, how they affect state localization,and how does gradient orientation fare against other competingorientation effects.Though extensively used for the study of nonequlibrium pat-tern formation, not many works are found before the turn ofthe century, addressing the orientation effect of gradients on astructure of stripes.One of the pioneers in the subject was Walton (1983)[9], a r X i v : . [ n li n . PS ] A ug who considered the onset of convection in the Rayleigh-Bénardproblem, with stress-free upper and lower surfaces, and per-fect insulating sidewalls. Walton assumed a Rayleigh numberweakly above the critical value at one of the sidewalls, andmonotonically decreasing towards the bulk of the convectioncell. Two cases were identified: if the hotter sidewall wasmantained sufficiently above the critical temperature a struc-ture of stripes perpendicular to the sidewall appears. The re-sult was latter canfirmed by the theoretical works of Cross(1982)[10, 11] and several others [12–14]. On the opposite,if the hotter sidewall was mantained slightly above or belowthe critical temperature, a week structure of stripes parallelto the walls emerges. The existence of this weak structurehad already been noted by the experimental work of Sruljes(1979)[15]. However, the following step, concerning the inter-action of the subcritical stripes with a supecritical structure, inthe case where a positive horizontal gradient of temperaure isimposed, was not accomplished. One of the purposes of thepresent work consists in analyzing this problem.At the end of the eighties and the beginning of the ninetiesa group at the Free University of Brussels, pointed to the ten-dency of stripes to align to the gradient of the control parameterand also, to competing effects that, in many cases, dominatethe orientation tendency created by the gradient [16, 17]. Acoexistence between hexagons with opposite phases, and be-tween hexagons with stripes perpendicular to the gradient ofthe control parameter was identified by Hilali et al. (1995)[18],using a Swift-Hohenberg equation with a quadratic term. Mal-omed & Nepomnyashchy (1993)[19] showed that the Lyapunovfunctional associated with a Newell-Whitehead-Segel equationgoverning the evolution of a structure of stripes depends onthe angle between the stripes and the gradient. These authorsshowed that the Lyapunov potential is minimized wheneverthe stripes are parallel to the gradient and attain a maximumwhen perpendicular to it. More recently, Hiscock & Mega-son (2015)[20] considered the orientation of stripe patterns inbiological systems, using the Swift-Hohenberg equation. Theyderived Ginzburg-Landau type equations for the amplitudes ofa structure of two perpendicular modes, one of them being ori-ented along the direction of gradient. Analytically, the steadystate amplitude of both modes was perturbed and it was foundthat the mode with stripes perpendicular to gradient is unsta-ble, and the resulting structure is made of stripes parallel tothe gradient. In addition, these authors added a reaction termto the Swift-Hohenberg equation and found patterns of stripesperpendicular to the gradient of the control parameter. Onlyperiodic boundary conditions (PBC) were considered by theauthors.The present paper extends Hiscock and Megason’s work.We initially address the orientation effect of the gradient of thecontrol paramenter, in absence of other competing effects. Weshow that whenever a preexisting structure of stripes perpendic-ular to the gradient is perturbed in the direction of the gradient,the original structure becomes unstable, and is eventually re-placed by stripes parallel to the gradient (Sec. III). The analysiswas made using weakly nonlinear theory, from which we de-rive a pair of Newell-Whitehead-Segel equations. In sequence,we report the results of our numerical simulations with the SH equation with different nonlinearities, in presence of nonuni-form forcings (Secs. IV and V). Forcings were assumed in theform of spatial ramps, sinusoidal and gaussian distributions ofthe control parameter. We considered Generalized Dirichlet(GDBC) and Periodic Boundary conditions (PBC), and the cu-bic (SH3), quadratic-cubic (SH23) and cubic-quintic versions(SH35) of the Swift-Hohenberg equation. While we investigateorientation effects due to local nonlinearities outside the forcingterm, there is also a recent interest [21] on the nonlocal termsthat could introduce a new lengthscale into the problem (shortor long-ranged). These nonlocal nonlinearities significantlymodify the coefficients of the amplitude equations, and openthe opportunity for future studies on their effects for patternselection and orientation in two-dimensional systems.The analytical results of Sec. III are compared with the nu-merical experiments. Now, due to restriction on finite or peri-odic domain, the gradient faces the competing effect of otherbulk and the boundaries effects, as well as possible discontinu-ities of ε (e.g. ramp in a periodic domain).Competing bulk effects appear in the level of forcing, andin the interaction with modes in all directions, either existingin pseudo-random initial conditions or generated by the non-linearity of the dynamics. If the forcing is sufficiently high,patterns with a high density of defects emerge and dominatethe orientation effect of the gradient. Another bulk effect ap-pears in nonuniform forcings, when somewhere in the domain ε =
0. This situation occurs, for example, in configurations ofthe control parameter resulting in domains where both subcriti-cal and supercritical regions coexist. In this case, the coherencelength of the structure, ξ ∼ ε − / , diverges at ε =
0, making thedomain “short”, in the sense of that critical regions are affectedby the boundaries, no matter how far they are.Boundary effects consist of the well known tendency ofstripes to approach supercritical sidewalls perpendicularly towalls [10–14], and of the less known effect of approachingsubcritical sidewalls in parallel to the walls [9, 17, 22]. We callthis last one by subcritical boundary effect or subcritical effect for short.Though none of these competing effects are taken into ac-count in the weakly nonlinear analysis, they can be investigatedby a complementary numerical study. For instance, we can askif, besides modes parallel to the gradient, there are other un-stable modes in finite or periodic domains, or if modes in any non-perpendicular direction are unstable. All these questionsare addressed in the present paper.The paper is organized as follows: Sec. II shows the Swift-Hohenberg equation and the associated Lyapunov functional.Sec. III presents our weakly nonlinear analysis, where we ad-dress the stability of stripes parallel and perpendicular to thegradient, with respect to perturbations in any direction. Sec. IVpresents the results of the numerical study with the SH3 equa-tion. Sec. V addresses the same effects with the quadratic-cubicand cubic-quintic equations, SH23 and SH35, respectively. Theconclusions are presented next, in Sec. VI. We included twoappendices describing the numerical procedures. Appendix A,describes the finite differences scheme adopted for solving theSwift-Hohenberg equation and lists the parameters used in thesimulations. Further details of the proceedure are given in ourrecent work [23]. Appendix B describes the spectral methodemployed for obtaining steady state amplitude profiles in crosssections of the patterns. II. THE SWIFT-HOHENBERG EQUATION
The SH equation has the so-called gradient dynamics, whichmeans there is a potential, known as a Lyapunov functional,associated with the order parameter field ψ ( x , t ) [with x = ( x , y ) ] that has the property of decreasing monotonically duringthe evolution [24]. It can be derived by using the L -gradientflow of the Lyapunov energy functional: F [ ψ ] = ∫ Ω d x (cid:26) − ε ( x ) ψ + α (cid:2) (∇ + q ) ψ (cid:3) − ζψ + − β ψ + γ ψ (cid:27) ; (1) ∂ t F [ ψ ] = − ∫ Ω d x ( ∂ t ψ ) (cid:54) . (2)where Ω represents the domain whose size is commensuratewith the length scales of the patterns. We consider a regulardomain Ω : { x ∈ [ , L x ] , y ∈ [ , L y ]} . As discussed, ε is acontrol parameter (forcing) that measures the distance from theonset of instability, and it may present a spatial dependency.The constants ζ , β and γ control the energy structure of thesystem, which describes the energy of phases of distinct sym-metry and their stability (with implications on the bifurcationdiagram). The constant α may present different physical in-terpretations (e.g. elastic constant), but it is typically scaledto α =
1. Also, Eq. 2 denotes the behavior of the Lyapunovfunctional, which monotonically decreases until steady state isreached [24]. By taking the variational derivative of Eq. 1 in L norm, the following Swift-Hohenberg equation is obtained: ∂ t ψ = − δ F [ ψ ] δψ ; (3) ∂ t ψ = ε ( x ) ψ − α (cid:16) ∇ + q (cid:17) ψ + ζψ + βψ − γψ . (4)Choosing different values for ζ , β and γ , we obtain the SH3,SH23 and SH35 equations, which we solve in two dimensionalsquare domains, subjected to Generalized Dirichlet Bound-ary Conditions (GDBC), and to Periodic Boundary Conditions(PBC). The numerical method consists of a finite differencessecond order in time and in space Crank-Nicolson scheme witha splitting of operators known as Stabilizing Correction. Thescheme preserves the monotonically decreasing nature of theassociated Lyapunov functional. A description of the schemeis given in [23] and summarized in the Appendix A. III. WEAKLY NONLINEAR ANALYSIS
The stability of preexisting patterns can be investigatedthrough the amplitude equations derived from the Swift-Hohenberg equation. These equations describe the motion of the amplitude that envelopes the oscillating order parameter ψ , and from this coarse-grained description we can evaluatehow perturbations evolve depending on the orientation of theinitial pattern and of the gradient of the control parameter[12, 14, 25]. Since the derivation of the amplitude equations isrelatively similar for SH3, SH23 and SH35, in terms of scalingand considerations between coefficients, we will only providederivation details for the SH3 case ( ζ = γ = ∂ t ψ = ε ψ − α (∇ + q ) ψ + β ψ (5)where β < α > α = x and another in y ψ ( x , t ) = A ( x , t ) e iq x + B ( x , t ) e iq y + c . c . where A and B are complex amplitudes. To performing amultiscale expansion, we introduce slow variables { X , Y , T } which are separated from the fast variables { x , y , t } . The am-plitudes are modulated along the slow variables as A ( X , Y , T ) and B ( X , Y , T ) , while oscillation of the order parameter in thevicinity of q lie in the scale of the fast variables.Assume we have an initial pattern of perfectly aligned stripesin the x direction, that is, stripes presenting a wavenumber q = q j (unit vector in the y direction). By introducing smallperturbations to the wavenumber in Eq. 5, such as q = q + δ q y and comparing terms, we find that for consistency the slowvariables should scale as X = ε / x , Y = ε / y , T = ε t . From this scaling, we note that derivatives in Eq. 5 shouldfollow the chain rule accounting for slow and fast variables, sothat ∂ x → ε / ∂ X , ∂ y → ∂ y + ε / ∂ Y , ∂ t → ε ∂ T . (6)For small ε , the order parameter ψ can be expanded aboutthe trivial solution as ψ = ε / ψ + ε ψ + ε / ψ + . . . By substituting the expanded ψ into Eq. 5 with derivativesacting in multiple scales as defined in Eq. 6, we are able tocollect terms from SH3 in powers of ε . This way, we obtainequations from each order of ε that should be independentlysatisfied, from which we can derive the functions ψ i and theequations governing the evolution of A ( X , Y , T ) and B ( X , Y , T ) .At orders O( ε / ) and O( ε ) , using the notation L c = (∇ + q ) ,we find O( ε / ) : L c ψ = ⇒ ψ = B e iq y + c . c . , A = O( ε ) : L c ψ + L c ψ = ⇒ ψ = B e iq y + c . c . , A = ψ = ε / ψ + ε ψ ψ + ε / ( ψ ψ + ψ ψ ) + . . . Therefore, at order O( ε / ) we find L c ψ = (− ∂ T + + α q ∂ Y − α∂ X − i α q ∂ Y ∂ X + β | B | ) B e iq y + ( . . . ) e i q y + c . c . (7)By rewriting Eq. 7 as L c ψ = θ , the solvability conditionassociated to this equation is that θ must be perpendicularto the null space of L ∗ c : θ ⊥ g ∈ Nu ( L ∗ c ) . This is the Fred-holm’s Alternative, the condition under which the inner product ( θ, g ) = ( ψ, L ∗ c g ) = e iq y , and e − iq y . Therefore, by enforcing solvability weobtain the amplitude equation for B ( X , Y , T ) , ∂ T B = B + α ( q ∂ Y − i ∂ X ) B + β | B | B . We can similarly find the amplitude equation for A ( X , Y , T ) .Since A = A =
0, we need to gather terms up to order O( ε / ) , so that from the Fredholm’s Alternative we obtain, ∂ T A = A + α ( q ∂ Y − i ∂ X ) A + β | B | A . In order to rewrite these amplitude equations in terms of theoriginal variables, we first note that both amplitudes A and B can be expanded in the same form as ψ , that is A = ε / A + ε A + . . . , B = ε / B + ε B + . . . Therefore, accounting for the possibility of a control parameterwith spatial dependence, the pair of amplitude equations forthe two-dimensional SH3 with stripes perpendicular to the y direction is ∂ t A = ε ( x ) A + α ( q ∂ y − i ∂ x ) A + β | B | A , (8) ∂ t B = ε ( x ) B + α ( q ∂ y − i ∂ x ) B + β | B | B . (9)That is, we obtain a system of coupled Newell-Whitehead-Segel(NWS) equations [19].The amplitude equations allow us to evaluate a preexistingpattern stability in the presence of perturbations, and the roleplayed by ∇ ε in such stability. Note that while A and B arecomplex amplitudes, it can be shown that the phase becomesa constant for steady state solutions of parallel stripes [14], sothat we focus on the equation for the real part of the amplitudein the following analysis.Assume ε = ε ( x ) is an increasing ramp in x only, and thatstripes are initially perpendicular to y , with a steady state am-plitude A =
0. The steady real amplitude B satisfies ε ( x ) B − α∂ x B + β B = . Away from the Turing point ( (cid:15) = B ≈ (cid:18) ε ( x )− β (cid:19) / . By introducing small perturbations δ A and δ B to the steadystate solutions, we find from Eqs. 8 and 9 that these perturba-tions evolve as ∂ t ( δ A ) = − ε ( x ) δ A + α ( q ∂ y − ∂ x ) δ A ∂ t ( δ B ) = − ε ( x ) δ B + α ( q ∂ y − ∂ x ) δ B . Since the existence of the solution ( A = , B ) requires ε >
0, this implies that when stripes are parallel to ∇ ε , thesolution ( , B ) is stable with respect to small perturbations δ A and δ B .For the case of stripes perpendicular to the ∇ ε , we keep theramp ε ( x ) in the x direction but change the preexisting patternto stripes aligned in the y direction. The consequence to Eqs. 8and 9 is that A , B and space derivatives swap, and the followingcoupled NWS equations are found ∂ t A = ε ( x ) A + α ( q ∂ x − i ∂ y ) A + β | A | A , (10) ∂ t B = ε ( x ) B + α ( q ∂ x − i ∂ y ) B + β | A | B . (11)Therefore, now we have a steady state solution B = A (cid:44)
0. Due to the ∂ x A derivative in Eq. 10, near the Turinginstability at ε = A will not behave as ( ε ( x )/− β ) / . The A solution satisfying ε ( x ) A + α q ∂ x A + β A = A ( x = ) = A ( x → ∞) = ( ε ( x )/− β ) / . Note that we set the Turing pointat x = ε , the solution is of thetype A ∼ √ ε tanh ( x √ ε / ) , which behaves as A ∼ x ε forsmall x and ε . Since for ε ( x ) = c x , where c is a positiveconstant, there is no analytic solution A , so we assume that A ∼ cx for small x .For stripes perpendicular to the control parameter gradient,perturbations δ B evolve as ∂ t ( δ B ) = ( ε ( x ) + β | A | ) δ B . (12)Taking into account ε ( x ) = cx , β < A ∼ cx for small x ,we conclude the solution ( A , B = ) is unstable. Therefore,while stripes with wavevector q ⊥ ∇ ε are stable with respect tosmall perturbations ( δ A , δ B ) , we find that stripes with q (cid:107) ∇ ε are unstable in the vicinity of the Turing point.In the following sections, we address to the asymptotic heightof the amplitude as h = (cid:18) ε ( x )− β (cid:19) / . (13)This quantity, which appears from the weakly nonlinear analy-sis for the two studied orientation of stripes, will be addressedto as an analytical result, and used for comparison with theattained steady state amplitudes of the numerical results. IV. COMPETITION BETWEEN THE GRADIENT,BOUNDARY AND BULK EFFECTS – SH3
In this section we perform a numerical study of the orienta-tion effect due to gradients of the control parameter in presenceof competing effects, which also complements the stabilityanalysis presented in Sec. III. The study was made througha numerical integration of the SH3 equation, using a finite-difference semi-implicit time splitting scheme that has beenpreviously adopted for Swift-Hohenberg [23, 26] and othernonlinear parabolic equations [27]. The common parametersused were α =
1, and q =
1. Unless otherwise noted, the com-putational domain consisted of 128 × × ∆ x ≈ . [ π /( q )] when using GDBC and ∆ x = π /( q ) when using PBC. We employed a time integra-tion scheme of second order accuracy using a Crank-Nicolsonapproach and the time steps used were ∆ t = . ∆ t = . ε along the x direction. Both Generalized Dirichlet(GDBC) and Periodic (PBC) boundary conditions were con-sidered. Cross sections of selected steady state patterns areshown, along with the envelopes obtained as the steady statesolution of Eqs. 9, 10 and 13.Subsection IV B presents the results of simulations startingfrom pseudo-random initial condition, forced with ramps of ε along the x direction. In Secs. IV C and IV D we describe theresults obtained with sinusoidal and gaussian forcings, respec-tively, both cases starting from pseudo-random initial condi-tions.All patterns shown are at the steady state, with the excep-tions of the initial conditions, and of those shown during theevolution. We consider a pattern at the steady state when itsvelocity of evolution, (cid:219) L , falls below 5 × − (see Sec. A).All simulations starting from pseudo-random distribution of ψ used the same distribution as an initial condition. A. Spatial ramps of ε and preexisting patterns This subsection presents results of eight simulations whoseinitial conditions consist of preexisting structures of straightstripes parallel or perpendicular to the ∇ ε . The forcing takesthe form of spatial ramps along the x direction, submittingthe dynamics to this spatial variation of ε . The results areshown in Figs. 1 and 3. Figure 1 presents the initial conditionsand the steady state of the simulations. Figure 3 shows theone dimensional profile of four patterns from Fig. 1, takenalong the x direction, at the middle height ( y -direction) ofthe domain. For these profiles, we compared the envelopesof modes either parallel or perpendicular to the gradient, toanalytic and numerical estimates based on the weakly nonlinearanalysis detailed in Sec. III.Results shown in the first two rows of Fig. 1 were run froman initial condition consisting of stripes parallel to the gradient,while the last two rows started from stripes perpendicular tothe gradient. GDBC was adopted for the simulations in the first and the third rows, and those in the second and fourth rows wereobtained adopting PBC. The preexisting structure of stripes isshown in the first column of Fig. 1. The configurations (steadystates) shown in the second and third column are numbered forreference.Configuration 1, with a ε ramp increasing from 0 to 0.1,evolved from stripes parallel to the x axis to a bent structure ofstripes approaching the upper and the right sidewalls, orientedperpendiculally to the walls. At the left sidewall, this structureis parallel to wall, a result that complies with the work of Wal-ton (1982)[22], who identified the onset of a weak structure ofstripes parallel to a slightly subcricrical or supercritical side-wall, in presence of a negative gradient of the Rayleigh numberpointing to the bulk of a Rayleigh-Bénard cell. A weak struc-ture of stripes perpendicular to the lower wall is visible closeto that wall, since GDBC favors this orientation due to the zeronormal derivative. Boundary effects dominate both the bulkeffects represented by the initial condition, and ∇ ε .A different situation occurs in the case of configuration 2.The gradient, along with orientation of the initial condition,force the preexisting structure to remain parallel to it, andstripes are kept straightly aligned up to the steady state. Thisresult is attributed to the increase in slope and magnitude ofthe ε ramp, which now increases from 0 up to 0.5.The cases represented by configurations 3 and 4 of Fig. 1were run with PBC. In both cases, the orientation effect ofthe gradient, along with the initial condition and the lack ofthe competition with boundary effects, result in stripes thatremained parallel to the gradient, independently of the forcingmagnitude.Configuration 5 presents a result similar to the one obtainedin configuration 1, while in configurations 6, 7 and 8, thepreexisting initial conditions persist, even with PBC, whereboundary effects are suppressed. The resulting orientationis, in these cases, dominated by the preexisting pattern andboundary effects, opposite to the orientation favored by ∇ ε .The results of configurations 6, 7 and 8 suggest that without acertain level of perturbation, the presence of the gradient is notenough to destabilize stripes in finite or periodic domains, inthe sense that no reorientation by the gradient is observed.While ∇ ε is not sufficient to reorient a preexisting patternwith initial wavevector q parallel to it (without perturbing thesystem), Fig. 2 shows through the evolution of the Lyapunovfunctional from Eq. 1 that the orientation of the pattern withrespect to the gradient strongly affects the relaxational dynam-ics and energy of the steady state structures. The top panel ofFig. 2 follows the energy for the dynamics leading to config-urations 3 and 7 in Fig. 1, and the bottom panel follows theenergy evolution leading to the patterns in configurations 4 and8. For both comparisons it is evident that the configuration ofstripes parallel to the ∇ ε is the one of minimum energy (amongthe two), whose amplitude quickly relaxes to satisfy the controlparameter ramp. However, for q (cid:107) ∇ ε , the relaxation towardsthe steady state is much slower, since this orientation is pe-nalized by the gradient and the steady state pattern presents ahigher associated energy.We also observe that for higher forcing levels and higherramp slope (bottom panel), the system achieves the steady state FIG. 1. The results of eight simulations with the SH3 equation,preexisting structures, and forced with a spatial ramp of the controlparameter. The first column presents the prescribed initial condition.Columns 2 and 3 show the attained steady state. The ramp of thecontrol parameter is given by the diagrams of first row. Rows twoand three correspond to simulations with rigid boundary conditions(GDBC), while results presented in rows four and five were obtainedwith periodic boundary conditions (PBC). faster than for lower forcing levels (top panel). Moreover, theenergy ratio between q ⊥ ∇ ε and q (cid:107) ∇ ε decreases from1 .
080 between the two steady states in the top panel to 1 .
015 inthe bottom panel. The previous two observations suggest thatthe orientational effects due to ∇ ε weaken as the forcing levelincreases.Figure 3 shows cross sections of configurations 1, 2, 3 and6 of Fig. 1. The cross sections were taken along the x direc-tion, at the middle height ( y -direction) of the domain. Theobtained profiles were superposed with the patterns envelope,estimated using two approaches, based on the weakly nonlin-ear analysis presented Sec. III. The first approach consists ofan analytic estimation of the envelope height, using the steady L y a punov f un c ti on a l Configuration 3: q ⊥ ∇ ε Configuration 7: q k ∇ ε L y a punov f un c ti on a l Configuration 4: q ⊥ ∇ ε Configuration 8: q k ∇ ε FIG. 2. Lyapunov functional curves of configurations with rampedforcings with each configuration indicated. It roughly corresponds tothe “normalized” modulus of the time derivative ∂ψ / ∂ t and thereforeis sensitive not only to the growth of the amplitude, but also to thepattern phase dynamics (Appendix A). state solution of the NWS equation given by Eq. 9, obtainedaway from the Turing point assuming small spatial variations.The amplitude solution is given by Eq. 13, which depends onlyparametrically on x in this approximation. Note that no sub-critical solutions are possible with this equation. The secondapproach consists in solving the NWS Eq. 9 at the steady statewith a pseudo-spectral method described in Appendix B.For case ( a ) of Fig. 3 we acquired the profile from configu-ration 1 of Fig. 1, which adopts GDBC. Strictly speaking, theenvelopes for this case should not necessarily fit the crests ofthe pattern, since these envelopes refer to single mode stripes,either parallel or perpendicular to the gradient, and not to bentstripes, where the angle of the wavevector with the gradientcontinuously varies across the domain. Nevertheless, the en-velope estimated with Eq. 13 matches well the pattern crests,except at the boundaries. The amplitude obtained as the steadysolution of Eq. 9 matches well the crests of the pattern. Thissatisfactory agreement suggests that the local amplitude de-pends, in this case, on the local value of ε , and not on the localorientation of the wavevector.Cases ( b ) and ( c ) of Fig. 3 refer to structures of stripes B parallel to the gradient, with GDBC (configuration 02) and PBC(configuration 03), respectively. The cross sections capture theamplitude of the stripes along the direction of the gradient. Theenvelope given by Eq. 13 fits well the amplitude of the patternaway from the boundaries. Moreover, an excellent matchingexists between the envelope obtained as the solution of the NWSEq. 9, and the one extracted directly from the pattern, with bothboundary conditions. Finally, case ( d ) of Fig. 3 was obtainedfrom a preexisting structure of stripes A , perpendicular to the (a) Configuration 01 (GDBC): 0 . (cid:54) ε ( x ) (cid:54) . . (cid:54) ε ( x ) (cid:54) . . (cid:54) ε ( x ) (cid:54) . . (cid:54) ε ( x ) (cid:54) . x -direction, taken at the middle of the height of the domain ( y -direction), and envelopes obtained by the twomethods based on the weakly analysis described in Sec. III. The pattern profiles extracted directly from the integration of the SH3 equationare represented by dotted lines. The envelopes estimated with Eq. 13 are shown in dashed lines, and the envelope obained as the steady statesolutions of the NWS Eqs.8 and 9 are represented in continuous lines. Eq. 13 gives a good estimation of the enveleope, except close to theboundaries, since this equation does not take boundary conditions into account. Note that the so obtained envelopes in case ( a ) shoud notnecessarily fit the crests of the pattern, since these curves refer to single mode stripes, either parallel or perpendicular to the gradient, and not tobended stripes, where the angle of the wavevector with the gradient continuous varies across the domain. Nevertheless, the envelope estimatedwith Eq. 13 matches well the pattern crests, except at the boundaries. gradient (configuration 06). As in case ( a ), the crests of thepattern fit well the envelopes, and reinforces the observationthat away from the boundaries ∇ ε dictates the behavior of theamplitude.In order to evaluate possible reorientation effects due to ∇ ε ,as suggested by Sec. III, we perturb steady state configurationthat originally have a monomodal pattern with q (cid:107) ∇ ε . Theorientation effect can be seen in some of the configurationsfrom Fig. 4, using PBC, so that we avoid effects from therigid imposition on the boundary in the GDBC case. From amonomodal initial condition, we first obtain steady states underdifferent ramps of ε , as seen in configurations 9, 10, 11, and12, which preserve the original preexisting vertical stripes. Byimposing a pseudo-random perturbation of δψ ∈ (− − , − ) to these steady states, the cases with a ramp of ε crossingthe Turing point at ε = ε . This resultagrees with the conclusion from the stability analysis (Eq. 12)as in face of perturbation, the neighborhood of the Turing pointshould become unstable when q (cid:107) ∇ ε , leading the pattern toreorient. However, when the ε ramp stays in the supercriticalregime, ε >
0, the initial pattern did not reorient, as shown inconfiguration 16.
B. Spatial ramps of ε and pseudo-random initial conditions In this subsection we present the results of four simulationsrun from the same pseudo-random initial condition, and thecross section for two of the obtained configurations. Two sim- ulations were performed with GDBC, and the other two withPBC. The results are sumarized in Fig. 5.In the case of configuration 17, the orientation effect ofthe gradient is dominated by the boundary and the subcriticaleffects: a bent structure of stripes perpendicular to the rightand to the upper sidewalls emerges. This structure persists inthe subcritical region, with weak stripes approaching the leftsidewall in parallel orientation. The result is in agreement withthe works of Sruljes (1979)[15] an Walton (1982)[22].The steady state pattern developed in configuration 17 is sim-ilar to the one appearing in configuration 1 of Fig. 1, with theorientation effect of the gradient dominated by the boundaryand the subcritical effects. The structure developed in configu-ration 18 is similar to the one of configuration 17, with stripesperpendiculars to the supercritical lower and right sidewalls,and paplelels to the critical left sidewall. Additionally an weakstructure of stripes perpendiculars to the upper supercriticalsidewall is visible.In the case of PBC, configurations 19 and 20, no boundaryeffects are present and the resulting structure selects a directionalmost parallel to the the gradient, even with the presence ofmodes in every directions in the initial condition. A Benjamin-Feir instblility appears at the left limit of the periodic structureof configuration 20.Figure 6 shows in the first row the cross section of configu-ration 17 from Fig. 5. In the second row, we move the Turingpoint from 1/2 to 1/4 of the domain length, in order to eval-uate if a translation of the forcing may result in any pullingof reorientation of the pattern. We used the steady pattern ofconfiguration 17 as the initial condition, for the configurationshown in the second row. In both cases, the cross section was (a)(b)FIG. 4. The results of eight simulations with the SH3 equation, preexisting structures, forced with a spatial ramp of the control parameter andPBC. The ramp of the control parameter is given by the diagrams of first row. In the first row are shown (a) steady patterns obtained from theindicated initial condition of preexisting vertical stripes without pseudo-random perturbation. In the second row are shown (b) Steady patternsobtained from the steady states above with the exact same pseudo-random perturbation (ranging from − − to 10 − ). Periodic domains (PBC)were chosen in order to observe “only” interactions between bulk and ∇ ε effects without the rigid imposition on the boundary (GDBC). Sincein the weakly nonlinear analysis we find that vertical stripes are unstable in the presence of the ∇ ε , a small perturbation was sufficient to thesystem evolve into a new steady pattern where stripes align to the ∇ ε direction (Configurations 13 and 15). Configuration 14 does not presentsuch behavior and the preexisting stripes orientation prevails. taken at the middle height ( y -direction) of the domain. Wecompare the ψ profile with the estimated amplitude of mode B , evaluated by the steady state solution of the NWS Eq.9, andalso with the analytic approximation from Eq. 13.As in configuration 1 profile from Fig. 3, the evaluatedenvelopes refer to modes parallel to the gradient and not tobent stripes of the pattern, where the angle between the localwavevector and the gradient continuously varies across the do-main. In spite of this fact, the envelopes qualitatively fit thepeaks of the pattern, suggesting that the height of the envelopefar of boundaries depends primarily on the local value of ε andnot on the direction of the wavevector.Note that the solution of the NWS Eq. 9 presents an esti-mate for the envelope on the subcritical region. The resultsshow that a translation of the forcing ramp expands the patterntowards the new location of the Turing point, keeping the orig-inal bent form of the pattern in configuration 17. Therefore,bulk and boundary effects still win over gradient effects, andno reorientation is observed. C. Sinusoidal forcings
This subsection presents the results of eigth simulations runfrom pseudo-random initial conditions with sinusoidal forc-ings in x . Both GDBC and PBC were prescribed. Sinusoidalforcings are of interest because they allow for multiple subcrit-ical and supercritical regions in a single domain, with multipleTuring points in x . Also, when compared to ramps, periodicforcings better accommodate PBC, so that undesired effectsdue to a jump of ε at the boundary are not an issue. The resultsare summarized in Fig. 7, where we present the resulting steadystate patterns (labeled as configurations 13-22). The first rowdisplays the distribution of ε forcings, the second patterns ob-tained by prescribing GDBC, and the third those obtained byprescribing PBC.Results shown in configurations 21 to 25 (GDBC) and 26to 30 (PBC) of Fig. 7 evolved either to patterns parallel togradient, or at least, with regions where the stripes are parallelto the gradient. The orientation effect clearly appears in thesesimulations, even when the restrictive GDBC are prescribed.To obtain configuration 21, we used a low amplitude sinu-soidal distribution of ε with subcritical regions ( − . ≤ ε ≤ . (a) GDBC(b) PBCFIG. 5. The results of four simulations with the SH3 equation, forcedwith a spatial ramp of the control parameter ε . All simulations startedfrom the same pseudo-random initial conditions, shown in the first col-umn. The remaining columns present the steady state. First row: theprescribed profile of the control parameter ε . Second and third rows:GDBC, and PBC boundary conditions, respectively. Simulations ofthe second row show that boundary effects dominate the orientaioneffect when GDBC are prescribed. In the third row, the absence ofboundary effects allows the dominance of the orientation effect of thegradient. cal regions close to the right and left walls where no patternemerges. A higher amplitude of the forcing, as depicted in con-figuration 30, leads to the emergence of stripes aligned to thegradient in all supercritical regions, so that a higher ε allowsto overcome energy penalizations due to boundary conditions.Configuration 23 of Fig. 7 was run with a sinusoidal forcingadded to a constant, so that (cid:15) (cid:62)
0. No subcritical regions arepresent in this simulation. A weak structure of small stripesperpendicular to the upper and lower walls emerges at the cen-ter of these walls. Several Benjamin-Feir instabilities are alsoobserved. Configuration 24 was run with a similar sinusoidalforcing, but of higher amplitude. Due to the higher forcing,the structure can more easily accommodate defects. A patternof winding stripes appear at the central region, with two focusdefects showing at the top and at the bottom of this region. Wenote that the winding form of these stripes comes from the factthat they anchor perpendicularly to the upper and lower walls,while approaching regions of ε close to zero with a parallelalignment. This is in agreement with our previous observationfrom Fig. 2, that orientational effects due to ∇ ε become lessprevailing as the magnitude of the forcing increases.Configuration 25 is obtained by increasing the minimum (cid:15) even further, so that we have a fluctuating forcing of highmagnitude and GDBC. The orientational effect of the gradient (a) configuration 17 with Turing point located at L / L / ε where the Turing point is located at 1/4 ofthe domain length. The same maximum value of ε prescribed forconfiguration 1 was adopted for the simulation shown in the secondrow of the present figure. Though associated to straight stripes alignedto the gradient, the envelopes fit well the crests of the patterns, whichconsist of bended stripes, with the angle between the local wavevectorand the gradient continuously varying across the domain. is fully dominated by the bulk and boundary effects. A patternof diagonal stripes with a high density of defects emerges.The results presented in the third row of Fig. 7 were run withPBC. In the case of configuration 26, a low forcing and the lackof boundaries prevent the emergence of defects. The orienta-tional effect of the gradient prevails and a structure of stripesparallel to the gradient emerges in all supercritical regions.Configuration 27 is similar to the previous, but the sinusoidalforcing presents a higher amplitude. This higher forcing insupercritical regions allows for stripes that deviate from thegradient alignment, and we observe columns of stripes that al-ternate between parallel and inclined alignments. The steadystate for this case strongly depends on the initial distribution of ψ , and once a column of inclined stripes is formed, it is unableto completely reorient in the gradient direction.Configuration 28 of Fig. 7 presents again a case where theforcing consists of a low amplitude with zero minimum ε . Theabsence of sidewalls and the relatively low forcing weakenscompeting effects and pattern aligns accordingly to the gra-dient. The resulting structure is parallel to the gradient andBenjamin-Feir instabilities are observed at the neighborhoodof the Turing point. For configuration 29 we use the same forc-ing as in configuration 24. The steady state pattern is similar theone of configuration 27, but with stripes occupying the entiredomain, as the forcing is non-negative. Lastly, configuration30 starts from the same forcing as configuration 25, and theresulting winding pattern shows that, even for PBC, ∇ ε fails toorient the stripes whenever the magnitude of the force remainsat large magnitude.Fig. 8 shows a cross section of the steady state pattern from0 (a) GDBC(b) PBCFIG. 7. The results of eight simulations with the SH3 equation, forced with a spatial sinusoidal profile of the control parameter ε . All simulationsstarted from the same pseudo-random initial conditions, shown in the first column. The remaining columns present the steady state, attainedwhen (cid:219) L ≤ × − . First row: the prescribed profile of the control parameter ε . Second and third rows: GDBC (a), and PBC (b) boundaryconditions, respectively.FIG. 8. Cross section along the x -direction, of the pattern shown inconfiguration 22 of Fig. 7 (dotted line). The cross section was taken atthe middle of the height ( y -direction). To the profile obtained directlyfrom the pattern we superposed the envelope of mode B , obtained bytwo methods: as the steady state solution of Eqs. 10 and 11, adopting asinusoidal distribution of ε along the x direction and as an estimationof the envelope profile, using Eq. 13. configuration 30 of Fig. 7. The cross section profile was takenat the middle of the height ( y -direction) and is represented bya a line with small circles. To this profile we superposed thepattern envelope of mode B estimated from two approaches:the first one consisted in the steady state solution of Eqs. 10and 11, for the amplitude B , using a sinusoidal distribution of ε . The so evaluated envelope is represented by a continuousline in Fig. 8. The second evaluation of the envelope was doneby using Eq. 13 also adopting a sinusoidal distribution of ε .The envelope is represented by a dashed line in Fig. 8. Thetwo estimated envelopes fit the central part of the pattern crosssection, and deviate at the boundaries, which are not taken into FIG. 9. Pattern evolution from pseudo-random initial conditions,subjected to PBC and a diagonal sinusoidal distribution of the controlparameter, given by: ε ( x ) = . [ q ( x + y )] . The structure evolveinto stripes aligned in the ∇ ε direction until the steady state is reached. account in the derivation of those two curves.Fig. 9 presents the result of a configuration consisting ofa pseudo-random distribution of ψ as the initial condition, asinusoidal forcing along the domain diagonal, and PBC. Lackof boundary effects along with the orientation effect of thegradient, and existence of modes along all direction in the initialcondition lead to a pattern of stripes parallel to the diagonal.Figure 10 shows the (cid:219) L × t curves of selected configurationsshown in Fig. 7. These curves present an irregular region at thevery begining of the simulations, when the patterns emergesfrom the pseudo-random initial condition. Most of the patterngrowth occurs at this phase. As a result, (cid:219) L decreases by someorders of magnitude. The evolution proceeds with changes inthe phase, and with the amplitude essentially saturated. (cid:219) L evolves irregularly at much lower level, with peaks occuring1 ˙ L [-] FIG. 10. Comparison between (cid:219) L curves of configurations 22(GDBC) and 27 (PBC) with sinusoidal forcings. It roughly corre-sponds to the “normalized” modulus of the time derivative ∂ψ / ∂ t andtherefore is sensitive not only to the growth of the amplitude, but alsoto the pattern phase dynamics (Appendix A). at the collapse of defects. This phase is followed, in all cases,by a linear (exponential) decrease of (cid:219) L . We assume that thepattern reached a steady state when (cid:219) L attains the value 5 × − . We mention that also the Lyapunov potential decreasesexponentially at this phase.Fig. 11 presents a case of a preexisting structure of stripesalong the y -direction and a sinusoidal profile of ε along thediagonal of the domain, using PBC. Despite lacking the restric-tive effect of boundary conditions, the gradient is dominatedby the initial condition, and the preexisting structure persists.Upon adding a noise δψ ∈ (− − , − ) to the initial condi-tion, the preexisting structures is destabilized and replaced bya sinusoidal distribution of stripes parallels to the gradient. D. Gaussian forcings, and pseudo-random initial conditions
By imposing gaussian forcings we observe another bulk ef-fect that competes with the gradient in orienting the stripes.Fig. 12 shows the steady state patterns obtained in four sim-ulations, two of them run with a sharper circular gaussiandistributions of the control parameter, centered at the middleof the domain and two, with a wider gaussian forcing. GDBCand PBC were considered for each configuration of ε . Thefour simulations started from the same pseudo-random initialcondition adopted in all cases presented in Secs. IV B and IV C.The adopted gaussian distribution is given by: ε ( x ) = Ae − R (( x − x ) + ( y − y ) ) . (14)In the first case, configurations 31 (GDBC) and 25 (PBC),we used a sharper distribution of ε with parameters R given inTab. II ( R and R ). For both, the resulting pattern takes theform of a target, with stripes presenting a wavevector parallelto the gradient, q (cid:107) ∇ ε , a completely opposite situation withrespect to several other cases run from pseudo-random initialconditions using forced with ramps or sinusoidal distributionsof ε . The orientation effect of the gradient does not appearin this case, and orientation is dominated by a geometric bulk (a)(b)(c) + ˙ L [-] (d)FIG. 11. Pattern evolution from preexisting vertical stripes, subjectedto PBC and a diagonal sinusoidal distribution of the control parameter,given by: ε ( x ) = F cos [ q ( x + y )] , where q = . q . For (a) F = . F = .
5, the preexisting structure persists, anddominates the orientation effect of the gradient. (c) Upon adding aperturbation in the form of a uniform distribution ranging from − − to 10 − to the initial condition, the preexiting structure collapsesand is replaced by stripes parallels to the gradient. (d) The timeevolution of (cid:219) L for the latter is shown and this simulation proceededuntil (cid:219) L (cid:54) × − . effect. Due to the disk-form of ε in the two-dimensional do-main, a target pattern is the one that fills the supercritical regionwhile minimizing defects, which is a geometric compatibilityeffects. Otherwise, if stripes were to orient accordingly to the ∇ ε , the resulting pattern would cointain a large ammount of de-fects (dislocations) to accomodate such orientation, increasingsignificantly the energy of the configuration. Therefore, target2 (a)(b)FIG. 12. The results of four simulations with the SH3 equation,forced with a gaussian distribution of ε . All simulations started fromthe same pseudo-random initial conditions, shown in the first column.Configurations 31 and 32 were run with GDBC, while PBC wereprescribed for configurations 33 and 33. patterns minimizes the Lyapunov potential, in spite of beingpenalized by the control parameter gradient.A second case with a wider gaussian forcing is shown inFig. 12, configurations 32 (GDBC) and 33 (PBC), with param-eters R given in Tab. II ( R and R ). This case corresponds to aforcing with a sharper distribution of Fourier modes, thereforea smaller range of modes persists in the steady state pattern. Weobserve that a pattern of stripes with wavevector aligned withthe diagonal of the domain appears. The pattern extends fora longer distance in this diagonal direction, even invading thesubcritical region, for both GDBC nad PBC. The effect is dueto the fact that the hardest direction for modulation of the am-plitude occurs in the direction of the wavevector, whereas theeasiest direction is the perpendicular one. This property of pe-riodic patterns results in amplitudes modulated in compliancewith Newell-Whitehead-Segel equations [28, 29]. V. COMPETITION BETWEEN THE GRADIENT,BOUNDARY AND BULK EFFECTS – SH23 AND SH35
In this section, we briefly explore two other forms of the SHequation with additional nonlinearities, in order to numericallyevaluate the effects of altering a bifurcation diagram on theorientation of resulting patterns. As mentioned before, Eq. 4can be addressed to as SH23 for γ =
0, and ζ, β (cid:44)
0. Anal-ogously, we refer to it as SH35 for ζ =
0, and γ, β (cid:44) ∇ ε . For small positive ε , hexagonal patterns are the minimum energy state, whichdestabilizes when ε is increased and stripe patterns becomeenergetically favored. The coexistence of both structures witha nonuniform forcing was addressed by Hilali et al. [18], wherestripes formed with q (cid:107) ∇ ε , contrary to our results for SH3. Inorder to clarify if such effect observed for SH23 was inducedby initial conditions (an initial ramp in ψ , in their case), andassess how ∇ (cid:15) interferes in the hexagon to stripe transition, weperform simulations for SH23 using different initial conditionsand ramps for (cid:15) . Numerical results are shown in Fig. 13, using ζ = .
65 and β =
1, in which we observe the possibility of co-existence between hexagons and stripes. The nonuniform forc-ings considered were the following ramps: − . (cid:54) ε ( x ) (cid:54) . . (cid:54) ε ( x ) (cid:54) .
5. Three configurations for each forcingwere considered, starting from pseudo-random initial condi-tions and preexisting patterns (horizontal and vertical stripes,respectively).Configurations 43 and 44 started from pseudo-random initialconditions, and configurations 37, 38, 39 and 40 had theirpreexisting condition perturbed with a pseudo-random noiseranging from − − to 10 − with an uniform distribution. Inthe absence of perturbations, the initial condition is preserved, i.e , hexagon patterns do not emerge. The results show that, inthe presence of a ∇ ε (cid:44)
0, the preexisting patterns are unstable,and we observe regions where stripes decay to ψ = ε = − . ε = . q ⊥ ∇ (cid:15) appear in the region of positive ε , with aweakly formed hexagonal structure close to ε =
0. Configura-tion 37 show that when a preexisting structure of stripes with q ⊥ ∇ (cid:15) is perturbed, stripes in the positive ε region remainperfectly aligned to the gradient, and no transition to hexagonsis observed. In configuration 39 we see that by perturbing aninitial condition of stripes with q (cid:107) ∇ (cid:15) , the remaining stripesdid not reorient according to ∇ ε , and (opposite to configuration37) a well formed column of hexagons appeared for regions ofsmall ε . This is a consequence of the higher energy associatedto stripes when q (cid:107) ∇ (cid:15) , as compared in Fig. 2, so that for small ε a stripe to hexagon transition is promoted. For configurations44, 38 and 40, where the ramp ranges from ε = ε = . ∇ ε this gradient has a weaker ef-fect on inducing stripe pattern alignment. In configuration 44,we do not see alignment starting from a pseudo-random initialcondition, while in configuration 38, even though the preexist-ing pattern was made of stripes with q ⊥ ∇ (cid:15) , hexagons stillemerged in the 0 < ε < .
25 region (opposite to configuration37).Finally, we present numerical results for the SH35 equationin the presence of a control parameter gradient, using β = FIG. 13. The results of six simulations with the SH23 equation,forced with a spatial ramp of the control parameter ε indicated in thefirst row. The parameters of the SH23 were: ζ = . β = γ =
1. Configurations 37, 38, 39 and 40 were perturbed so that theremaining stripes share the domain with hexagons, as expected forthis range of parameters.
3, and γ =
1. The ψ = ε c = − β / γ for which bothstripes and ψ = ε > ε c , stripes are energetically favored, whilefor ε < ε c the equilibrium state is ψ =
0. A consequenceof the subcritical bifurcation is that even in ε > ε c regions,stripes do not form from a pseudo-random initial condition forfinite values of ε c . Therefore, in the SH35 case we adopt asquare pattern as initial condition, of the type ψ = A cos ( q x ) + B cos ( q y ) , in order to evaluate the ∇ ε effect on filtering thepattern.For obtaining configurations 41 and 42, we use GDBC, whilefor configurations 43 and 44, we use PBC. For the chosen setof parameters ε c ≈ − .
52. With a ramp ranging from ε = − ε =
0, for both GDBC (away from the boundary) andPBC the y direction mode was filtered, and the resulting struc-tures present stripes with q ⊥ ∇ ε for ε > ε c , and ψ = ε < ε c . By changing the ramp to ε = − . ε =
0, we seein configuration 44 (PBC) that away from the boundaries thestructure perfectly aligns according to ∇ ε . However, in con- (a) GDBC(b) PBCFIG. 14. The results of four simulations with the SH35 equation,forced with a spatial ramp of the control parameter ε . All simulationsstarted from the same preexisting square pattern ( ψ = A cos ( q x ) + B cos ( q y ) ), shown in the first column. The remaining columns presentthe steady state. First row: the prescribed profile of the controlparameter ε . Second and third rows: results for GDBC, and PBC,respectively. Results of the second row show that boundary effectsdominate the orientaion effect when GDBC are prescribed. In thethird row, the absence of boundary effects allows the dominance ofthe orientation effect of the gradient. figuration 42 we note that by decreasing the ramp inclination,orientation effects due to boundary conditions become strongerin comparison to ∇ ε , as they mostly dictate the resulting pat-tern. From our simulations using SH35 with PBC, we observedthat patterns tend to orient more strongly according to ∇ ε thanin the case of SH3 or SH23, which is presumably associated tothe subcritical nature of the bifurcation. VI. CONCLUSIONS
This work addressed the stripe orientation effect due to gra-dients of the control parameter ε in the Swift-Hohenberg dy-namics, and how it fares against competing effects. A weaklynonlinear analysis confirms and extends existing results show-ing that stripes tend to align to the gradient. We show thatstripes with wavevector q perpendicular to ∇ ε are stable andcorrespond to a lower energy state than stripes with q parallelto ∇ ε , which are unstable near the Turing point. However, ourmain numerical results show that the orientation effect of thegradient, though existing, does not always prevail when facingcompetition with other bulk, boundary, geometric, and periodiceffects due to computational domains. This competition leadsto the emergence of a rich dynamics, as apparent in our results,4which strongly depends on the magnitude of the forcing andinitial conditions (preexisting patterns or pseudo-random). Inthis sense, the various forms of forcing and initial conditionsaddressed in this work, while extensive, do not cover the fullrange of possible cases.Among our most representative results, we cite simulationswhere, starting from a structure of stripes with wavevector q perpendicular to ∇ ε , the initial pattern is preserved, and caseswhere this orientation is compromised by competing effects.Even starting from a structure of stripes with q parallel tothe gradient, we found cases where the initial pattern persists.However, upon introducing a finite pseudo-random perturba-tion on the initial conditions, either ∇ ε fully reorients the stripesor competing effects interfere, reshaping the pattern (e.g. in-troducing defects). Forcing with a gaussian distribution of ε resulted, in one case, in a target structure of stripes perpendic-ular to the gradient ( q (cid:107) ∇ ε ) and, in other, in straight stripesperpendicular to the domain main diagonal. In the former, thetarget pattern results from a geometric constraint, as a patternwith q ⊥ ∇ ε would imply in a large density of defects (higherenergy). In the latter, the elongated pattern along the domainmain diagonal complies with a well known effect, which statesthat the most difficult direction for modulation is the one of thewavevector [28, 29].An estimation for the local amplitude of the patterns couldwas derived in the form of a NWS equation, function ofthe nonuniform forcing (cid:15) ( x ) . Its steady state solution, h = [ ε ( x )/− β ] / (Eq. 13), neglects diffusion of the amplitude andstill represents a satisfactory approximation of its envelope farfrom the boundaries and the Turing point. This estimation ismore accurate when the pattern preserves a dominant mode inthe x (or y ) direction, but even when there are other modespresent, the approximation is satisfactory.In this work, we observed that the orientation effect of thegradient prevailed in many, but not all, of the explored config-urations, either starting from a preexisting structure of stripesparallel or perpendicular to the gradient, or with a pseudo-random initial condition. The orientation effect is relevant formany physical systems presenting periodic patterns, such as indevelopmental biology [20, 39], smectic mesophases [8], andlocalized sand patterns [40], whose dynamics have been studiedby Swift-Hohenberg type equations, but present mechanismsof stripe orientation that are not well understood. ACKNOWLEDGMENTS
The authors thank FAPERJ (Research Support Foundationof the State of Rio de Janeiro) and CNPq (National Councilfor Scientific and Technological Development) for the finan-cial support. Daniel Coelho acknowledges a fellowship fromthe Coordination for the Improvement of Higher EducationPersonnel-CAPES (Brazil). A FAPERJ Senior Researcher Fel-lowship is acknowledged by J. Pontes. The authors dedicatea special thanks to prof. D. Walgraef, from the Free Univer-sity of Brussels, who early pointed to the orientation effects ofnonuniform forcings, on patterns of stripes. This endeavor is based on his ideas.
APPENDIX A: FINITE DIFFERENCE SCHEME FORSOLVING THE BIDIMENSIONAL SWIFT-HOHENBERGEQUATION
The cubic, quadratic-cubic and cubic-quintic (SH3, SH23,SH35, respectively) equations were solved by the finite differ-ences method in uniform and structured meshes, with a semi-implicit Crank-Nicolson scheme. The scheme is uncondition-ally stable, features truly second order representation of timeand space derivatives, and strict representation of the associatedLyapunov functional, which monotonically decays. Details ofthe scheme are given by Christov and Pontes (2002) [24] andby Coelho et al. (2020)[23]. The scheme is summarized below,for the sake of completeness.
Target scheme
In order to construct a Crank-Nicolson second order in timenumerical scheme, we adopt the following discrete representa-tion of Eq. 4: ψ n + − ψ n ∆ t = (cid:2) ε ( x ) − α q − α q ∂ x − α q ∂ y − α∂ x − α∂ x ∂ y − α∂ y + ζ (cid:0) ψ n + (cid:1) + ( ψ n ) + β (cid:0) ψ n + (cid:1) + ( ψ n ) − γ (cid:0) ψ n + (cid:1) + ( ψ n ) (cid:35) (cid:18) ψ n + + ψ n (cid:19) . (A1)The superscript ( n + ) refers to variables evaluated at the nexttime step, and n , to the ones evaluated at the current one. TheRHS of Eq. A1 consists of a nonlinear operator (in braces)actuating on the order parameter ψ , evaluated as the meanvalue at the middle of the time steps ( n ) and ( n + ) . The termsin Eq. A1 are grouped in the target scheme, as follows: ψ n + − ψ n ∆ t = (cid:16) Λ n + / x + Λ n + / y (cid:17) ψ n + + ψ n + f n + / . (A2)For the SH23 and SH3 equations ( γ = Λ n + / x , Λ n + / y and f n + / are defined as: Λ n + / x = (cid:34) − α (cid:32) ∂ x + q (cid:33) − β (cid:0) ψ n + (cid:1) + ( ψ n ) (cid:35) ; Λ n + / y = (cid:34) − α (cid:32) ∂ y + q (cid:33) − β (cid:0) ψ n + (cid:1) + ( ψ n ) (cid:35) ; f n + / = (cid:104) ε ( x ) − α (cid:16) q ∂ x + q ∂ y + ∂ x ∂ y (cid:17) + ζ (cid:0) ψ n + (cid:1) + ( ψ n ) (cid:35) (cid:16) ψ n + + ψ n (cid:17) , ζ = Λ n + / x = (cid:34) − α (cid:32) ∂ x + q (cid:33) − γ (cid:0) ψ n + (cid:1) + ( ψ n ) (cid:35) ; Λ n + / y = (cid:34) − α (cid:32) ∂ y + q (cid:33) − γ (cid:0) ψ n + (cid:1) + ( ψ n ) (cid:35) ; f n + / = (cid:104) ε ( x ) − α (cid:16) q ∂ x + q ∂ y + ∂ x ∂ y (cid:17) + β (cid:0) ψ n + (cid:1) + ( ψ n ) (cid:35) (cid:16) ψ n + + ψ n (cid:17) . Terms are assigned to Λ n + / x and Λ n + / y in order to constructnegative definite operators that assure the uncondional stabilityof the scheme. The scheme is nonlinear and requires internaliterations at each time step. Terms with superscript ( n + ) inthe function f n + / and in both operators, are evaluated at thecurrent internal iteration, whereas the remaining ones ( ψ n + in Eq. A2) are the truly implicit ones, evaluated at the nextinternal iteration. The splitting
The original scheme is replaced by two following equations,using the Stabilizing Correction scheme. The two equations areequivalent to the target scheme up to a second order error in theevaluation of the time derivative, and does not change the steadystate solution. In addition, the two-equations scheme minimizememory requirements and truncations errors [23, 24]: (cid:101) ψ − ψ n ∆ t = Λ n + / x (cid:101) ψ + Λ n + / y ψ n + f n + / + ( Λ n + / x + Λ n + / y ) ψ n ; ψ n + − (cid:101) ψ ∆ t = Λ n + / y ( ψ n + − ψ n ) , where (cid:101) ψ is an intermediary estimation of ψ at the new time step.The second equation provides a correction and thus we obtain ψ for the new time. Variables in the RHS of both equations areevaluated as the mean value between the time step n and theestiamtion for the new one, at the previous internal iteration.All spatial derivatives are represented by second order centraldifferences formulæ. Steady state criteria
We follow the structure evolution by assessing the rate ofchange in time of the pattern during the simulation by moni-toring (cid:219) L , the relative norm rate of change defined as: (cid:219) L = ∆ t (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) n x (cid:205) i n y (cid:205) j | ψ n + i , j − ψ ni , j | n x (cid:205) i n y (cid:205) j | ψ n + i , j | (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) . (A3) TABLE I. Final times required to reach the steady state criteria foreach configuration.
Config. Steady State Config. Steady State t = t = t = t = t =
297 25 t = t =
82 26 t = t = t = t = t = t = t = t = t = t = t = t = t = t = t = t = t = t = t = t = t = . t = t = t = t = t = t = t = t = . t = . t = t = t = t = t = . t = t = . It roughly corresponds to the ratio between the spatial averageof the modulus of time derivative ∂ψ / ∂ t and the spatial averageof the modulus of the function itself. Futhermore, it is sensitivenot only to the growth of the amplitude, but also to the patternphase dynamics. The simulations proceeded until (cid:219) L (cid:54) × − , which is our criterion for reaching the steady state. The discrete implementation of the Lyapunov fuctional
The Lyapunov functional associated to the SH equation isimplemented through the discrete formula derived by Chris-tov & Pontes (2001) [24] for the SH3, and extended for theSH35 by Coelho et al. (2020) [23]. The formula presents a O (cid:0) ∆ t + ∆ x + ∆ y (cid:1) approximation of the functional given byEq. 2:6 F n + − F n ∆ t = − n x (cid:213) i = n y (cid:213) j = (cid:32) ψ n + i , j − ψ ni , j ∆ t (cid:33) ; F n = n x (cid:213) i = n y (cid:213) j = (cid:34) − (cid:15) (cid:16) ψ ni , j (cid:17) − ζ (cid:16) ψ ni , j (cid:17) − β (cid:16) ψ ni , j (cid:17) + γ (cid:16) ψ ni , j (cid:17) + α q (cid:16) ψ ni , j (cid:17) (cid:35) − α q n x (cid:213) i = n y (cid:213) j = (cid:34) ψ ni + , j − ψ ni , j ∆ x (cid:35) + (cid:34) ψ ni , j − ψ ni − , j ∆ x (cid:35) + (cid:34) ψ ni , j + − ψ ni , j ∆ y (cid:35) + (cid:34) ψ ni , j − ψ ni , j − ∆ y (cid:35) + α n x (cid:213) i = n y (cid:213) j = (cid:34) ψ ni + , j − ψ ni , j + ψ ni − , j ∆ x + ψ ni , j + − ψ ni , j + ψ ni , j − ∆ y (cid:35) . (A4) TABLE II. Parameters adopted in the numerical study throughout thiswork.
Parameter Formulæ Value Description q - 1.0 Critical wavenumber λ π / q π Critical wavelength w x , w y - 8 Wavelengthsper domain length g r - 16 Grid resolution n x , n y w x × g r , 128 Nodes per mesh side w y × g r N n x × n y ×
128 Total mesh nodes L x , L y w x λ , w y λ ≈ . L ) ∆ x , ∆ y L /( n − ) ≈ . ∆ x , ∆ y L / n ≈ . ∆ t - 0.5 Time step (SH3) ∆ t - 0.1 Time step (SH23/SH35) A - 0.2 Gaussian maximumvalue (peak) R n − x - Configs. 04 and 09 R . n − x - Configs. 05 and 10 x , y L x / L y / The monotonic decay of the finite differences version is en-forced, provided that the internal iterations converge [24].
Parameters adopted in the simulations
The simulations shared common parameters regarding thewavenumber, domain sizes and time steps. Those are pre-sented in the following table along with the parameters usedin particular simulations, such as the gaussians distributionsparameters.Space and time steps were chosen following [23] for a goodcompromise between accuracy and computational cost (perfor-mance). The grid resolution represents the number of nodesper critical wavelength.
APPENDIX B: PSEUDO-SPECTRAL SCHEME FORSOLVING THE ONE DIMENSIONALNEWELL-WHITEHEAD-SEGEL EQUATION
In the weakly nonlinear analysis section, a pair of coupledNewell-Whitehead-Segel (NWS) equations (Eqs. 8 and 9) isderived for the modes A ( x , t ) and B ( x , t ) via the multiple scaleformalism. In order to compare the amplitude envelopes fromSH simulations and the ones described by those amplitudeequations, we develop numerical solutions for the NWS. Theone-dimensional ( x -direction) NWS equation has the form, ∂ t u = ε ( x ) u − α∂ x u + β u , (B1)where α = β = − t (cid:62) u ≡ u ( x , t ) is a real func-tion described in the regular domain Ω : { x ∈ [ , L x ]} withperiodic boundary conditions (PBC) and generalized Dirichletboundary conditions (GDBC). Since an analytical study is nottrivial for this equation, we develop a numerical study using asemi-implicit pseudo-spectral method with first-order accuracyin time. The Fourier approach was adopted for the configura-tions with periodic boundary conditions (PBC). The Chebyshevapproach was adopted for the configurations subjected to gen-eralized Dirichlet boundary conditions (GDBC), for dealingwith the imposed boundary conditions. Both approaches arebriefly discussed in the following sections. Fourier pseudo-spectral scheme
The eigenfunctions of the fourth-order differential operatorover the domain with periodic boundary conditions (PBC) arethe Fourier modes e ik · x ( k ∈ Z N ). Since ∂ x e ik · x = | k | e ik · x ,equation B1 can be written as ∂ t ˆ u k = ε ( x ) ˆ u k − α | k | ˆ u k + β ˆ u k , (B2)where ˆ u k is the Fourier coefficient associated with the mode k .The fourth order derivative term is treated implicitly since it hasa numerical stabilizing property (denoted by index n + n ). The nonlinear terms arecomputed in real space in order to avoid computing Fouriermode convolutions (higher computational effort) and thereforeare treated explicitly. Since we are interested in the steady statesolution, a semi-implicit first order accurate in time schemewas employed, that can be expressed as follows:ˆ u n + k = µ k (cid:16) ˆ u nk + ∆ t ˆ f n (cid:17) , (B3)where µ k = (cid:0) + α ∆ tk (cid:1) − and ˆ f n = F ( f ( u n )) is theFourier transform of the nonlinear and variable coefficientterms f ( u n ) = ε ( x ) u n + β ( u n ) . This transformation is per-formed via a fast Fourier transform (FFT)-based code withoutdealiasing, using Octave FFT library. Chebyshev pseudo-spectral scheme
Equation B1 is solved using Chebyshev spectral collocationmethod [30]. Chebyshev polynomials of degree n have n zerosin the interval ξ ∈ [− , ] that should be mapped to the physicaldomain x ∈ [ , L x ] . For that purpose, a simple mapping ischosen: x = . ( ξ + ) L x . The numerical scheme can beexpressed as (cid:16) I + α ∆ tD x (cid:17) u n + = u n + ∆ t f n , (B4)where f n = ε ( x ) u n + β ( u n ) and D x is the Chebyshev collo-cation for the fourth-order differential operator on the mappeddomain. This system of linear equations is solved subjected tothe boundary conditions for the amplitude B . For such purposewe adopted B = ∂ B / ∂ x =
0, for the boundary points x = x = L x . These conditions are consistent with the GDBC usedin the SH equation simulations, for the order parameter ψ . [1] J. Swift and P. C. Hohenberg, Physical Review A (1977),10.1103/PhysRevA.15.319.[2] D. Walgraef, Spatio-Temporal Pattern Formation: With Exam-ples from Physics, Chemistry, and Materials Science , 1st ed.,Partially Ordered Systems (Springer-Verlag New York, 1997).[3] N. Provatas, M. Greenwood, B. Athreya, N. Goldenfeld, andJ. Dantzig, International Journal of Modern Physics B , 4525(2005).[4] K. R. Elder, M. Katakowski, M. Haataja, and M. Grant, Physicalreview letters , 245701 (2002).[5] K. R. Elder and M. Grant, Physical Review E , 051605 (2004).[6] H. Sakaguchi and H. R. Brand, Physica D , 274 (1996).[7] J. Burke and E. Knobloch, Physical Review E , 056211 (2006).[8] E. Vitral, P. H. Leo, and J. Viñals, Physical Review E ,032805 (2019).[9] I. Walton, Journal of Fluid Mechanics , 455 (1983).[10] M. Cross, The Physics of Fluids , 936 (1982).[11] M. Cross, Physical Review A , 1065 (1982).[12] M. C. Cross and P. C. Hohenberg, Rev. Modern Phys. , 851(1993).[13] H. S. Greenside and W. M. Coughran, Physical Review A ,398 (1984).[14] P. Manneville, Non-equilibrium Structures and Weak Turbulence (Academic Press: San Diego, CA, USA, 1990).[15] J. A. Sruljes,
Zellularkonvection in Behältern mit HorizontalenTemperauregradienten , Ph.D. thesis, Fakultät für Machinenbau,Univ. Karlsruhe, Karlsruhe (1970).[16] J. Pontes,
Pattern formation in spatially ramped Rayleigh-Bénardsystems , Ph.D. thesis, Free University of Brussels, Brussels, Bel-gium (1994).[17] J. Pontes, D. Walgraef, and C. I. Christov, Journal of Computa-tional Interdisciplinary Sciences , 11 (2008).[18] M. F. Hilali, S. Métens, P. Borckmans, and G. Dewel, PhysicalReview E , 2046 (1995).[19] B. A. Malomed and A. A. Nepomnyashchy, Europhysics Letters(EPL) , 195 (1993).[20] T. W. Hiscock and S. G. Megason, Cell systems , 408 (2015).[21] D. Morgan and J. H. Dawes, Physica D , 60 (2014). [22] I. Walton, Studies in Applied Mathematics , 199 (1982).[23] D. L. Coelho, E. Vitral, J. Pontes, and N. Mangiavacchi, Sub-mitted to Journal of Computational and Applied Mathematics(2020).[24] C. Christov and J. Pontes, Mathematical and Computer Mod-elling (2002), 10.1016/S0895-7177(01)00151-0.[25] R. B. Hoyle, Pattern formation: an introduction to methods (Cambridge University Press, 2006).[26] C. Christov, J. Pontes, D. Walgraef, and M. Velarde, ComputerMethods in Applied Mechanics and Engineering (1997),10.1016/S0045-7825(96)01176-0.[27] E. Vitral, D. Walgraef, J. Pontes, G. Anjos, and N. Man-giavacchi, Computational Materials Science (2018),10.1016/j.commatsci.2018.01.034.[28] A. C. Newell and J. A. Withehead, J. Fluid Mech. , 279 (1969).[29] L. A. Segel, Journal of Fluid Mechanics , 203 (1969).[30] J. P. Boyd, Chebyshev and Fourier spectral methods (CourierCorporation, 2001).[31] G. Dewell and P. Borckmans, Physics Letters A , 189 (1989).[32] V. E. Podneks and I. W. Hamley, Journal of Experimental andTheoretical Physics Letters , 617 (1996).[33] D. Boyer and J. Viñals, Physical Review E (2002),10.1103/PhysRevE.65.046119.[34] N. N. Yanenko and M. Holt, The Method of Fractional Steps (Springer, 1971).[35] J. Douglas and H. H. Rachford., Trans. Amer. Math. Soc. (1956),10.1090/S0002-9947-1956-0084194-4.[36] N. Ghoniem and D. Walgraef,
Instabilities and Self-organizationin Materials (Oxford Univ. Press, 2008).[37] U. M. B. Marconi and P. Tarazona, Journal of Physics: Con-densed Matter , A413 (2000).[38] J. S. Langer, Reviews of modern physics , 1 (1980).[39] M. Ruppert, F. Ziebert, and W. Zimmermann, New Journal ofPhysics , 052001 (2020).[40] A. Auzerais, A. Jarno, A. Ezersky, and F. Marin, Physical ReviewE94