A new resonance mechanism in the Swift--Hohenberg equation with time-periodic forcing
AA new resonance mechanism in the Swift–Hohenberg equationwith time-periodic forcing
Punit Gandhi ∗ and Edgar Knobloch Department of Physics, University of California, Berkeley CA 94720, USA
C´edric Beaume
Department of Aeronautics, Imperial College London,South Kensington, London SW7 2AZ, UK (Dated: October 15, 2018)
Abstract
The generalized Swift–Hohenberg equation with a quadratic-cubic nonlinearity is used to studythe persistence and decay of localized patterns in the presence of time-periodic parametric forcing.A novel resonance phenomenon between the forcing period and the time required to nucleate onewavelength of the pattern outside the pinning region is identified. The resonance generates distinctregions in parameter space characterized by the net number of wavelengths gained or lost in oneforcing cycle. These regions are well described by an asymptotic theory based on the wavelengthnucleation/annihilation time near the boundaries of the pinning region. The resulting theory leadsto predictions that are qualitatively correct and, in some cases, provide quantitative agreementwith numerical simulations. ∗ punit [email protected] a r X i v : . [ n li n . PS ] N ov . INTRODUCTION Spatially localized structures in physical, chemical and biological systems [1, 2] oftenconsist of a time-independent spatial pattern embedded in a homogeneous background. Thetheory behind the origin and properties of such structures is well understood, at least in onespatial dimension. In this theory the localized structures are described in terms of heterocliniccycles connecting the homogeneous state to the patterned state and back again [3, 4]. Suchcycles may be structurally stable, and so persist over an interval of parameter values locatedwithin the region of bistability between the homogeneous and spatially periodic states. Thisso-called snaking region [5, 6] typically contains two or four families of homoclinic solutionsconnecting the homogeneous state to itself, organized within a snakes-and-ladders structure[7, 8]. These correspond to spatially localized states of ever greater length and accumulateon the heteroclinic cycle as their length increases. Examples of this behavior have beenidentified in both gradient and nongradient systems, including buckling of slender structures[6, 9], shear flows [10], doubly diffusive convection [11–13], porous media convection [14] androtating convection [15], among others.Localized states are also encountered in systems with a fluctuating or noisy background[16–19] as well in periodically driven systems [20–25]. Temporal focing has, in general, anumber of consequences. In extended systems it may destabilize existing patterns or leadto resonant excitation of new patterns and a variety of phase-locking phenomena [26, 27].In addition, new structures may be generated by rapid switching between two coexistingattractors [28, 29]. Localized structures may be impacted in two different ways. First,the temporal forcing may render existing localized structures time-dependent, and second,it may generate bistability between a homogeneous state and an extended parametricallydriven spatially periodic pattern. The latter case creates a parameter regime where spatiallylocalized time-dependent patterns may be found [30].In the present paper we focus on time-independent systems supporting spatially localizedstates and study the effect of time-periodic forcing on these states. Our interest in this typeof problem is motivated in part by recent studies of the growth of vegetation patterns nearthe transition to desertification [31–35]. Simple models of this process predict the presenceof patchy patterns [35] and the properties of such patterns are yet to be examined whenseasonal variation in growth conditions is included.2ost of the systems mentioned above can be modeled using equations of Swift–Hohenberg-type despite the fact that this equation is of gradient type. This is because the snakes-and-ladders structure of the snaking region in gradient and nongradient systems is identical(although the stability properties of the solutions may differ [36]). Consequently, we adopthere a model of this type and investigate the effect of time-dependent forcing on the existenceand stability of localized states within this model. We identify a number of new structuresin this system, including time-dependent breathing states and structures that grow or shrinkin an episodic manner. In particular, we identify a novel resonance phenomenon betweenthe forcing period and the time required to nucleate a new wavelength, triggered wheneverthe forcing parameter falls outside the pinning region. This resonance leads to a complexpartitioning of the parameter space whose structure can be understood qualitatively, and insome cases quantitatively, using appropriate asymptotics.This paper is organized as follows. In the next section we summarize pertinent resultsconcerning the autonomous Swift–Hohenberg equation with competing quadratic and cubicnonlinearities. In section III we consider the effect of high frequency temporal forcing ofthis equation, and then in section IV we focus on the different breathing states present forintermediate frequencies. This section forms the bulk of the paper. Low frequency forcing isconsidered in section V, followed by brief conclusions in section VI.
II. THE SWIFT–HOHENBERG EQUATION
The quadratic-cubic Swift–Hohenberg equation (SHE) serves as a model for patternformation in a broad range of physical systems. This equation which, in one dimension, takesthe form u t = ru − (cid:0) ∂ x (cid:1) u + bu − u , (1)describes the dynamics of a real field u ( x, t ) in time. The parameter r specifies the strengthof the forcing while the parameter b > (cid:112) /
38 determines the extent of the bistability regionbetween the homogeneous state u h ≡ u p ( x ), u p ( x ) = u p ( x + 2 π )for all x . The equation can be written in terms of a Lyapunov functional F [ u ], referred to as3he free energy , u t = − δ F [ u ] δu , F [ u ] = 1Γ Γ / (cid:90) − Γ / − ru + (cid:2) (1 + ∂ x ) u (cid:3) − b u u dx. (2)Thus, on a domain of finite spatial period Γ all initial conditions approach a steady statecorresponding to a local minimum of the free energy.In order to study the effects of time-dependence, we write r = r + ρ sin ωt , where r , ρ , and T = 2 π/ω define the offset, amplitude, and period of the oscillation. Note thatthis type of parametric forcing leaves the homogeneous state u h = 0 unchanged. In thefollowing we take b = 1 . π (i.e. 40 characteristic wavelengths), unless otherwise noted. In addition we imposethe symmetry x → − x of eq. (1) on all solutions of the system thereby focusing on even solutions. This procedure allows us to perform computations on the half domain. Weintegrate the equation forward in time using a fourth order exponential time differencingscheme [39] on an equidistributed mesh. Our calculations are performed in Fourier spaceand fully dealiased. In cases where a larger domain was necessary, the spatial density of gridpoints was kept constant. Steady state solutions of the constant forcing case were computedusing the numerical continuation software AUTO [40]. A. Stationary localized states
For b = 1 . r ≡ r , a spatially periodic solution u p bifurcates subcritically from u h at r = 0. The periodic state passes through a saddle-node at r sn ≈ − . u h in − . < r <
0. The Maxwellpoint is located at r M ≈ − . F [ u p ] = F [ u h ] = 0. The pinning or snaking region r − < r < r + straddles this point( r − ≈ − . r + ≈ − . L ) or minima (hereafter L π ) at x = 0 as described in [7]. In a finite domain snaking continues until the domain is (almost)filled with pattern; thereafter the solution branches exit the pinning region and terminate onbranches of periodic states near their saddle-node. Throughout this study we present ourresults in terms of the amplitude of the pattern, A = max x ( u ), and the location x = f of4 || u || A + D + P ± D − A − r + r − r sn (a) f A sn A − A + A / (b) x u (c) FIG. 1. (a) Bifurcation diagram showing the normalized L -norm || u || = (cid:113) (cid:82) Γ / − Γ / u dx of time-independent solutions of eq. (1) as a function of the forcing parameter r . Vertical dashed linesdelimit the amplitude regime A − , the depinning regimes D ± and the pinning region P ± . Thecharacteristics of each regime are described in the text. (b) The same as (a) but projected onthe amplitude A = max x ( u ) and the position f > u ( x ) corresponding to the red circles in (a) and (b), with black dashed lines indicatingthe locations x = ± f of the fronts. the front connecting the pattern to the homogeneous state relative to the axis of symmetry x = 0 of the pattern, f = 2 (cid:82) Γ / xu dx (cid:82) Γ / u dx . (3)As the amplitude A of snaking localized solutions is comparable to that of the periodicstate at the same parameter values, larger values of f indicate broader localized structures.However, between the pinning region and r = 0 the solutions broaden to fill the availabledomain as their amplitude decreases to zero. Thus f increases without bound as A → . Temporal dynamics of localized initial conditions Spatially localized initial conditions of eq. (1) eventually settle on a steady state, but thetype of steady state and the transient leading to it depend on r and the initial condition.The relevant regimes organized around the presence of steady spatially localized states canbe identified in fig. 1 and are summarized below: • Regime A − : r < r sn . Only the trivial state is stable. The dynamics is dominated by anoverall amplitude (or body) mode and the amplitude of any localized initial conditiondecays homogeneously to zero. • Regime D − : r sn ≤ r ≤ r − . Two stable states are present: u h and u p , with F [ u h ] < F [ u p ]. Spatially localized initial conditions evolve via a depinning (or edge) moderesponsible for the progressive loss of spatial periods while keeping their amplitudeconstant. The solution collapses to the trivial state only when its extent becomescomparable to one wavelength. • Regime P ± : r − ≤ r ≤ r + . There is a large number of coexisting stable and unstablestates: trivial, spatially periodic and spatially localized with different numbers ofperiods. The long-time behavior of the system is determined by the basins of attractionof the stable states and hence by the initial conditions provided. • Regime D + : r + ≤ r ≤
0. The situation is similar to that in D − but this time F [ u p ] < F [ u h ]. Spatially localized initial conditions nucleate additional wavelengthsunder the influence of the depinning mode, and in periodic domains evolve into thespatially periodic state. • Regime A + : r >
0. The only stable state is the spatially periodic state.These regimes are depicted in the phase portraits in fig. 2. The computations use L localizedsolutions from the snaking region, hereafter u ( x ), as initial conditions. These evolve first in A to the appropriate amplitude, followed by depinning if r falls outside the pinning region.If r = r ± + δ , where | δ | (cid:28)
1, the resulting front propagates at an overall constantspeed determined by the nucleation time T dpn ∝ | δ | − / computed in reference [7]. In thiscalculation, the solution takes the form u ( x, t ) = u ( x ) + (cid:112) | δ | a ( t ) v ± ( x ) + O ( | δ | ), where6 a) A − , r = − . D − , r = − . P + , r = − . D + , r = − . FIG. 2. Space-time plots (left panels) and sample phase space trajectories (right panels) illustratingthe dynamics of localized solutions of L type in the different parameter regimes in fig. 1, initializedwith different values of r . Green dots indicate stable periodic states for the given forcing while bluedots indicate stable localized states. The purple region shows the pinning region. a is the time-dependent amplitude associated with the eigenmode v ± that is responsible7or triggering a nucleation (+) or annihilation ( − ) event. The equation that governs thedynamics of a is α ˙ a = (cid:112) | δ | ( α sgn( δ ) + α a ) , (4)where sgn( δ ) is the sign of δ and the coefficients α j for each of the two cases ( r ± ) are computednumerically from the following integrals: α = Γ / (cid:90) v ± ( x ) dx, α = Γ / (cid:90) v ± ( x ) u ( x ) dx, α = Γ / (cid:90) v ± ( x ) ( b − u ( x )) dx. (5)We now discuss the solutions describing a nucleation event near r + , where α > r − , where α <
0. Within the pinning region, δ <
0, a pair of stable and unstable steady state solutions u is present, corresponding to thevicinity of a fold on the right of the snaking branch L (fig. 3(b)) and all initial conditionsapproach the stable state or diverge. Outside of the pinning region, δ >
0, there are nostable solutions and the amplitude a → ∞ for all initial conditions. The upper panels offig. 3(a,c) show typical trajectories a ( t ) corresponding to the dynamics represented by thearrows in fig. 3(b), with the right panel showing three successive nucleation events. Weapproximate the time between depinning events, i.e. the nucleation time T dpn+ , as the timeinterval between successive asymptotes where a diverges.Near the fold of the snaking branch L and to its right (0 < δ (cid:28)
1) the system undergoesdynamics on the time scale δ − / (eq. (4)) and thus T dpn+ = O ( δ − / ). Upon leaving thevicinity of the fold (i.e. when a → ∞ ), the system transitions towards the next fold on thesnaking branch (fig. 3) before slowing down again. This transition corresponds to a nucleationevent that adds a wavelength to each side of u ( x ) and the process repeats at successivefolds. Since the structure of v + ( x ) is almost independent of the length 2 f of the localizedstate (it is an edge mode) the resulting process is periodic, a fact that can be highlightedby introducing the Riccati variable b defined by a = − α ˙ b/α (cid:112) | δ | b . In terms of b , eq. (4)becomes the oscillator equation ¨ b + δ Ω b = 0, where Ω ≡ α α /α >
0. The lower panels offig. 3(a,c) show the oscillator amplitude b ( t ) corresponding to the a ( t ) solutions shown justabove. However, care must be taken in the interpretation of this equation, because b → a → ∞ while eq. (4) breaks down already when a = O (1). Despite this caveat weshall find the use of the variable b useful since it highlights the possibility of a temporalresonance when the system (1) is forced with a time-periodic forcing. A similar discussion8 a t b (a) δ < t a t b (c) δ > FIG. 3. Amplitude of the depinning mode in terms of a ( t ) (upper panels) and b ( t ) (lower panels)near r + for (a) δ < δ >
0; (b) shows the corresponding bifurcation diagram. Figure (c)shows three successive nucleation events corresponding to times where a ( t ) → ∞ (top panel) or b ( t ) = 0 (lower panel). applies to annihilation events near r = r − .Figure 4 compares the leading order theoretical prediction T dpn+ = π/ √ δ Ω + obtained fromeq. (4) (dashed lines) with numerical simulations of eq. (1). The theory works well for δ (cid:28) T dpn ) − = (cid:80) σ n δ n/ , n ≥
1, andcomputing the σ n using a least squares fit. Given the results in fig. 4 we determined thata fifth order calculation was sufficiently accurate even when r ∼ r sn . The figure shows thenucleation time in D + (red) and annihilation time in D − (blue). The symbols representresults from simulations, the dashed lines represent the prediction from the leading orderasymptotic theory, while the solid lines represent the fifth order numerical fit. The times T col for a marginally stable periodic state at r sn (black crosses) and a localized state at r − (blackdiamonds) to reach the trivial state by amplitude decay in A − are shown in black. The blackdashed line represents the leading order asymptotic theory applied to the periodic state near r sn . The coefficients σ n for both the asymptotic theory and the fifth order numerical fit forour choice of parameters are summarized in table I. We will find that the numerical fit isrequired for quantitative agreement between the theory presented in section V and numerical9imulations presented in section IV. However, the theory cannot be applied to localized statesin A − since these states undergo both amplitude decay and depinning. Ω σ σ σ σ σ T dpn+ T dpn − T colper T colloc - 0.2081 0.4431 2.962 -34.15 79.52 TABLE I. Values of the coefficients σ j determined from a least squares fit of the depinning/collapsetime to simulations with constant forcing of the form T − = (cid:80) n =1 σ n | r − r ± ,sn | n/ . The frequencyΩ is calculated numerically from the integrals in eq. (5) in each case. III. THE HIGH FREQUENCY LIMIT
We begin our study of the effects of a time-periodic forcing on localized states by consideringthe limit of fast oscillations. We first consider the case when the frequency of the forcingcycle is so fast that the motion of the fronts does not permit nucleation/annihilation ofadditional/existing periods. We then increase the amplitude of the forcing cycle so that thestructure remains unpinned for an appreciable amount of time.
A. The averaged system
The qualitative behavior of eq. (1) is unchanged when the forcing frequency is high enoughthat insufficient time is spent outside of the pinning region for depinning to occur. The effectof the periodic forcing in this case is small, producing rapid amplitude fluctuations of theexisting localized states. We introduce the effective Maxwell point ¯ r M using the relation (cid:104)F [ u ] (cid:105) = 0 , (6)where the brackets indicate an average over the forcing cycle. We assume that the periodicforcing occurs at a high frequency, ω → ω/(cid:15) , where (cid:15) (cid:28) φ = ωt/(cid:15) .We seek solutions of eq. (1) in the form u ( x, t ) = u ( x, t, φ ) + (cid:15)u ( x, t, φ ) + (cid:15) u ( x, t, φ ) + . . . ,10 r T r sn r − r + T dpn+ T dpn − T col (a) √ δ T − (b)(c) FIG. 4. (a) The time between nucleation events ( T dpn+ , red crosses), annihilation events ( T dpn − ,blue circles) and the time to collapse to the trivial state ( T colper , black crosses) as functions of theparameter r , starting from marginally stable L solutions at r = r + , r = r − and the periodicstate at r = r sn , respectively. The symbols show results from direct numerical simulations, thesolid lines are fits to this data, and the dashed lines are predictions from the leading order theoryin reference [7]. The corresponding results for the collapse time for a localized state at r − arealso shown ( T colloc , black diamonds). (b) Comparison between numerical data (circles/crosses) andthe leading order theory (dashed lines), showing ( T dpn ± ) − as a function of the square root of thedistance δ from the pinning region. The corresponding fifth order fits are shown using solid lines.(c) A space-time representation of a simulation at r ≈ − . δ ≈ . r + , with red representing high values and blue low values of u ( x ). The solid black line shows the instantaneous front position f ( t ). ωu φ = (cid:15) (cid:104) ( r + ρ sin φ ) u − (cid:0) ∂ x (cid:1) u + bu − u − u t (cid:105) , (7)where t is the original timescale on which the averaged dynamics take place, and assumethat ρ, r , b, ω = O (1). The leading order equation ( u ) φ = 0 gives u ( x, t, φ ) = A ( x, t ). Atorder O ( (cid:15) ), we obtain ω∂ φ u = ( r + ρ sin φ ) u − (cid:0) ∂ x (cid:1) u + bu − u − ∂ t u . (8)The solvability condition requires that the integral over a single period of the fast oscillationof the right side of eq. (8) vanishes. This condition yields the governing equation for A : ∂ t A = r A − (cid:0) ∂ x (cid:1) A + bA − A . (9)Thus, in the limit of a high frequency forcing cycle with order unity amplitude, the leadingorder behavior follows the time-independent SHE. Corrections arise at second order, as wenow show.Equations (8) and (9) show that the O ( (cid:15) ) correction to the leading order behavior is givenby u ( x, t, φ ) = − ρω cos φA ( x, t ) + A ( x, t ) . (10)At O ( (cid:15) ) we obtain ω∂ φ u = ( r + ρ sin φ ) u − (cid:0) ∂ x (cid:1) u + 2 bu u − u u − ∂ t u (11)leading to the solvability condition ∂ t A = r A − (cid:0) ∂ x (cid:1) A + 2 bA A − A A . (12)Similarly, the second order correction to the solution takes the form u = ρ ω cos 2 φA ( x, t ) − ρω sin φ (cid:2) bA ( x, t ) − A ( x, t ) (cid:3) (13) − ρω cos φA ( x, t ) + A ( x, t ) , and we obtain, at O ( (cid:15) ), ω∂ φ u = ( r + ρ sin φ ) u − (cid:0) ∂ x (cid:1) u + 2 bu u − u u + bu − u u − ∂ t u , (14)yielding ∂ t A = r A − (cid:0) ∂ x (cid:1) A + b (2 A A + A ) − A A + A A ) − (cid:0) ρω (cid:1) A . (15)12e can define an averaged variable with error at order O ( (cid:15) ) that describes the dynamics onthe long timescale: A ≡ π π (cid:90) (cid:0) u + (cid:15)u + (cid:15) u (cid:1) dφ = A + (cid:15)A + (cid:15) A . (16)On summing the solvability conditions, we obtain the following equation for the dynamics ofthe averaged variable ∂ t A = r A − (cid:0) ∂ x (cid:1) A + bA − (cid:2) π ( ρT ) (cid:3) A + O ( T ) , (17)where, for clarity, we have introduced the period of the forcing cycle T ≡ π(cid:15)/ω . The resultis a SHE with a modified cubic term.We find that the averaged Maxwell point of the system, defined by eq. (6), is in fact theMaxwell point of the averaged system (17). This can be checked explicitly by noting that¯ F [ A ] = (cid:104)F [ u + (cid:15)u + (cid:15) u ] (cid:105) + O ( (cid:15) ) , (18)where ¯ F is the free energy of the averaged system with periodic forcing, F is the free energyof the system with a constant forcing r , and the average is over a forcing cycle. Furthermore,¯ F [ A ] = F [ u ] + ρ T π Γ Γ / (cid:90) − Γ / u dx + O ( T ) , (19)implying that the free energy in the fluctuating system is greater than that of the systemwith constant forcing r . We can use this expression to calculate the frequency-inducedshift of the Maxwell point explicitly by finding the value of r where ¯ F [ A ] = 0. Because theperiodic forcing has increased the energy of the spatially periodic state, the Maxwell pointof the averaged system necessarily shifts to the right (¯ r M > r M ) to compensate while theboundaries of the pinning region also shift to the right. Following [7] we obtain¯ r ± = r ± + ρ T π (cid:82) Γ / − Γ / u v ± dx (cid:82) Γ / − Γ / u v ± dx . (20)Here u is the marginally stable solution of the constant forcing system at r ± , and v ± are theeigenmodes of the linearized problem at r ± responsible for wavelength nucleation/annihilation.Both integrals are positive and we find that2¯ p ≡ ¯ r + − ¯ r − ≈ p − . ρT ) , (21)13here p = ( r + − r − ) / shrinks the width of the pinning region andshifts it to the right. B. Large amplitude forcing
We may repeat the above calculation in the case ρ → ρ/(cid:15) so that we are now dealing withfast oscillations with a large amplitude. This allows the system to spend enough time outsideof the pinning region for depinning to take place. However, in this limit, a large fraction ofthe forcing cycle is actually spent in the amplitude decay regimes A ± , and thus the leadingorder dynamics will not be comprised simply of nucleation and annihilation events. Thisregime is described by the equation ωu φ − ρ sin( φ ) u = (cid:15) (cid:104) r u − (cid:0) ∂ x (cid:1) u + bu − u − u t (cid:105) , (22)where r , b, ρ = O (1), and we look for solutions in the form u = u + (cid:15)u + (cid:15) u + . . . . Atleading order we obtain ( ωu ) φ − ρ sin( φ ) u = 0, with solution u ( x, t, φ ) = e − ( ρ/ω ) cos φ A ( x, t ) . (23)At O ( (cid:15) ), the governing equation becomes ω∂ φ u − ρ sin( φ ) u = r u − (cid:0) ∂ x (cid:1) u + bu − u − ∂ t u . (24)Imposing the requirement that u is periodic on the fast timescale leads to the solvabilitycondition (cid:82) RHS e ( ρ/ω ) cos φ d φ = 0, where RHS stands for the right-hand-side of eq. (24), andan evolution equation for A : ∂ t A = r A − (cid:0) ∂ x (cid:1) A + bI ( ρω ) A − I ( ρω ) A . (25)Here I ( x ) is the modified Bessel function of the first kind. Thus the slowly evolving amplitude A of the leading order solution u satisfies a SHE with modified coefficients. A higher ordercalculation shows that there is no additional correction at O ( (cid:15) ) and the averaged dynamicsfollows SHE with the modified nonlinear coefficients of eq. (25) up to O ( (cid:15) ).In contrast to the ρ ∼ O (1) case, the cubic nonlinearity is now dramatically increasedrelative to the quadratic one. This results in a rapid decrease in the region of bistability as14 /ω increases. Indeed, at ρ/ω ≈ .
02 the region of bistability disappears in a codimensiontwo point where the bifurcation that creates the periodic state transitions from subcritical tosupercritical.
IV. INTERMEDIATE FREQUENCIES: BREATHING LOCALIZED STRUCTURES
We now move away from the high frequency limit and investigate parameter combinationsthat permit depinning. For this purpose we consider parameter excursions that allow thesystem to traverse P ± and spend a significant time in both D + and D − , i.e., we take r − < r < r + , and ρ > p , where p ≡ ( r + − r − ) / A. The fate of stable localized initial conditions
Figure 5 shows sample results for ρ = 0 . r = − .
28, in each case starting from the samestable spatially localized L solution of the time-independent problem r ≡ r . The figureshows that, depending on T , the solution can undergo growth/decay through a depinning-likeprocess (fig. 5(a,c)), decay to the trivial state via an amplitude mode (fig. 5(d)), or takethe form of a periodic orbit corresponding to a localized solution with no net motion of thefronts (fig. 5(b)). Moreover, the growth/decay of new wavelengths can occur regularly fromone period of the forcing to the next, or in a seemingly irregular way. In particular, fig. 5(a)shows a growth scenario for T = 50 in which the solution grows in length by one wavelengthon each side after approximately three cycles of the forcing. This process is irregular in thesense that the number of nucleation and decay events is not constant from one period ofthe forcing to the next. It is also interesting to note that this simulation does not reach thespatially periodic state, but instead approaches an oscillating state with a defect at the edgeof the periodic domain. In contrast, fig. 5(c), obtained for T = 250, shows a very regularpattern of five nucleations events followed by six decay events during the course of eachforcing cycle, resulting in an overall decay of the state. Finally, fig. 5(d) for T = 350 showsan initial growth phase followed by abrupt amplitude decay to the trivial state.Although the wavelength of a localized solution depends on the forcing parameter r ,15 A sn A − A + / (a) T = 50 f A sn A − A + / (b) T = 150 f A sn A − A + / (c) T = 250 f A sn A − A + / (d) T = 350 FIG. 5. Space-time plots (left panels) and the corresponding phase space trajectories (right panels)for solutions of eq. (1) with r ( t ) = − .
28 + 0 . πt/T ), b = 1 .
8, initialized using an L solutionat r = − .
28. The red dashed lines in the right panels correspond to evolution past the time windowrepresented in the left panels, while the green lines represent spatially periodic solutions of thetime-independent case. The period T is indicated below each plot. The trajectory in (a) terminateson a time-periodic defect state.
16t is always near the preferred wavelength 2 π , and thus f undergoes abrupt jumps byapproximately 2 π whenever the fronts depin. These jumps are most evident during thegrowth phase since the time between nucleations is longer than the time between annihilations.The results shown in fig. 5 are independent of the length of the initial stable state selected forthe simulation, provided that 6 π (cid:46) f (cid:46) Γ / − π throughout the simulation, i.e., providedthe structure remains well localized. We likewise report that the initial phase of oscillationdoes not affect the dynamics.We distinguish periodic orbits from growing and decaying orbits using the instantaneousfront velocity V f defined by the relation V f ≡ ˙ f . We look at the averaged front velocity (cid:104) V f (cid:105) over a cycle period (calculated for oscillation periods for which 6 π (cid:46) f (cid:46) Γ / − π anddisregarding the first oscillation period), and consider an orbit periodic if |(cid:104) V f (cid:105)| < × − (corresponding to no net nucleation or annihilation event within a time period of about3 × units). For decaying and growing orbits, the average change in the front position (cid:104) ∆ f (cid:105) over one cycle helps distinguish the regular behavior in fig. 5(c) from the irregulardynamics in fig. 5(a): regular dynamics translate into (cid:104) ∆ f (cid:105) close to an integer number ofnucleation/decay events ( ≈ nπ ). B. Spatially localized periodic orbits
We now investigate the existence of periodic orbits like the state exemplified in fig. 5(b).For ρ = 0 .
1, we do a parameter scan varying the mean forcing amplitude r − ≤ r ≤ r + in steps of ∆ r = 10 − and the oscillation period 10 ≤ T ≤
400 in steps of ∆ T = 1. Ateach point a simulation is run to calculate (cid:104) V f (cid:105) initialized with a steady state localizedsolution at r ≡ r . In most cases the simulations were run for 2000 units of time (4000units were necessary for the longer oscillation periods). The results are shown in fig. 6. Theregion where (cid:104) V f (cid:105) = 0 is labeled P O and corresponds to parameter values at which periodicorbits are found. For short periods T for which there is insufficient time for nucleation orannihilation within a cycle, the region of periodic orbits spans nearly the whole pinning region(bottom part of fig. 6). With increasing T the range of existence of periodic orbits narrows aspredicted by the theory in section III for fast oscillations but does not do so monotonically.The figure reveals that sweet spots where the range is larger than in the pinched regions aboveand below occur at regular intervals of the forcing cycle period. For ρ = 0 .
1, the pinched17 T P O −0.33 −0.31 −0.29 −0.2710100200300400
FIG. 6. Diagram showing the region of existence of periodic orbits
P O (shaded region) foroscillation amplitude ρ = 0 .
1. The simulations are initialized with a L solution that is stablefor constant forcing r . The range of r shown corresponds to the pinning interval ( r − , r + ) in thetime-independent case. The three dots (lower right corner, 2nd and 7th sweet spots from bottom)indicate parameter values for the periodic orbits shown in fig. 7. regions are separated by ∆ T ≈
43. The region of existence of periodic orbits is asymmetricowing to major differences in the depinning dynamics in regimes D − and D + . Moreover,the region slants to higher values of the forcing as the period T increases, a property thatis related to the additional time spent in regime A − during the decay phase. Region P O eventually asymptotes to r ≈ − . T → ∞ , the threshold for entering regime A − ,where amplitude decay takes over from depinning as the leading mode of decay. In contrastto the high frequency case, here the Maxwell point determined from the time-averaged freeenergy moves to lower values of r and is no longer a good predictor of the region of periodicorbits.Figure 7 shows three different stable periodic orbits from different sweet spots, corre-sponding to T = 10, 100 and 300. The left panels indicate that these solutions are convergedto machine precision and do not seem to suffer from slow instabilities, while the remaining18 −30 forcing cycle || ∆ u || f A sn A − A + π π (a) T = 10, r = − . −30 forcing cycle || ∆ u || f A sn A − A + π π (b) T = 100, r = − . −30 forcing cycle || ∆ u || f A sn A − A + π π (c) T = 300, r = − . FIG. 7. Periodic orbits for three parameter combinations: The normalized L norm of the differenceof two solutions exactly one period apart showing convergence to machine precision (left panels),space-time plots of the corresponding converged solution over one cycle period (middle panels), andthe ( A, f ) trajectory of the converged solution (right panels). panels provide insight into the balance between growth and decay over the course of thecycle period. Figure 7(a), for T = 10, shows a periodic pulsation in amplitude but no frontmotion. Figure 7(b), for T = 100, reveals a periodic orbit characterized by both amplitudeand front oscillation as does fig. 7(c) for T = 300. The T = 100 orbit, which is located in the second fully formed sweet spot from the bottom of fig. 6, undergoes two nucleation eventsfollowed by two decay events during the course of each forcing cycle. The T = 300 exampleis from the 7th sweet spot and undergoes 7 nucleation/decay events per cycle; the exampleshows that the nucleations occuring during the growth phase of the forcing between t ≈ t ≈
150 are significantly slower than the decay between t ≈
200 and t ≈ A − , something that is only possible because ofthe short amount of time spent in this regime.We also examined the dependence of the region P O on the amplitude of oscillation, ρ . The results for T ≤
200 (fig. 8) show that as ρ increases, the sweet spots span anincreasingly smaller interval in the period T and thus a larger variety of periodic orbits canbe observed within a given range of T as ρ increases. In addition, the whole sweet spotstructure asymptotes more quickly towards r = − . ρ increases. r T −0.30 −0.2810100200 (a) ρ = 0 . r T −0.30 −0.2810100200 (b) ρ = 0 . r T −0.30 −0.2810100200 (c) ρ = 0 . FIG. 8. Region of existence of periodic orbits when (a) ρ = 0 .
06, (b) ρ = 0 .
08, (c) ρ = 0 .
1, usingthe same color code as in fig. 6.
C. Structures undergoing net growth or decay
The existence of the periodic orbits discovered above is closely related to the dynamicsof the fronts connecting the localized patterned state to the background state, suggestingthat we can use the net displacement (cid:104) ∆ f (cid:105) of the fronts within a forcing cycle to classify thegrowing/decaying solutions. We therefore calculated (cid:104) ∆ f (cid:105) on the same grid as that usedto find the periodic orbits P O , starting from narrower localized states to the right of
P O and broader localized states to the left of
P O , all stable. In some cases (e.g. for
T >
P O calculations was necessaryto capture enough oscillations.The results are summarized in fig. 9. The different colored regions are determined bythe conditions ( n − . π < (cid:104) ∆ f (cid:105) < ( n + 0 . π , n = ± , ± , . . . and represent regionswhere regular behavior is observed. The zones between these regions (shown in gray) are“transition zones” that will be discussed below. The figure shows that the region of existence20 T P O −0.33 −0.31 −0.29 −0.2710100200300400
FIG. 9. The number of spatial periods added/lost per cycle for an oscillation amplitude ρ = 0 .
1. Thesimulations were initialized with L localized solutions that are stable at r in the time-independentsystem. The central black region corresponds to the P O region (cf. the shaded region in fig. 6). Thelight blue region to the left corresponds to decay by one wavelength on each side of the localizedstate per cycle, the next to decay by two wavelengths per cycle etc. The regions to the right of
P O correspond instead to net growth by one wavelength, two wavelengths etc. on each side of thelocalized state per cycle. The large white region to the left indicates the location of decay to thetrivial state within one cycle period. Transition zones where irregular behavior is observed areshown in gray. The dots indicate the location in parameter space of the solutions plotted in fig. 5while the horizontal line refers to a region that is studied in fig. 12. of the periodic orbits is surrounded by regions of decay (to the left) and growth (to theright). Beginning in the periodic orbit region
P O and moving to the right (increasing r ),the first region encountered ( O +1 ) corresponds to growth by one wavelength on either sideof the pattern per cycle. The next region ( O +2 ) corresponds to growth by two wavelengthson either side, and so on for the subsequent regions which we refer to as O n , where n isa positive integer. The regions to the left of P O correspond to decay instead of growth.The closest region to
P O , O − , exhibits one wavelength decay on either side of the pattern21er cycle and so on for O n , n < −
1. Each of these regions is separated from its neighborby a transition zone where irregular dynamics are observed and displays the same sweetspot–pinched structure as the
P O region: the regions expand and contract successively as T increases. Some insight into this structure can be gained by looking at the number q ( m ) ofwavelengths gained (lost) on each side of a localized structure during an excursion of thetrajectory into regime D + ( D − ). A sketch of the corresponding results in fig. 10 shows thatareas corresponding to the gain (loss) of a fixed number of wavelengths during a D + ( D − )excursion form bands, and that the intersections of these bands define subregions labeled − m O + qn , where n = q − m , corresponding to net gain or loss of n wavelengths per cycleresulting from the annihilation of m wavelengths followed by the nucleation of q wavelengths(fig. 10(c)). This procedure allows us to assign a unique label to each subregion in theparameter plane (excluding the transition zones in between). Spending more time or goingfarther into D + ( D − ) will result in more nucleations (annihilations) over a forcing cyclebecause more time is spent outside of the pinning region. This explains the evolution of the − m O + qn structure as r increases: the time spent in D + increases and the time spent in D − decreases. Similarly, as the period T of the forcing increases, more time is spent in both D + and D − , resulting in an increase of both q and m . This translates into larger oscillations inthe location of the fronts of the localized structures for longer periods. Finally we can gainintuition about the cliff beyond which localized solutions collapse to the trivial state withina single forcing cycle by considering the length of time spent in the regime A − . This regionis characterized by the time required for solutions to decay to the homogeneous state. Asthe cycle period increases, the center of oscillation that allows the system to just reach thisthreshold time is pushed further to the right. So the edge of the cliff moves to increasingvalues of r as the cycle period T increases.The transition zones narrow as the period T increases, and a closer look reveals a complexstructure resembling a devil’s staircase, a characteristic of mode-locking. Figure 11 shows (cid:104) f (cid:105) within the transition zone between P O and O +1 at T = 80 calculated on a domain of160 spatial periods using a grid of r values with spacing ∆ r = 10 − . The results revealthe presence of increasingly thin regions in which n wavelengths are gained/lost from eitherside of the localized structure within N cycles. These regions thus correspond to fractionalgrowth/decay of the solutions suggesting a complex structure on all scales. Whether there areregions of nonperiodic dynamics corresponding to irrational numbers cannot be determined22 +4 +5 +6 +7 +8 +9 (a) -0 -2 -1 -4 -3 -5 -6 (b)(c) FIG. 10. Sketch of the sweet spot classification scheme. The lines indicate transitions between thenumber of wavelengths gained in (a) and lost in (b) on each side of the localized pattern during onecycle period. These lines are superimposed over the data from numerical simulations in (c) as ameans to classify the regions of growth and decay. through simulations, but the asymptotic results of section IVE seem to indicate that theyform a dense subset of the transition zone.
D. Amplitude decay
In addition to the depinning-like dynamics observed in the colored regions outside of
P O in fig. 9, amplitude decay occurs in the white region. In this region the initial localizedsolution collapses to the trivial state within a single forcing cycle. The boundary of this23 ∆ f i r P O O +1 FIG. 11. The quantity (cid:104) ∆ f (cid:105) in the transition zone between P O and O +1 at T = 80 exhibits astructure characteristic of a devil’s staircase. region is formed by the accumulation of the depinning bands identified in fig. 10(b). −0.293 −0.290 −0.287 −0.284−6−30 N r P O (a) −7 −5 −3 −110 −5 −3 N ∆ r (b) FIG. 12. (a) The quantity N ≡ (cid:104) ∆ f (cid:105) / π , representing the shift in the front location averaged overa cycle period, as a function of r in regime D − when T = 100 (fig. 9, horizontal line). (b) Length ofthe plateaus (green crosses) and of the transition regions between them (red circles) as determinedfrom (a), shown in a semilog plot, together with linear approximations to the data (straight lines,given in the text). The plateaus are plotted at integer values of N while the transition zonesbetween plateaus N = n and N = n + 1 are taken to correspond to N = n + 0 . We look more carefully at the accumulation point of the decay bands in fig. 12. The plateauscorrespond to the loss of integer numbers of wavelengths per forcing cycle. Figure 12(a) showsthat the width of the plateaus as well as of the transition zones between them decreases asone approaches the accumulation point, while fig. 12(b) shows the width ∆ r of the plateausand the transition zones as a function of N , the number of wavelengths lost per cycle ( N < | N | r P ( N P ) = 0 . e . N P , (26)∆ r T ( N T ) = 0 . e . N T , (27)where ∆ r P (resp. ∆ r T ) denotes the width of the plateau corresponding to the loss of N P wavelengths per cycle (resp. width of the transition zone between the pair of closest integersto N T ). To obtain these fits we used all the data in fig. 12(b) on the plateau widths but onlythe transition zones between N = − . N = − .
5. To consider even smaller values of N would have required considerably more numerical effort without improving substantiallythe accuracy of the fit, while values of N closer to 0 lead to departures from the asymptoticregime. Both formulas show similar exponential decrease, thereby confirming the presence ofan abrupt “cliff” at the accumulation point (fig. 12(a)). Furthermore, we see that the widthof the transition zones tends to about 1 / | N | increases. E. Asymptotic theory: small oscillations
To understand the structure of the parameter plane in fig. 9 we need to understandthe process of depinning in the time-dependent system. For this purpose we will considerparameter excursions that take the system outside of the pinning region long enough for anucleation or annihilation event to occur. We therefore suppose that r → r + + (cid:15) ( δ + ρ sin (cid:15)ωt )for which the oscillation period is of the same order as the nucleation time. An analogouscalculation near r − would produce similar results. In this regime the problem is governed bythe equation u t = (cid:0) r + + (cid:15) ( δ + ρ sin (cid:15)ωt ) (cid:1) u − (cid:0) ∂ x (cid:1) u + bu − u . (28)Since the dynamics takes place on an O ( (cid:15) − ) timescale we define the slow timescale τ = (cid:15)t and write ∂ t → (cid:15)∂ τ . We look for a solution in the form u = u + (cid:15)u + (cid:15) u + . . . , obtaining,at leading order, r + u − (cid:0) ∂ x (cid:1) u + bu − u = 0 . (29)As a result we pick u to be a localized solution at a saddle-node bifurcation of the snakingbranch in the time-independent case. In this case u is stationary but only marginally stable.At O ( (cid:15) ), we obtain ∂ τ u = r + u − (cid:0) ∂ x (cid:1) u + 2 bu u − u u . (30)25ince u is stationary, u must be of the form of a zero eigenvector of the Swift–Hohenbergequation linearized about the saddle-node solution. The relevant eigenvector v + correspondsto wavelength addition and is symmetric with respect to x → − x . Since we focus on statesthat do not drift we can ignore the marginal but antisymmetric eigenvectors correspondingto translation and phase. Thus u = a ( τ ) v + .To determine the amplitude a we must go to O ( (cid:15) ). At this order, the equation is ∂ τ u = r + u − (cid:0) ∂ x (cid:1) u + 2 bu u − u u + ( δ + ρ sin ωτ ) u + bu − u u . (31)The solvability condition for u is [7] α ˙ a = α ( δ + ρ sin ωτ ) + α a , (32)with the coefficients α j calculated from the integrals defined in eq. (5).We can turn this equation into a Mathieu equation using the Riccati transformation a = − α ˙ b/α b , obtaining ¨ b = − Ω ( δ + ρ sin ωτ ) b, δ > , (33)where Ω ≡ α α /α > α j evaluated at r + (cf. section IIB). Thus Ω + ≈ . δ from theboundary, r = r − + δ : ¨ b = Ω − ( δ + ρ sin ωτ ) b, δ < , (34)where we have set a = α ˙ b/α b . The integrals α j are now evaluated at r = r − andΩ − = (cid:112) − α α /α ≈ . ω and the characteristicdepinning frequency √ δ Ω + . The properties of eq. (33) are summarized in the standardstability diagram for the Mathieu equation [41] shown in fig. 13 in terms of the scaled distancefrom the boundary of the pinning region δ/ρ and the scaled oscillation period Ω + √ ρT /π .The shaded zones indicate that the solutions of (33) are bounded for all time, while thesolutions are unbounded in the white bands. In terms of the amplitude a in eq. (32), theshaded areas correspond to transition zones where a non-integer number of nucleation eventsoccurs during each cycle of the forcing. In fact, nonperiodic dynamics occur for irrational26alues of the associated Mathieu characteristic exponent within these zones. The first whiteband on the far left corresponds to stable periodic orbits that do not undergo nucleation.The state undergoes one nucleation per oscillation in the white band immediately to theright, and the number of nucleations per oscillation increases by integer values within eachsubsequent white band. - - ∆ (cid:144) Ρ H W Ρ (cid:144) Π L T (a) Π Π- ΩΤ a (b) Π Π- ΩΤ a (c) Π Π- ΩΤ a (d) FIG. 13. (a) The stability diagram for eq. (33). The white bands correspond, from left to right, tostates that undergo exactly 0,1,2,3,... nucleation events per forcing cycle. A non-integer numberof nucleations per cycle occurs in the gray transition zones in between. (b) Sample solution a ( τ )within the 0 region when δ = − ρ/ T = 2 π/ Ω + √ ρ . (c) Sample solution a ( τ ) within the +2 regionwhen δ = 0, T = 6 π/ Ω + √ ρ . (d) Sample solution a ( τ ) within the transition region between regions+3 and +4 when δ = ρ , T = 4 π/ Ω + √ ρ . Nucleation events correspond to divergences in a ( τ ). We remark (cf. section IIB) that care must be taken in interpreting the solutions toeqs. (33) and (34) since the zeros of b ( τ ) correspond to solutions of eq. (32) that diverge to ±∞ .During this process higher order nonlinearities enter eq. (32) with the result that the Riccatitransformation no longer yields a linear equation. Thus the solutions of eqs. (33) and (34)27n fact fail to describe the depinning process near the zeros of b ( τ ), and the correspondingsolution a ( τ ) is determined by “gluing” together a series of individual nucleation events.However, as suggested by the description in eqs. (33) and (34), the resulting nucleationprocess is indeed periodic, albeit in the frame moving with the front x = f ( τ ). F. Asymptotic theory: large oscillations
With intuition gained from the small amplitude theory we now consider the case of largeparameter oscillations that take the system just outside of the pinning region, but with along enough period that there is time for depinning to occur. Because the system spendsonly a small fraction of the forcing cycle outside of P ± , we require the forcing cycle period tobe yet longer, T = O ( (cid:15) − ), in order that depinning takes place. We therefore define the slowtimescale T = (cid:15) t and write the forcing parameter in the form r = r c + (cid:15) r + ( p + (cid:15) δ ) sin (cid:0) (cid:15) ωt (cid:1) , (35)where r c ≡ ( r + + r − ) / p ≡ ( r + − r − ) / r is included for greater generality, as is illustrated in the schematic infig. 14(a). A periodic orbit from a simulation, colored according to the value of the forcingparameter r is also shown.We anticipate that in the above setup nucleation will occur on the faster timescale τ = O ( (cid:15) − ) and so look for solutions in the form u = u + (cid:15)u + (cid:15) u + ... , where u j ≡ u j ( x, τ, T ).Writing ∂ t = (cid:15)∂ τ + (cid:15) ∂ T , we obtain at leading order (cid:2) r c + p sin ( ω T ) − (1 + ∂ x ) (cid:3) u + bu − u = 0 . (36)Thus we can choose u ( x, T ) to be a stable localized solution of the time-independent Swift–Hohenberg equation within the pinning region, with T determining the value of the forcingparameter within this region. The solutions follow the corresponding segment of the L branch (fig. 14(b)) as long as π ( n − / < ω T < π ( n + 1 / n . As we willsee, special care must be taken near the extrema of the forcing cycle when the system leavesthe pinning region and the dynamics take place on the faster timescale τ .28 - r c r + p r (( r p (a) f A − A M A + π π (b) FIG. 14. (a) A schematic of the forcing function r ( t ) used in the asymptotic theory. (b) A periodicorbit with ρ = p + 10 − , T = 5000, and r = − . r :purple corresponds to P ± , orange to D + , and blue to D − . The O ( (cid:15) ) correction reads (cid:2) r c + p sin( ω T ) − (1 + ∂ x ) + 2 bu − u (cid:3) u = 0 . (37)When ω T = (2 n + 1 / π (resp. (2 n + 3 / π ), the quantity u solves the linearized Swift–Hohenberg equation at r + (resp. r − ) with symmetric solution v + (resp. v − ), i.e. thedepinning mode responsible for growth (resp. decay) of the localized pattern. In contrast,when π ( n − / < ω T < π ( n + 1 / r + we take the slow time to be ω T = π/ (cid:15)θ ; asimilar procedure can be carried out near r − by taking ω T = 3 π/ (cid:15)θ , and subsequentcycles of the forcing can be handled in the same way. The time derivative now becomes ∂ t = (cid:15)ω∂ θ and the O ( (cid:15) ) equation (37) becomes (cid:2) r + − (1 + ∂ x ) + 2 bu − u (cid:3) u = ω∂ θ u . (38)29ince u is the marginally stable localized solution at r + it follows that ∂ θ u = 0 and hencethat u = a ( θ ) v + ( x ).At O ( (cid:15) ) we obtain (cid:2) r + − (1 + ∂ x ) + 2 bu − u (cid:3) u = ω∂ θ u − ( b − u ) u − ( r + δ − pθ ) u , (39)for which the solvability condition is α ωa (cid:48) = α ( r + δ − pθ ) + α a , (40)where the prime denotes the θ derivative and the coefficients α j are determined by theintegrals in eq. (5).Using the transformation a = − α ωb (cid:48) /α b , we obtain a linear oscillator problem with atime-dependent frequency, b (cid:48)(cid:48) = − p Ω ω (cid:0) θ − θ (cid:1) b, (41)where θ = 2( r + δ ) /p and, as before, Ω = α α /α . The system exits P + when r + δ > − θ + , θ + ] corresponds to the time interval spent in D + . We now use a matchingprocedure to connect this solution to the case when ω T (cid:54) = π/
2, noting that u ( x, T ) → ω T → π/
2, so that the solution remains stable as it approaches the boundary of the pinningregion. Since the leading order solution for large | θ | is given by a ( θ ) ≈ (cid:112) pα / α | θ | werequire that a ( θ ) → (cid:112) pα / α θ < θ → −∞ . The solution of eq. (40) satisfying thisrequirement can be written in terms of parabolic cylinder functions [42], b = D ν ( z ), where ν = √ p Ω + θ √ ω − , z = − (2 p ) / (cid:114) Ω + ω θ. (42)Each zero z = z of b = D ν ( z ) corresponds to a nucleation event, since a diverges to ±∞ as z → z . To determine the outcome of such an event we must consider the limit as θ → ∞ tomatch the ˙ r > r < π/ < ω T < π/
2. The parabolic cylinder functionbehaves like D ν ( z ) → √ π Γ[ − ν ] | z | − ( ν +1) e z / (43)for real z → −∞ , where Γ[ ν ] is now the Gamma function. Examining the sign of Γ[ ν ] showsthat as long as ν is not a positive integer, a → (cid:112) pα / α θ > ν is apositive integer correspond to the solution landing exactly on an unstable solution.30e can get a very simple expression for the number of nucleation events that occur near r + by counting the number of real zeros of D ν ( z ) for a given ν . We start by noting that for z → ∞ ( θ → −∞ ), D ν ( z ) approaches zero from above. When ν <
0, there are no real zerosand D ν ( z ) > z . Equation (43) shows that in the limit that z → −∞ , the sign of D ν ( z ) depends on the sign of Γ[ − ν ]. For ν <
0, the sign is positive and D ν ( z ) → + ∞ as z → −∞ without crossing zero. At ν = 0, Γ[ − ν ] = ∞ and D ν ( z ) → < ν <
1, there will be one zero crossing as D ν ( z ) → −∞ as z → −∞ . The number of zeros continues to increase by one each time there is a sign changein Γ[ − ν ] so that for n − < ν < n , there will be n > D ν ( z ). Therefore there willbe n + nucleation events if n + − < Ω + T π √ p ( r + ρ − r + ) < n + + , (44)where we have re-expressed the condition in terms of the amplitude ρ ≡ p + (cid:15) δ , offset r ≡ r c + (cid:15) r , and the period T ≡ π/(cid:15) ω . A similar relation applied for the number ofannihilations n − < n − − < Ω − T π √ p ( r − ρ − r − ) < n − + . (45)The above conditions also reveal the presence of bifurcation delay, as expected of anonautomous bifurcation problem. This delay manifests itself in the shift of the criticalvalue r + δ = 0 for the presence of a fold to the threshold value determined by ν = 0, viz., r + δ = ω/ √ p Ω + : the system enters D + by as much as ω/ √ p Ω + without triggering anucleation event. Figure 15(b,c) shows the amplitude a as a function of the scaled time(2 p ) / Ω + θ/ω just before and after this threshold. The discontinuous jump in fig. 15(c)represents a nucleation event and is obtained by gluing together two separate asymptoticcalculations near different, but adjacent saddle-nodes on the same snaking branch. The same“inertial” effect is observable even when the system does not leave P + , i.e., δ + r <
0. Usingthe property D ν ( z ) → √ ν π (cid:32) − ν ) / − √ z Γ[ − ν/ − (1 + 2 ν ) z − ν ) / (cid:33) , z → , (46)we see that even when the system just barely reaches the boundary of the pinning region, r + δ = 0, the perturbation a ( θ ) remains finite (fig. 15(a)). Indeed the minimum value of31 a r a (a) ν = − . Θ a r a (b) ν = − . Θ a r a (c) ν = 0 . FIG. 15. A plot of the amplitude a ( θ ) of the O ( (cid:15) ) correction u ( x ) to the solution that is marginallystable at r = r + as a function of a scaled O ( (cid:15) ) time near the boundary of the pinning region fordifferent values of ν = √ p Ω + θ √ ω ( r − ρ − r − ) − . The time θ = 0 corresponds to the peak of theforcing cycle where “inertial” effects are expected. The thin solid (dashed) line shows the amplitude a for the stable (unstable) localized solution as functions of r but replotted in terms of the time θ .Below each frame is a schematic representation of the trajectory of the amplitude a as a functionof the forcing parameter r . The stable (unstable) steady state branches of the constant forcingcase are shown in solid (dashed) lines for reference. (a) ν = − .
5: the system does not leave thepinning region, but there are still deviations from the stable state. (b) ν = − .
1: the system exitsthe pinning region, but not far enough for nucleations to occur. (c) ν = 0 .
1: the system penetratesinto D + past the threshold for a nucleation to occur (represented by a discontinuous jump). a ( θ ) occurs for θ > θ = 0. In fact a ( θ ) can be calculated explicitly in terms ofparabolic cylinder functions using the relation dD ν ( z ) dz = zD ν ( z ) − D ν +1 ( z ) (fig. 15).In the parameter regime analyzed here, the system tracks a given stable localized statewith a delay in r of up to O ( (cid:15) ) for most of the forcing cycle. All of the interesting dynamicsoccur within a small time interval when the system visits the vicinity of the boundary of thepinning region. If it ventures far enough outside of this boundary, nucleation/annihilationevents begin to take place after a delay. Once the system reenters the pinning region, thesystem settles on the nearest stable but longer/shorter localized structure. The settlingprocess also happens within the vicinity of the boundary of the pinning region and there may32r may not be an additional nucleation/annihilation event during this settling, depending onwhere in the process the system was upon re-entering the pinning region P ± .To understand the structure of growing, steady-state, and decaying solutions in thislimit, we need only compare the growth from depinning that occurs near ω T = π/ ω T = 3 π/
2. The formation of the pinched zonesand the sweet spot structure of the stationary solutions can be predicted by balancing thegrowth on the right of the pinning region to the decay on the left as we shall now see. Theresulting prediction is compared with numerical simulation in fig. 16 for ρ = p + 10 − , andshows excellent quantitative agreement. Specifically, the colored regions are determined bynumerical simulation and these match the coding scheme of fig. 9 while the red and blue linesare predictions for the transitions between the various regions of the classification schemedetailed in fig. 10. Note, however, that the values of r in fig. 16 span only about 1 / ρ . P Or T −0.3000 −0.2995 −0.2990 −0.29855010002000300040005000 FIG. 16. A comparison of the asymptotic theory (red/blue lines) with numerical simulations (colors)for ρ = p + 10 − ≈ .
04. The dark region corresponds to the region
P O . The red (blue) linesindicate transitions in the number of nucleations (decays) that occur during one forcing cycle. r , ρ ) parameter plane for T = 100 and T = 200, along with anextension of the predictions from the above asymptotic theory (eqs. (44) and (45)). Theextension has been computed by replacing √ p in the denominator of the expressions by √ ρ as a means for correcting for the cases when ρ (cid:54)≈ p . The modified theory is able to accuratelypredict the location of the transition between the zero and ± r ρ −0.33 −0.31 −0.29 −0.2700.050.10 (a) T = 100 r ρ −0.33 −0.31 −0.29 −0.2700.050.10 (b) T = 200 FIG. 17. The number of spatial periods gained/lost in one forcing cycle when (a) T = 100 and (b) T = 200. All simulations were initialized with stable L solutions at the corresponding r in theconstant forcing case. The purple region labeled P O indicates the location of periodic orbits andcorresponds to the blue region in fig. 6. The light blue region immediately to the left indicatesdecay by one wavelength on each side of the localized structure per forcing cycle, the next region tothe loss of two wavelengths per cycle, and so on. The solution grows by one wavelength on eachside of the localized structure per forcing cycle in the region immediately to the right of the darkregion, and so on. The white region to the left indicates solutions that collapse to the trivial statewithin one cycle.
V. THE LOW-FREQUENCY LIMIT: ADIABATIC THEORY
In this section we consider the remaining case, that of low-frequency forcing. In thisregime we may neglect inertial effects that can cause delays in the onset of depinning or34llow for the completion of depinning events within the pinning region. We note that, byapplying the technique of matched asymptotics (see, for example, [43]), we can estimate thedepinning delay to be ∼ | dr/dt | − / where the derivative is evaluated at r ± . A. Sweet spot structure
Using the adiabatic approximation described above, the number n ± of nucleation/annihilationevents over the course of a forcing cycle can be estimated from the expression n ± = ± (cid:90) T ± dtT dpn ± ( t ) , (47)where T ± is the time spent outside of the pinning region and T dpn ± ( t ) is the time betweennucleation/annihilation events of the constant forcing problem with parameter r ( t ). Thesuper/subscript + (resp. − ) refers to regime D + (resp. D − ). We assume that the dynamicswithin the pinning region allow the system to either complete the nucleation process (corre-sponding to rounding n ± up), or settle back down to the state already reached (correspondingto rounding n ± down). We also suppose that the threshold between completing a nucleationevent or settling back corresponds to n ± + 1 / n ± ], to denote thenearest integer. We recall that leading order asymptotics near the edge of the pinning regionpredict that ( T dpn ± ) − = Ω ± δ / ± /π [7] and use this prediction together with the assumption r ( t ) = r + ρ sin 2 πt/T to obtain n ± = ± √ ρ Ω ± Tπ (cid:2) K (cid:0) − η ± (cid:1) − η ± E (cid:0) − η ± (cid:1)(cid:3) , (48)where η ± = | r − r ± | /ρ < K (m) = π/ (cid:90) (cid:112) − m sin θ dθ E (m) = π/ (cid:90) (cid:112) − m sin θ dθ. (49)are the complete elliptic integrals of the first and second kind [42].The predictions of the adiabatic theory in eq. (48) are shown in fig. 18 for ρ = p + 10 − ≈ .
04 and ρ = 0 .
1. The red and blue lines indicate transitions between adjacent values of [ n + ]and [ n − ], respectively. The plot in fig. 18(a) is colored according to the simulation resultsto emphasize the quantitative accuracy of the adiabatic prediction for ρ = p + 10 − . The35ccuracy of these predictions diminishes with increasing ρ as shown in fig. 18(b) for ρ = 0 . P O based on the theory (eq. (48)), i.e.,
P O is the region where [ n + ] + [ n − ] = 0. The predictions r T −0.3000 −0.2995 −0.2990 −0.29855010002000300040005000 (a) ρ = p + 10 − , theory and simulation (b) ρ = 0 .
1, theory only
FIG. 18. Predictions from adiabatic theory using the asymptotic approximation eq. (48) for thedepinning time when (a) ρ = p + 10 − ≈ .
04 and (b) ρ = 0 .
1. The colors in (a) refer to thenumerical simulation results for the locations of
P O (dark), O + n (alternating yellow and orange),and O − n (alternating shades of blue). The dark region in (b) is the adiabatic theory prediction ofthe P O region. for the ρ = 0 . P O does not slant as in the simulations (fig. 9). The qualitative disagreement occurs becauseof a breakdown of the asymptotic prediction for T dpn ± when the system enters too far intoregions D ± . In addition to the quantitative disagreement of the depinning times, the theoryomits the amplitude mode that destroys the localized states in A − . We can account for (i)by making use of numerical fits in place of the asymptotic theory for T dpn ± as described infig. 4 and can also extend the theory to include predictions about the cliff mentioned in(ii) as described in the next section. The theory cannot, however, account for (iii) as theslanting is a result of the coupling between the amplitude and depinning modes which wehave neglected. 36 . The cliff FIG. 19. Adiabatic prediction in the ( r , T ) plane of the decay vs nucleation dynamics of a spatiallylocalized initial condition of eq. (1) with time-periodic forcing. Positive (resp. negative) numbersrepresent [ n + ] (resp. [ n − ]), the change in the number of wavelengths due to nucleation (resp.annihilation) events during one cycle. The results are obtained using the relation (48) with ρ = 0 . b = 1 .
8. The figure is plotted over the same r interval as fig. 9. We can approximate the dynamics of the overall amplitude decay of a localized state in A − by computing the time T col ( r ), r < r sn , for a solution initialized with the periodic stateat r sn to decay to the trivial state. This calculation mirrors the asymptotic calculation for T dpn ± . For r = r sn + δ col , | δ col | (cid:28)
1, we find( T col ) − ≈ Ω col π | δ col | / . (50)Substituting T col into eq. (47) in place of T dpn ± leads to an equation analogous to eq. (48)for n col . Since n col is at most one we assume that the threshold for the solution to decayirrevocably is n col = 1 /
2. This procedure yields a prediction for the location of the cliff inparameter space. For improved numerical accuracy for forcing cycles that penetrate far into A − , a numerical fit is useful (fig. 4). Figure 19 reveals the dramatic improvement in the37 r , T ) phase diagram that results from this procedure applied to both T dpn ± and T col (boldblack line). The hybrid adiabatic theory augmented with numerical fits from simulationsof the system under constant forcing is in remarkably good agreement with the simulationseven for ρ = 0 .
1. The predicted extent in T of the sweet spots is (cid:52) T ≈
45, independently of T , which is within 5% of the value computed from simulations, viz., (cid:52) T ≈
43. In addition,the predicted period at which the cliff occurs for a given value of of r deviates from thevalue computed from simulations by (cid:52) T (cid:46)
10, with the maximum deviation occurring forperiods below T = 100. As expected, the agreement improves for larger periods and awayfrom the cliff: the top of the subregion − O +10+7 is predicted to be at ( r , T ) ≈ ( − . , r , T ) ≈ ( − . , VI. CONCLUSION
We have considered the effects of parametric time-periodic forcing on the dynamicsof localized structures in the Swift–Hohenberg equation with competing quadratic-cubicnonlinearities. In the high frequency limit averaging theory yields an averaged system that isalso of Swift–Hohenberg-type. When oscillations are large enough, the time-varying forcingaffects the averaged dynamics by modifying the coefficients of the nonlinear terms, therebyreducing the region of existence of spatially localized states (the pinning region in the caseof constant forcing) and displacing it to larger values of r . For intermediate frequencies,the dynamics become more complex owing to depinning of the fronts bounding the localizedstructure over a significant fraction of the forcing cycle, resulting in breathing localizedstructures exhibiting behavior analogous to pinning, depinning, and amplitude collapsefamiliar from the constant forcing case. Of particular significance is the observation of a newresonance phenomenon between the forcing period and the time required to nucleate a newwavelength of the pattern. The presence of this resonance is responsible for the complexstructure of the parameter space, which breaks up into regions labeled by a pair of integers( m, q ) denoting the number of wavelengths lost ( m ) per cycle and the number gained ( q ).When n ≡ q − m = 0 the resulting state is periodic in time, and corresponds to a state thaton average neither expands nor shrinks. We have described the resulting structure of theparameter space in terms of sweet-spots favoring the existence of such “pinned” states andpinched zones where the resonance was destructive and periodic localized structures absent.38e found that these properties could be understood on the basis of appropriate asymptotics,valid either when the forcing cycle did not penetrate far into the depinning regions, or forlow frequency forcing. In both cases we showed that the number of nucleation/annihilationevents can be computed by adapting existing theory of the depinning process, and used theseresults to partition the parameter space. A similar approach was successful in obtaining theaccumulation point of the decay regions beyond which all initial conditions collapse to thetrivial state within one forcing cycle. Our calculations suggest that this accumulation isexponential and involves regions of frequency locking corresponding to all rational numbers.We found that asymptotic theory provided an excellent qualitative description of theresonance phenomenon, and moreover that quantitative agreement could often be obtainedby augmenting the leading order nucleation theory with numerical fits to the nucleation timesadopted from the time-independent case, thereby greatly extending the range of validity ofthe theory.In view of the success of the Swift–Hohenberg equation in modeling localization in agreat variety of systems with bistability between a homogeneous and a patterned state weexpect that the model studied here, eq. (1), captures faithfully the phenomenology arisingfrom a resonance between the forcing period and the nucleation time in systems undergoingtemporary depinning as a result of the forcing. As such we envisage numerous applications ofthe theory presented here to temporally forced systems such as models of vegetation growthin arid regions subject to seasonal forcing [44, 45].In future work it will be of interest to extend the present analysis to systems withtime-periodic forcing that is not purely sinusoidal and in particular to periodic forcingwith asymmetry between the rise and fall phases. In addition, many experimental systemsexhibiting spatially localized states, including binary fluid convection [46] and plane Couetteflow [10], possess an additional midplane reflection symmetry whose effects are well modeledby the Swift–Hohenberg equation with a cubic-quintic nonlinearity [47, 48]. In future work wewill investigate the properties of temporally-forced cubic-quintic Swift–Hohenberg equationwith a view to elucidating the behavior expected when such systems are forced periodically.39 CKNOWLEDGMENTS
This work was supported in part by the National Science Foundation under GrantsDMS-1211953 and CMMI-1232902. [1] A. R. Champneys, “Homoclinic orbits in reversible systems and their applications in mechanics,fluids and optics,” Physica D , 158–186 (1998).[2] E. Knobloch, “Spatially localized structures in dissipative systems: open problems,” Nonlin-earity , T45–T60 (2008).[3] N. J. Balmforth, “Solitary waves and homoclinic orbits,” Annu. Rev. Fluid Mech. , 335–373(1995).[4] P. D. Woods and A. R. Champneys, “Heteroclinic tangles and homoclinic snaking in theunfolding of a degenerate reversible Hamiltonian–Hopf bifurcation,” Physica D , 147–170(1999).[5] Y. Pomeau, “Front motion, metastability and subcritical bifurcations in hydrodynamics,”Physica D , 3–11 (1986).[6] G. W. Hunt, M. A. Peletier, A. R. Champneys, P. D. Woods, M. A. Wadee, C. J. Budd, andG. J. Lord, “Cellular buckling in long structures,” Nonlinear Dynamics , 3–29 (2000).[7] J. Burke and E. Knobloch, “Localized states in the generalized Swift-Hohenberg equation,”Phys. Rev. E , 056211 (2006).[8] J. Burke and E. Knobloch, “Snakes and ladders: localized states in the Swift–Hohenbergequation,” Phys. Lett. A , 681–688 (2007).[9] G. W. Hunt, H. M. Bolt, and J. M. T. Thompson, “Structural localization phenomena andthe dynamical phase-space analogy,” Proc. R. Soc. London A , 245–267 (1989).[10] T. M. Schneider, J. Gibson, and J. Burke, “Snakes and ladders: localized solutions of planeCouette flow,” Phys. Rev. Lett. , 104501 (2010).[11] I. Mercader, O. Batiste, A. Alonso, and E. Knobloch, “Convectons, anticonvectons andmulticonvectons in binary fluid convection,” J. Fluid Mech. , 586–606 (2011).[12] C. Beaume, A. Bergeon, and E. Knobloch, “Convectons and secondary snaking in three-dimensional natural doubly diffusive convection,” Phys. Fluids , 024105 (2013).
13] C. Beaume, E. Knobloch, and A. Bergeon, “Nonsnaking doubly diffusive convectons and thetwist instability,” Phys. Fluids , 114102 (2013).[14] D. Lo Jacono, A. Bergeon, and E. Knobloch, “Three-dimensional binary fluid convection in aporous medium,” J. Fluid Mech. , R2 (2013).[15] C. Beaume, H.-C. Kao, E. Knobloch, and A. Bergeon, “Localized rotating convection withno-slip boundary conditions,” Phys. Fluids , 124105 (2013).[16] B. Sch¨apers, M. Feldmann, T. Ackemann, and W. Lange, “Interaction of localized structuresin an optical pattern-forming system,” Phys. Rev. Lett. , 748–751 (2000).[17] A. Prigent, G. Gr´egoire, H. Chat´e, O. Dauchot, and W. van Saarloos, “Large-scale finite-wavelength modulation within turbulent shear flows,” Phys. Rev. Lett. , 014501 (2002).[18] D. Barkley and L. S. Tuckerman, “Computational study of turbulent laminar patterns inCouette flow,” Phys. Rev. Lett. , 014502 (2005).[19] M. G. Clerc, C. Falcon, and E. Tirapegui, “Additive noise induces front propagation,” Phys.Rev. Lett. , 148302 (2005).[20] P. B. Umbanhowar, F. Melo, and H. L. Swinney, “Localized excitations in a vertically vibratedgranular layer,” Nature , 793–796 (1995).[21] J. M. G. Vilar and J. M. Rubi, “Spatiotemporal stochastic resonance in the Swift–Hohenbergequation,” arXiv preprint cond-mat/9702205 (1997).[22] O. Lioubashevski, Y. Hamiel, A. Agnon, Z. Reches, and J. Fineberg, “Oscillons and propagatingsolitary waves in a vertically vibrated colloidal suspension,” Phys. Rev. Lett. , 3190–3193(1999).[23] P. Binder, D. Abraimov, A. V. Ustinov, S. Flach, and Y. Zolotaryuk, “Observation of breathersin Josephson ladders,” Phys. Rev. Lett. , 745–748 (2000).[24] A. M. Rucklidge and M. Silber, “Design of parametrically forced patterns and quasipatterns,”SIAM J. Appl. Math. , 298–347 (2009).[25] J. V. I. Timonen, M. Latikka, L. Leibler, R. H. A. Ras, and O. Ikkala, “Switchable staticand dynamic self-assembly of magnetic droplets on superhydrophobic surfaces,” Science ,253–257 (2013).[26] A. Yochelis, C. Elphick, A. Hagberg, and E. Meron, “Frequency locking in extended systems:The impact of a Turing mode,” Europhys. Lett. , 170–176 (2005).
27] B. Marts, A. Hagberg, E. Meron, and A. L. Lin, “Resonant and nonresonant patterns inforced oscillators,” Chaos , 037113 (2006).[28] J. Buceta, K. Lindenberg, and J. M. R. Parrondo, “Stationary and oscillatory spatial patternsinduced by global periodic switching,” Phys. Rev. Lett. , 024103 (2001).[29] I. Belykh, V. Belykh, R. Jeter, and M. Hasler, “Multistable randomly switching oscillators:The odds of meeting a ghost,” Euro. Phys. J. Special Topics , 2497–2507 (2013).[30] J. Burke, A. Yochelis, and E. Knobloch, “Classification of spatially localized oscillations inperiodically forced dissipative systems,” SIAM J. Appl. Dyn. Sys. , 651–711 (2008).[31] J. A. Sherratt, “An analysis of vegetation stripe formation in semi-arid landscapes,” J. Math.Biol. , 183–197 (2005).[32] M. Tlidi, R. Lefever, and A. Vladimirov, “On vegetation clustering, localized bare soil spotsand fairy circles,” in Dissipative Solitons: From Optics to Biology and Medicine , Lecture Notesin Physics, Vol. 751 (Springer, Berlin, 2008) pp. 1–22.[33] J. A. Sherratt, “Pattern solutions of the Klausmeier model for banded vegetation in semi-aridenvironments I,” Nonlinearity , 2657–2675 (2010).[34] A. Y. Kletter, J. von Hardenberg, and E. Meron, “Ostwald ripening in dryland vegetation,”Comm. Pure Appl. Anal. , 261–273 (2012).[35] E. Meron, “Pattern-formation approach to modelling spatially extended ecosystems,” Ecol.Model. , 70–82 (2012).[36] J. Burke and J. Dawes, “Localized states in an extended Swift–Hohenberg equation,” SIAM J.Appl. Dyn. Sys. , 261–284 (2012).[37] J. Burke, S. M. Houghton, and E. Knobloch, “Swift–Hohenberg equation with broken reflectionsymmetry,” Phys. Rev. E , 036202 (2009).[38] H.-C. Kao, C. Beaume, and E. Knobloch, “Spatial localization in heterogeneous systems,”Phys. Rev. E , 012903 (2014).[39] S. M. Cox and P. C. Matthews, “Exponential time differencing for stiff systems,” J. Comp.Phys. , 430–455 (2002).[40] E. J. Doedel, “Auto: A program for the automatic bifurcation analysis of autonomous systems,”Congr. Numer. , 265–284 (1981).[41] N. W. McLachlan, Theory and application of Mathieu functions (Oxford Press, London, 1951).
42] M. Abramowitz and I. A. Stegun,
Handbook of mathematical functions , Vol. 1 (Dover, NewYork, 1972).[43] J. Hinch,
Perturbation methods (Cambridge University Press, Cambridge, 1991).[44] V. Guttal and C. Jayaprakash, “Self-organization and productivity in semi-arid ecosystems:Implications of seasonality in rainfall,” J. Theor. Biol. , 490–500 (2007).[45] K. Siteur, E. Siero, M. B. Eppinga, J. D. M. Rademacher, A. Doelman, and M. Rietkerk,“Beyond Turing: The response of patterned ecosystems to environmental change,” EcologicalComplexity , 81–96 (2014).[46] O. Batiste, E. Knobloch, A. Alonso, and I. Mercader, “Spatially localized binary-fluidconvection,” J. Fluid Mech. , 149–158 (2006).[47] S. M. Houghton and E. Knobloch, “Swift–Hohenberg equation with broken cubic-quinticnonlinearity,” Phys. Rev. E , 016204 (2011).[48] I. Mercader, O. Batiste, A. Alonso, and E. Knobloch, “Travelling convectons in binary fluidconvection,” J. Fluid Mech. , 240–266 (2013)., 240–266 (2013).