A discrete complex Ginzburg-Landau equation for a hydrodynamic active lattice
AA discrete complex Ginzburg-Landau equation for ahydrodynamic active lattice
Stuart J. Thomson , Matthew Durey , and Rodolfo R. Rosales School of Engineering, Brown University, Providence, RI 02912, USA Department of Mathematics, Massachusetts Institute of Technology,Cambridge, MA 02139, USA
Abstract
A discrete and periodic complex Ginzburg-Landau equation, coupled to a discretemean equation, is systematically derived from a driven and dissipative oscillator model,close to the onset of a supercritical Hopf bifurcation. The oscillator model is inspiredby recent experiments exploring active vibrations of quasi-one-dimensional lattices ofself-propelled millimetric droplets bouncing on a vertically vibrating fluid bath. Oursystematic derivation provides a direct link between the constitutive properties of thelattice system and the coefficients of the resultant amplitude equations, paving theway to compare the emergent nonlinear dynamics—namely discrete bright and darksolitons, breathers, and traveling waves—against experiments. Further, the amplitudeequations allow us to rationalize the successive bifurcations leading to these distinctdynamical states. The framework presented herein is expected to be applicable to awider class of oscillators characterized by the presence of a dynamic coupling potentialbetween particles. More broadly, our results point to deeper connections betweennonlinear oscillators and the physics of active and driven matter. a r X i v : . [ n li n . PS ] O c t Introduction
The celebrated complex Ginzburg-Landau equation (CGLE) [1, 2] is a generic model de-scribing the dynamics of spatially extended, dissipative systems near a Hopf bifurcation.In contrast to the potentially complex, high-dimensional microscopic equations regulating aparticular physical system, amplitude equations [3, 4] such as the CGLE are typically castin terms of only a few macroscopic variables, or order parameters [5, 6, 7, 8, 9]. In gen-eral, the form of such effective models may be posited on phenomenological grounds, theirstructure determined through a combination of linear stability and symmetry arguments[3]. This universal approach can, however, obfuscate the connection between the coeffi-cients of the amplitude equation and the physical parameters of the system under study.A more robust approach sacrifices derivational simplicity in favor of obtaining the ampli-tude equation(s) directly from the underlying microscopic equations of the system, typicallycontinuous nonlinear partial differential equations [10, 11, 12, 13]. However, for systemswhich are fundamentally discrete—for example, nonlinear oscillators—amplitude equationsare typically posited as discretized versions of their continuous counterparts, seldom derivedin a systematic manner from the original governing equations [14, 15, 16].We herein present a rigorous framework to systematically derive a discrete and periodiccomplex Ginzburg-Landau equation (dCGLE) for a driven and dissipative nonlinear oscilla-tor, close to the onset of a supercritical Hopf bifurcation. The oscillator model is inspired byrecent experiments exploring the active vibrations of a hydrodynamic lattice of self-propelledmillimetric droplets [17, 18]. The coefficients appearing in our dCGLE are directly relatedto the constitutive properties of the physical lattice system, paving the way to compare thenumerical results of the resultant amplitude equations—namely the emergence of discretebright and dark solitons, breathers, and travelling waves—against experiments. Further, theamplitude equations allow us to rationalize the successive bifurcations leading to these dis-tinct dynamical states. Although we present the case of the hydrodynamic lattice, we proposethat the framework presented herein is applicable to a wider class of oscillators characterizedby the presence of a dynamic coupling potential between particles. On a fundamental level,our results suggest deeper connections between nonlinear and nonlocal oscillators and the2hysics of active and driven matter [8, 19, 20].
This study is motivated by experiments of quasi-one-dimensional lattices of millimetricdroplets, bouncing synchronously and periodically on the surface of a vertically vibratingfluid bath and confined to an annular channel [17]; see Figure 1. (For a broader perspectiveof the physics of bouncing droplets, see [21, 22, 23] and references therein.) Upon successiveimpacts, each droplet excites a field of standing waves whose decay time, T M , increases withthe vertical acceleration of the bath and diverges at the Faraday threshold [24, 25]. Thesuperposition of the wave fields generated by each droplet forms the global lattice wave field,which acts as an inter-droplet potential, mediating the spatiotemporal coupling of the lattice.This wave-mediated coupling represents a distinguishing feature of this new class of coupledoscillator: the waves produced at each droplet impact give rise to an effective self-generated, dynamic coupling potential between droplets, one that evolves continuously with the dropletmotion.For sufficiently weak vibrational forcing, the droplets exhibit stationary bouncing in a cir-cular, equispaced lattice. Above a critical vertical acceleration of the bath (or, alternatively,critical decay time, T M ), the droplets destabilize to small lateral perturbations, oscillatingabout their equilibrium position (Figures 1b and 1c). Physically, these oscillations emergedue to the competition between droplet self-propulsion—arising through the propulsive forceenacted on each droplet by the local slope of the lattice wave field—and wave-mediated, non-local coupling between droplets. Oscillations of the lattice are further offset by dissipativeeffects due to drag. That self-propulsion is achieved and sustained by the continual exchangeof energy of the droplet with its environment—in this case, the vibrating bath—renders thishydrodynamic lattice a novel example of an inertial active system [26, 27, 28, 29].As shown in experiments [17], oscillations of the lattice follow the onset of either a super-critical or subcritical Hopf bifurcation, depending primarily on the proximity of neighboringdroplets. Our focus here is on the supercritical case, for which periodic, small-amplitude,out-of-phase oscillations arise beyond the bifurcation point, initially uniform over all droplets.(When the bifurcation is subcritical, the dynamics are profoundly different: in experiments,3igure 1: (a) Overhead perspective of a chain of 40 equispaced, millimetric droplets of siliconeoil, confined to an annular channel and surrounded by a shallow layer of fluid. The reflectedcolor in the channel emphasises the deformation of the fluid surface as droplets impact thebath and excite subcritical Faraday waves. (b) A subset of droplet polar positions obtainedfrom experiments for a lattice consisting of 20 droplets [17]. Each droplet undergoes out-of-phase oscillations with respect to its neighbor, following a supercritical (Hopf) bifurcation[18]. (c) The instantaneous positions of all 20 droplets in the lattice for the same experimentas (b). The net result of the instability is the out-of-phase oscillations of two decahedralsub-lattices, colored red and blue.the system approaches a distant attractor and the emergent dynamics manifest as a self-sustaining, nonlinear solitary-like wave [17].) The dependence of the form of these bifur-cations on the parameters of the lattice system, and the ensuing dynamics of the uniform,periodic state, was characterized via a weakly nonlinear analysis of a mathematical modeldescribing the droplet lattice [18]. Upon further increase of the vibrational forcing, thisperiodic state can itself destabilize, leading to spatial modulations of the droplet oscilla-tion amplitude, a phenomenon not captured by the analysis presented in [18]. To exploreand rationalize the onset and resultant dynamics of these spatial modulations, we present ageneralized weakly nonlinear theory of the lattice, in the vicinity of the supercritical Hopfbifurcation. We proceed to briefly summarize the results of [18] as they pertain to ourderivation of the governing amplitude equations presented herein.4 .2 Lattice model and linear theory Model —The principal assumption underpinning the hydrodynamic lattice model [18] is thatthe horizontal motion of each of the N droplets in the lattice may be averaged over onebouncing period, which we denote T F . This stroboscopic approximation [30] effectively elim-inates the droplets’ synchronous vertical motion from consideration. To further simplifymatters, we assume that the droplets lie on a circle of constant radius, R , which, in exper-iments, is determined by the inner and outer radii of the annular channel. Combining thismotion with the stroboscopic approximation yields the following equation of motion for thecircumferential position, x n ( t ), of each droplet in the lattice [18]: m ¨ x n + ¯ D ˙ x n = − mg ∂h∂x (cid:12)(cid:12)(cid:12)(cid:12) x = x n . (1a)Dots denote differentiation with respect to time, t , and the space variable, x ∈ [0 , L = 2 πR ],is directed along the circumference of the circle on which the droplets lie. According toequation (1a), the motion of each droplet of mass m is thus governed by a balance betweeninertia, a linear drag with drag coefficient ¯ D , and the time-averaged propulsive wave forceenacted on each droplet by the local slope of the lattice wave field, h ( x, t ). By periodicity, h ( x, t ) = h ( x + L, t ). The remaining parameter, g , is acceleration due to gravity. It is tobe understood that h ( x, t ) is the stroboscopic global lattice wave field—the time-averaged superposition of wave fields generated by each individual droplet in the lattice—projectedonto the circle.A distinguishing feature of the hydrodynamic lattice is that the propulsive wave forceenacted on each droplet depends explicitly on the prior trajectory of every droplet in thelattice [18, 24, 25, 30]. The time-dependent evolution of the lattice wave field, h , may bedescribed by ∂h∂t + 1 T M h = 1 T F N (cid:88) m =1 H ( x − x m ) , (1b)where the wave kernel, H ( x ), is the quasistatic wave field generated by stationary bouncing ofeach individual droplet, time-averaged over T F and projected onto the circle [18]. Promptedby the fundamental aspects of the fluid system [17], our theory requires only that H ( x )be sufficiently smooth and periodic with H ( x ) = H ( x + L ), exhibiting variations over a5haracteristic wavelength, λ , and an exponential spatial decay with lengthscale l d ; see Figure2(b) for a prototypical example. (A candidate wave kernel satisfying these properties ispresented in § h , wave dissipation, and the superposition of waves generated about the instantaneousposition of each droplet in the lattice.The dynamical system (1) is non-dimensionalized via the scalings t = m ¯ D ˆ t = t ˆ t, x = λ ˆ x, h = λ gt ˆ h = h ˆ h, H = h T F t ˆ H . Upon dropping the carets, we arrive at the dimensionless system [18]¨ x n + ˙ x n = − ∂h∂x ( x n , t ) , (2a) ∂h∂t + νh = N (cid:88) m =1 H ( x − x m ) , (2b)where ν = t /T M is the dimensionless dissipation rate of the wave field. Recalling that thedecay time, T M , may be regarded as a proxy for the vertical vibrational acceleration of thebath [24, 25], ν plays the role of a control, or bifurcation, parameter in the dimensionlesssystem (2). While ν is convenient for algebraic manipulations, we will interpret our resultsin terms of the dimensionless memory parameter , M = 1 /ν , where the influence of priordroplet dynamics increases with M . For future reference, we provide a list of the salientdimensionless parameters related to the lattice model (2) in Table 1. Linear theory —We consider a static, equispaced lattice with droplet positions x n = nδ and h = ν − (cid:80) Nm =1 H ( x − mδ ), where δ = 2 πR/λN is the droplet arc-length separation alongthe circle of dimensionless radius R/λ . The critical value of M above which the wave forcepromotes sustained self-propulsion of the droplets, and hence oscillations of the lattice, isdetermined from the linear stability of the lattice system (2). We summarize the key results;full details are presented in [18].Coinciding with experiments [17], we consider small perturbations of the droplet positionsof the form x n ( t ) = nδ + η [ A exp(i knα + λ k t ) + c . c . ] , η (cid:28) , (3)with a concomitant perturbation to the wave field, h . Here α = 2 π/N is the angular spacingof the droplets, A is a complex amplitude, i is the imaginary unit, and c . c . denotes complex6arameter DefinitionLattice model (2) x n droplet positions N number of droplets h stroboscopic lattice wave field H single-droplet quasistatic wave kernel with characteristicwavelength, λ , decay length, l d , and amplitude, A α = 2 π/N ; δ = αR/λ = αr droplet angular separation; droplet arc-length separation ν ; M = 1 /ν dissipation rate of wave field; memory parameter ε = √ ν c − ν ; ν c = 1 /M c perturbation parameter; instability threshold for super-critical Hopf bifurcation k c ; ω c critical wavenumber and angular frequency at supercriti-cal Hopf bifurcationAmplitude equations (9) µ = α/ε ; µ c bifurcation parameter for amplitude equations; thresholdfor onset of spatial modulations A n , B n n th droplet oscillation amplitude and rotational drift c g ; σ i ; γ i ; D i group speed; growth coefficients; coupling coefficients; dif-fusion coefficientsTable 1: List of salient parameters in the lattice model (2), amplitude equations (9), andstability analysis thereof. 7onjugation of the preceding term. The eigenvalues, λ k , for each integer wavenumber, k ,satisfy the dispersion relation D k ( λ k ; ν ) = 0 [18], where D k ( λ k ; ν ) = λ k + λ k + c ν − c k λ k + ν (4)and the real constants c k are defined as c k = N (cid:88) n =1 cos( knα ) H (cid:48)(cid:48) ( nδ ) . By symmetry considerations, we need only consider discrete wavenumbers in the interval k ∈ [0 , N ], where N = (cid:98) N/ (cid:99) . Notably, the coefficients c k depend on both the number ofdroplets, N , and the droplet separation, δ .After rearranging D k ( λ k ; ν ) = 0 and writing ν = 1 /M , the eigenvalues, λ k , describingthe asymptotic linear stability of the equispaced lattice are roots of the cubic polynomial q k ( λ k ; M ), where q k ( λ k ; M ) = M λ k + ( M + 1) λ k + ( c M + 1) λ k + M ( c − c k ) . (5)As the memory parameter, M , is increased, the lattice becomes asymptotically unstable if,for some wavenumber, k , there is at least one eigenvalue, λ k , for which Re( λ k ) >
0. Wenote, however, that a fundamental property of the lattice system is its rotational invariance,characterized by the neutrally stable eigenvalue λ = 0. As we shall see in § δ , through thecoefficients c k . Specifically, the lattice can destabilize via two distinct mechanisms, dependingprimarily on the droplet separation, δ : (i) an oscillatory instability, wherein the real part ofa pair of complex-conjugate eigenvalues transitions from negative to positive as M increasesthrough M = M c = 1 /ν c , or (ii) a so-called “geometrical instability,” independent of thememory parameter, M , brought on by geometrical frustration of the lattice wave field. Wefocus our attention in this paper on case (i), which arises when c k ≤ c for all k [18]. InFigure 2(a), we present an oscillatory instability arising for N = 20 droplets: as the control8arameter, M , is increased, a single wavenumber k c = N/ λ k c ) = 0 forsome critical value of M = M c , while Im( λ k c ) = ω c (cid:54) = 0, corresponding to an oscillatoryinstability of the lattice. In terms of the coefficients c k and critical memory, M c , the angularfrequency, ω c , satisfies ω c = c M c + 1 /M c [18]. Further, that the system destabilizes to apair of complex-conjugate eigenvalues points to a Hopf bifurcation [31].The existence of a Hopf bifurcation was confirmed in [18] by an accompanying weaklynonlinear stability analysis of the equispaced lattice. The analysis presented in [18] demon-strates that, just beyond the instability threshold (0 < M − M c (cid:28) x n = nδ + [ D ( T ) + O ( ε )] + ε (cid:2) A ( T )e i( k c nα + ω c t ) + c . c . (cid:3) + O ( ε ) , (6)where 0 < ε = √ ν c − ν (cid:28)
1. The slowly varying complex amplitude A ( T ), a function of theslow timescale T = ε t , is governed by a Stuart-Landau equationd A d T = σ A − ¯ σ | A | A (7a)with an accompanying equation governing the evolution of the rotational drift, D , of thelattice, namely d D d T = ¯ γ | A | . (7b)The origin of this rotational drift may be traced back to the k = 0 mode of the dispersionrelation (5), corresponding to rotational invariance. As alluded to earlier, we will find ananalogous equation in the amplitude equations presented in § σ , ¯ σ , and ¯ γ in equations (7) are defined in terms of theparameters of the lattice system (2) [18]. We shall see how they are related to the coefficientsof the amplitude equations derived in this paper in § σ satisfies Re( σ ) >
0, corresponding to exponential growth of the oscillation amplitude. Whether the Hopfbifurcation is super- or sub-critical depends on the sign of Re(¯ σ ): For Re(¯ σ ) >
0, thebifurcation is supercritical, and subcritical when Re(¯ σ ) <
0. Further, when k c = N/
2, it isfound that ¯ γ = 0, in which case D = constant, corresponding to an arbitrary shift of thedroplet equilibrium positions. 9he stability of the equispaced lattice and its complex dependence on the system param-eters, is summarized in Figure 2(c). Each region of parameter space is color-coded accordingto the type of instability that arises for a lattice of N = 20 droplets as δ and the dimension-less spatial decay length of the wave kernel, l = l d /λ , are varied independently. Of particularinterest here are the green regions, indicating a supercritical Hopf bifurcation. A crucial fea-ture of the stability analysis just described is the value of the critical wavenumber, k c [18].When the Hopf bifurcation is supercritical, it is found that k c = (cid:98) N/ (cid:99) (except near theboundaries with subcritical Hopf bifurcations), a fact that we will exploit in § N is odd, a symmetry-breaking, oscillatory-rotary motion of the lattice arises, whichwill manifest in § ν < ν c or, equivalently, M > M c ) reveal the onset of a secondbifurcation where spatially uniform, small-amplitude oscillations of the droplet positionsgive way to spatio-temporal modulations in the droplet oscillation amplitude, arising via along-wavelength instability. To capture this second bifurcation, and the ensuing dynamics,the weakly nonlinear analysis leading to the amplitude equations (7) must be generalized toaccount for spatial, as well as temporal, effects. Starting from the lattice system (2), we use the asymptotic method of multiple scales [32]to derive a generalized set of amplitude equations, accounting for both spatial and temporalmodulations of the droplet oscillation amplitude and rotational drift. Specifically, we showthat each droplet position evolves according to x n = nδ + ε (cid:2) A n ( T )e i( k c nα + ω c t ) + c . c . + B n ( T ) (cid:3) + O ( ε ) , (8)where A n is the slowly varying complex amplitude of the n th droplet in the lattice, B n isthe rotational drift, and T = ε t is the slow time scale. As is typical in the method ofmultiple scales, eliminating secular terms at successive orders of ε yields a coupled systemof equations for A n and B n , resulting in a dCGLE coupled to a discrete mean equation:d A n d T + µ c g ∇ A n = σ A n − σ | A n | A n + µγ A n ∇ B n + µ D ∆ A n , (9a)10 .5 11.5 22.5123 Sup. HopfSub. HopfGeometric -15 -10 -5 0 5 10 15-0.0500.050.10 5 10 15 20-0.20
Figure 2: (a) Behavior of the eigenvalues of (4) in the case of a supercritical Hopf bifurcationfor N = 20 droplets. At a critical threshold M = M c , a single wavenumber k c = 10first touches Re( λ k c ) = 0. (b) The prototypical wave-field kernel, H ( x ), used in numericalsimulations presented herein (see § N = 20 droplets as the parameters l and δ are variedindependently [17, 18], as derived from (5) and (7). We delineate regimes of super- andsub-critical Hopf bifurcations, as well as geometric frustration of the equispaced lattice. Thediamond indicates the parameter values used in (a) and (b), specifically δ = l = 1 . B n d T = µ D ∆ B n + 2 µ Re [ γ A n ∇ A ∗ n ] + µγ | A n | . (9b)The operators ∇ and ∆, defined in §
2, are discrete gradient and Laplacian operators, re-spectively. In terms of the dissipation rate, ν , and the number of droplets, N , the controlparameter, µ , in the amplitude equations (9) that arises from our analysis is defined µ = αε = 2 πN √ ν c − ν . Our theory is valid for µ = O (1) or, equivalently, α ∼ ε . (We note that this conditionrequires that N is sufficiently large.) Spatially uniform oscillations for which A n ( T ) = A ( T )11nd B n ( T ) = B ( T ) are also captured by our theory, in which case the amplitude equations(9) reduce in form to the Stuart-Landau equations (7) with D = εB . Table 1 provides areference list of relevant parameters appearing in the amplitude equations (9).We pause to emphasize a few aspects of the system (9). With regards to the coefficients,the group speed, c g , growth coefficients, σ i ( i = 1 , γ i ( i = 1 , , D i ( i = 1 , D > γ are both real, while the remainingparameters are all complex with Re( σ ) > D ) >
0. Notably, σ = ¯ σ + O ( α ),which is consistent with (7a) as α → B n , appears at O ( ε ) in the expansion (8), whereas the driftis an O (1) term in (6). As discussed in Appendix A, this change arises since ¯ γ = αγ when k c (cid:46) (cid:98) N/ (cid:99) and α ∼ ε (cid:28)
1, resulting in γ = O (1); hence, the drift, D , appearing in (6) ispromoted to higher order, specifically O ( ε ).We note that, in our system, discreteness originates at the level of the underlying micro-scopic equations (2), and thus is a connate characteristic of the resultant amplitude equations(9). This feature is in contrast with discrete versions of the CGLE motivated by a directdiscretization of the continuous CGLE using standard finite difference operators [14, 15, 16].Similarly, periodicity arises from the periodicity of the lattice, rather than being imposed ex post facto through periodic boundary conditions [33, 34]. Interestingly, the system (9) isthe discrete and periodic analogue of the amplitude equations describing a host of disparatephysical phenomena presented elsewhere [35, 36, 37, 38]. This paper is organized as follows. For the interested reader, we provide further details of themultiple scales procedure involved in deriving the amplitude equations (9) from the latticemodel (2) in §
2. We then proceed to analyze the successive bifurcations of the system (9) in §
3, rationalizing the onset of the second bifurcation leading to a long-wavelength instabilityof time-periodic oscillations of the lattice, and eventually to spatiotemporal chaos. Numericalsolutions of the amplitude equations (9) beyond the second bifurcation are presented in § §
5, along with a discussion of future directions.
The weakly nonlinear analysis leading to the derivation of the amplitude equations (9) as-sumes slowly varying spatial and temporal modulations of the oscillation amplitude of eachdroplet. For a supercritical Hopf bifurcation, stable, small-amplitude oscillations arise when ν is only slightly below the critical threshold, ν c (corresponding to M slightly above M c );specifically we consider 0 < ε = √ ν c − ν (cid:28) x n = nδ and h = ν − c (cid:80) Nm =1 H ( x − mδ ). We then pose the followingasymptotic multiple-scales expansions x n ∼ nδ + ∞ (cid:88) i =1 ε i x ( i ) n ( t, T ) , h ∼ ν c N (cid:88) m =1 H ( x − mδ ) + ∞ (cid:88) i =1 ε i h ( i ) ( x, t, T ) , (10)where the slow time-scale is T = ε t .Following the typical recipe for the method of multiple scales, inserting (10) into thelattice system (2) leads to a hierarchy of problems for x ( i ) n and h ( i ) at successive orders of ε . A series of solvability conditions must then be satisfied at each order of ε to guaranteethe suppression of secular terms that would otherwise lead to unbounded solutions and aviolation of the multiple-scales ansatz . Satisfying the solvability condition arising at O ( ε )leads to the amplitude equations (9). Before arriving there, however, there are several aspectsof the current problem that complicate the weakly nonlinear, multiple-scales analysis of thelattice system (2). Our derivation falls into three stages: (i) we first solve for the wave fieldperturbations, h ( i ) , allowing us to project the wave force onto the droplet trajectories, x ( i ) n ;(ii) we next exploit the spatial decay of the wave kernel, H , and the assumed slowly varyingspatial effects to approximate the nonlocal inter-droplet dynamics by p -nearest-neighborinteractions, where p is determined by the decay length of the wave and the spacing of thelattice; (iii) finally, we consider the limit of weak asymmetry when the number of dropletsis large (equivalently, when k c departs slightly from N/ key deas below, with a full account provided in Appendix A. Our first point of interest is at O ( ε ), where we have the system ∂ x (1) n ∂t + ∂x (1) n ∂t = − x (1) n ν c N (cid:88) m =1 H (cid:48)(cid:48) (( n − m ) δ ) − ∂h (1) ∂x (cid:12)(cid:12)(cid:12)(cid:12) x = nδ , (11a) ∂h (1) ∂t + ν c h (1) = − N (cid:88) m =1 x (1) m H (cid:48) ( x − mδ ) . (11b)At this stage in conventional multiple-scales analyses of nonlinear oscillators [31, 32], oneis typically interested in solving for the perturbed oscillator position, x (1) n , alone. In thepresent case, however, we must also contend with equation (11b) governing the free surfaceperturbation, h (1) . In order to project the dynamics entirely onto the droplet trajectories, x (1) n , our first key idea is to use the form of (11b) to define the auxiliary variables, X n ,satisfying [18] ∂X n ∂t + ν c X n = x (1) n . (12)A particular solution of (11b) is then h (1) = − N (cid:88) m =1 X m H (cid:48) ( x − mδ ) . (13)We neglect the homogeneous solution of (11b), which decays exponentially on the fast time-scale, t . Now that h (1) is expressed in terms of the auxiliary variables, X n , through equation(13), the linear system (11) may be recast as a dynamical system for x (1) n and X n . Specifically,substituting (13) into (11a) yields L n x (1) = 0, where x (1) = (cid:16) x (1)1 , . . . , x (1) N (cid:17) and the linearoperator, L n , is defined as L n x (1) = ∂ x (1) n ∂t + ∂x (1) n ∂t + N (cid:88) m =1 (cid:32) x (1) n ν c − X m (cid:33) H (cid:48)(cid:48) ( δ ( n − m )) . (14)Informed by the ansatz (3), we now seek a solution to L n x (1) = 0 of the form x (1) n = A n ( T )e i( k c nα + ω c t ) + c . c . + B n ( T ) , X n = A n ( T ) ν c + i ω c e i( k c nα + ω c t ) + c . c . + 1 ν c B n ( T ) , (15)14here we recall that the critical wavenumber of instability, k c , and angular frequency, ω c ,are determined from linear theory ( § n inboth the complex amplitude, A n , and mean, B n , generalizes the spatially uniform expansion(6), which simply leads to a Stuart-Landau equation [18]. By inserting (15) into (14), wefind that L n x (1) = (cid:40) e i φ n ν c + i ω c N (cid:88) m =1 ( A n − A n − m )e − i k c mα H (cid:48)(cid:48) ( mδ ) (cid:41) + c . c . + 1 ν c N (cid:88) m =1 ( B n − B n − m ) H (cid:48)(cid:48) ( mδ ) , (16)where φ n = k c nα + ω c t . We note that, en route to obtaining (16), we first write A n − m =( A n − m − A n ) + A n and then simplify the resultant expression using the properties of thedispersion relation, namely D k ( λ k ; ν c ) = D (0; ν c ) = 0. We now arrive at the second key idea , which lies at the heart of our analysis: approximatingthe discrete convolutions arising in equation (16). When A n and B n are spatially uniform(i.e. independent of n ), the right-hand side of (16) is identically zero. To allow for spatialvariations, we approximately satisfy (16) at this order by approximating the discrete con-volutions. Recalling the angular spacing parameter, α = 2 π/N , we assume a distinguishedlimit between the relative sizes of α and ε . Specifically, we assume α ∼ ε and thus set α = µε , where µ = O (1) is the control parameter arising in the amplitude equations (9).This limiting case allows us to approximate convolutions by spatially local terms.The following approximation technique applies to any 2 π -periodic function F ( θ ) thatis slowly varying (that is, the derivatives of F are all of size O (1)) and to any periodicfunction G ( x ) = G ( x + L ) exhibiting exponential spatial decay. By denoting F m = F ( mα ),we approximate convolutions of the form N (cid:88) m =1 F n − m G ( mδ ) , (17)by first introducing an interpolating polynomial of degree 2 p passing through the points15 n − p , . . . , F n + p . As justified in Appendix B, for | m | ≤ p = O (1), we may then express F n − m = F n − αm ∇ F n + α m F n + O ( α ) , (18)where ∇ and ∆ are the central finite difference operators approximating the first and secondderivatives of F ( θ ) using 2 p + 1 points spaced α apart. The difference stencils for p = 1 and p = 2 are listed in Appendix B.We then exploit the assumed exponential decay of G to extend the polynomial approxi-mation (18) outside of the interval m = − p, . . . , p , only incurring exponentially small errorswhen substituting into the convolution, provided that p is sufficiently large. It follows thatthe convolution (17) may be approximated by N (cid:88) m =1 F n − m G ( mδ ) = F n N (cid:88) m =1 G ( mδ ) − α ∇ F n N (cid:88) m =1 a m G ( mδ ) + α ∆ F n N (cid:88) m =1 b m G ( mδ ) + O ( α ) , (19)where the periodic coefficients a m = a m + N and b m = b m + N are odd and even, respectively,defined as a m = m and b m = 12 m for | m | < N/ , with a N/ = 0 and b N/ = N / N even.The coupling integer p > | G ( ± pδ ) | (cid:28)
1; due to the exponential decayof G over the length scale l , this condition is equivalent to exp( − pδ/l ) (cid:28)
1. Moreover, if p becomes too large then the neglected terms in (18) may also be large since αp may not besufficiently small: hence, p must also satisfy αp (cid:28)
1. For the case where N is large ( α (cid:28) l ∼ δ , p = 1 and 2 adequately satisfy both of these conditions. We note that the largerthe value of p , the greater the level of coupling between droplets in the resultant amplitudeequations (9), where p = 1 corresponds to nearest-neighbor interactions.Employing the slowly varying approximation (19) in (16) yields the sought-after approx-16mation to the discrete convolutions. Specifically, we obtain L n x (1) + (cid:40) − α ∇ A n e i φ n ν c + i ω c N (cid:88) m =1 a m e − i k c mα H (cid:48)(cid:48) ( mδ ) + α ∆ A n e i φ n ν c + i ω c N (cid:88) m =1 b m e − i k c mα H (cid:48)(cid:48) ( mδ ) (cid:41) + c . c . + α ∆ B n ν c N (cid:88) m =1 b m H (cid:48)(cid:48) ( mδ ) = O ( α ) . (20)We note that there is not a ∇ B n term in (20) as the symmetry of the wave kernel, H ,determines that its coefficient vanishes, specifically (cid:80) Nm =1 a m H (cid:48)(cid:48) ( mδ ) = 0. Recalling ourassumption that α ∼ ε , terms of O ( α n ) in equations (20) are consequently promoted to O ( ε n +1 ), appearing as secular terms (either those constant in t or proportional to e i φ n ) onthe right-hand of the expansion of equation (2a). Thus, the ∆ A n and ∆ B n terms in equation(20) are destined to become the diffusion-like terms in the amplitude equations (9). Finally, our third key idea is to use the observed property of the supercritical Hopf bifur-cations, outlined in Figure 2(a), that k c (cid:46) (cid:98) N/ (cid:99) . We first define χ = N/ − k c , where wetypically find that χ = 0 or χ = 1 /
2. Then, by recalling that α = 2 π/N , we use the form k c = N/ − χ to recast the ∇ A n coefficient in (20) as αν c + i ω c N (cid:88) m =1 a m e − i k c mα H (cid:48)(cid:48) ( mδ ) = i αν c + i ω c N (cid:88) m =1 a m ( − m sin( mαχ ) H (cid:48)(cid:48) ( mδ ) , which is an even function of α for χ (cid:54) = 0, and vanishes otherwise. In the former case, thisterm is expected to be of size O ( α ) as α →
0, which may be verified numerically. Further,we are prompted to define the O (1) group velocityˆ c g = 1 α N (cid:88) m =1 a m e − i k c mα H (cid:48)(cid:48) ( mδ ) ν c + i ω c , which vanishes when k c = N/
2. We then acknowledge—as a corollary of our second keyidea—that the term α ˆ c g ∇ A n appears as a secular term at O ( ε ), ultimately resulting in theadvective term in our dCGLE (9a). (We note that ˆ c g is divided by a further coefficient toyield the c g in (9); see Appendix A for clarification.) In this final step, equation (20) reducesto L n x (1) = O ( ε ), verifying our ansatz for x (1) n .17s shown in Appendix A, a combination of the foregoing three key ideas is used to sys-tematically derive the amplitude equations (9) from the lattice system (2). The remainingterms involving time derivatives and the nonlinear coupling terms in (9) are obtained byeliminating higher-order secular terms, both those that are promoted from O ( ε ) and othersthat appear at O ( ε ). We now proceed to analyze the stability of the amplitude equations(9), elucidating the second bifurcation leading to spatiotemporal modulation of the dropletoscillation amplitude and drift as the control parameter, µ , is decreased from infinity (cor-responding to ν < ν c ). As discussed in §§ ν < ν c , or equiv-alently, for µ < ∞ . In this section, we elucidate the mechanism leading to a modulationalinstability of the spatially uniform solution of the amplitude equations (9). The onset ofspatial amplitude modulations in the canonical complex Ginzburg-Landau equation is theeponymous Benjamin-Feir-Newell (BFN) instability, after its discovery in describing the in-stability of periodic surface gravity (Stokes) waves [39, 40]. As we shall see, in our system,this instability takes a slightly different form due to the coupling of the complex amplitude, A n , with the mean, B n . In what follows, we conduct a linear stability of the amplitude equa-tions (9), the computations for which are standard [2], but lengthy. We therefore highlightonly the key features here. In the case of a supercritical Hopf bifurcation, where Re(¯ σ ) ∼ Re( σ ) > α (cid:28) A (0) n ( T ) = ρ exp(iΩ T ) and B (0) n ( T ) = µγ ρ T + constant , (21)where the modulus and angular frequency of the complex amplitude are ρ = (cid:112) Re( σ ) / Re( σ ) and Ω = Im( σ σ ∗ ) / Re( σ ) .
18e now consider small perturbations about the spatially uniform state of the form A n ( T ) = A (0) n ( T )(1 + A n ( T )) , B n ( T ) = B (0) n ( T ) + B n ( T ) , where 0 < | A n | ∼ | B n | (cid:28)
1. Substituting this ansatz into (9) and neglecting terms ofquadratic order and higher, we obtain a linear system governing the perturbations A n and B n , supplemented by an additional equation for C n = A ∗ n . The resulting linear systemmay then be diagonalized by considering a discrete Fourier transform in n . Specifically, weconsider solutions of the form A n = N − (cid:88) ξ =0 ˆ A ξ e i ξnα , B n = N − (cid:88) ξ =0 ˆ B ξ e i ξnα , C n = N − (cid:88) ξ =0 ˆ C ξ e i ξnα , where the wavenumber, ξ , is an integer. Under this transformation, we obtaind ˆ A ξ d T = M ξ ( µ ) ˆ A ξ , (22)where ˆ A ξ = ( ˆ A ξ , ˆ B ξ , ˆ C ξ ) T andM ξ ( µ ) = − µ c g ∇ ξ − ρ σ + µ D ∆ ξ µγ ∇ ξ − ρ σ µρ ( γ ∗ ∇ ξ + γ ) µ D ∆ ξ µρ ( γ ∇ ξ + γ ) − ρ σ ∗ µγ ∗ ∇ ξ − µ c ∗ g ∇ ξ − ρ σ ∗ + µ D ∗ ∆ ξ . (23)For the case of nearest-neighbor interactions ( p = 1), the Fourier symbols, ∇ ξ and ∆ ξ , ofthe difference operators, ∇ and ∆, are defined as ∇ ξ = i sin( ξα ) α and ∆ ξ = 2(cos( ξα ) − α . The eigenvalues, Λ ( j ) ξ (for j = 1 , , ξ ( µ ) determine the stability of the spatiallyuniform state (21), where the dependence of Λ ( j ) ξ on µ is presented in Figure 3. The onset ofinstability is determined by the eigenvalue (or one of a pair of complex-conjugate eigenvalues)of M ξ ( µ ) with maximal real part; we denote this eigenvalue as Λ ξ ( µ ) for each wavenumber, ξ ∈ { , , . . . , N − } . At a given value of µ >
0, the perturbed system (22) is neutrallystable if Re(Λ ξ ( µ )) ≤ ξ , and unstable otherwise. We first note, however, that therotational and temporal invariance of the spatially uniform state (21) implies that Λ = 0has multiplicity two. The remaining eigenvalue for ξ = 0 is − ρ Re( σ ) = − σ ) < N = 36 (left) and N = 37 (right)computed for nearest-neighbor interactions ( p = 1). ( a ) The case ξ = 1. The real part ofthe eigenvalues Λ ( j )1 for j = 1 , , µ , where the asymptotic behavior for µ (cid:29) µ decreases determines the instability threshold for ξ = 1 (red circle). For N = 37, spuriousinstability is predicted by the amplitude equations (9) for µ (cid:29)
1, a regime inconsistent withthe µ = O (1) assumption. ( b ) The instability threshold, µ ξ , as a function of the wavenumber, ξ , for ξ ≤ N/
2. The long-wave mode, ξ = 1 (red circle), is the first to destabilize as µ isdecreased.corresponding to a stable perturbation from the periodic state for all µ >
0. Hence, if aninstability to the periodic state arises, then it is for ξ (cid:54) = 0, corresponding to the emergenceof a nontrivial spatial pattern. Moreover, for µ (cid:29) ξ (cid:54) = 0, the diagonal elements of thematrix M ξ ( µ ) dominate, from which we infer that the eigenvalues of M ξ ( µ ) are approximatedby Λ (1) ξ ∼ µ D ∆ ξ , Λ (2) ξ ∼ µ ( D ∆ ξ − c g ∇ ξ ) , Λ (3) ξ ∼ µ ( D ∗ ∆ ξ − c ∗ g ∇ ξ ) , (24)corresponding to advection and diffusion of perturbations. We recall that Re( D ) > D >
0, and ∆ ξ < ξ (cid:54) = 0. Hence, for c g = 0 (arising when k c = N/ µ (cid:29) N = 36). However,20or c g (cid:54) = 0, there are regimes in which the spatially uniform state is spuriously predicted tobe unstable for µ (cid:29) µ = O (1) assumption), as presentedin Figure 3 for N = 37. To sidestep this spurious prediction, we therefore define µ ξ > µ at which Re(Λ ξ ( µ ξ )) = 0 and Re(Λ ξ ( µ )) is a decreasing function of µ at µ = µ ξ (see Figure 3( a ) for reference). In the µ = O (1) regime of interest, the system istherefore unstable to perturbations for µ < µ c = max ξ µ ξ . A crucial feature of the system (9) is the coupling of the complex amplitude, A n , with themean, B n , which acts to drive the instability of the spatially uniform state (21). We observethat variations in B n act as a source term in the amplitude, A n , thus promoting spatialvariations in A n . For large µ , the diffusion term in (9a) counteracts the growth of spatialvariations of A n , but this smoothing effect may become subdominant to the source term when µ is sufficiently small. Likewise, a similar competition between diffusion and the excitationof spatial variations in B n driven by variations in A n is apparent in equation (9b). As aconsequence, this coupling provides a positive feedback loop for the emergence of spatialvariations, whereas a BFN-like mechanism instead relies on sufficiently small dissipation. Infact, a true BFN instability only arises in the system (9) when the coupling coefficients, γ i ,and the group speed, c g , vanish. A necessary condition for instability in this artificial caseis Re( σ D ∗ ) < ξ ( µ ), we observe that, similar to the BFNinstability, the amplitude equations (9) exhibit a long-wave instability at ξ = 1 (see Figure3); hence, in the cases considered here, we have µ c = µ . To capture this phenomenonanalytically, we consider the special where N is even and k c = N/
2, resulting in c g = 0 and γ = 0. In this case, the spatially uniform state destabilizes via a real eigenvalue, for whichthe corresponding value of µ ξ satisfies det M ξ ( µ ξ ) = 0. This condition is satisfied by µ ξ = 0or when µ ξ has a positive solution to (cid:104) | D | D ∆ ξ (cid:105) µ ξ − ρ ∆ ξ (cid:104) D ∆ ξ Re( σ D ∗ ) + ∇ ξ Re( γ ∗ γ D ) (cid:105) µ ξ − ρ ∇ ξ Im( σ γ ∗ ) Im( γ ) = 0 . As ∇ ξ = −∇ N − ξ and ∆ ξ and ∆ N − ξ , the coefficients of this equation are invariant under the21apping ξ → N − ξ , so µ ξ = µ N − ξ for all ξ . Moreover, to elucidate the dependence of µ ξ on ξ , we consider the limiting case ξα (cid:28)
1. In this limit, we observe that ∇ ξ ∼ ∆ ξ ∼ − ξ ,yielding µ ξ = κρ √ /ξ , where κ is a real and positive root of the quartic polynomial P ( κ ) = (cid:104) | D | D (cid:105) κ + (cid:104) D Re( σ D ∗ ) + Re( γ ∗ γ D ) (cid:105) κ − Im( σ γ ∗ ) Im( γ ) . If P has a real and positive root then the form µ ξ ∼ ξ − suggests that the first wavenumberto destabilize is ξ = 1 (or ξ = N −
1, by symmetry); this long-wave instability is thusconsistent with our slowly varying approximation of convolutions (see § We proceed to explore the dependence of µ c on the inter-droplet spacing, δ , and spatialdecay length, l , of the wave kernel, H . Motivated by the form of the wave field arisingin the bouncing-droplet system [24], a candidate wave kernel that satisfies the assumedperiodicity, exponential decay and quasi-monochromaticity may be derived by projectingthe form of the dimensionless radially symmetric wave F ( r ) = A J (2 πr )sech( r/l ) onto acircle of circumference 2 πr = N δ , where r = R/λ [18]. Here J is the Bessel function ofthe first kind with order zero and A is the amplitude of the wave. The resultant algebraicform of the wave kernel is H ( x ) = F (cid:18) r sin x r (cid:19) , (25)where an example of this wave kernel is given in Figure 2( b ). For the numerical resultspresented herein, we consider A = 0 . A ; increasing A simply serves to increase the amplitudeof the wave field accompanying the equispaced lattice, the main consequence of which is aconcomitant decrease in the critical memory, M c .As presented in Figure 4, the dependence of the onset of spatial modulations on thesystem parameters, l and δ , can be quite intricate. Near the boundaries between super- andsub-critical Hopf bifurcations (where geometric and subcritical Hopf instabilities arise in thewhite regions in Figure 4), µ c can be very large—a feature that appears to be correlated22ith k c departing from N/ µ = O (1) assumption under which theamplitude equations (9) were derived. Away from these boundaries, we observe regions inwhich µ c = O (1); indeed, simulation of the lattice system (2) reveals favorable agreement ofthe numerical and theoretical instability threshold. (This agreement improves as α = 2 π/N becomes smaller, consistent with our assumptions; see supplementary material.) Near themiddle of each ‘band’ in which supercritical Hopf bifurcations arise, we observe that µ c = 0,corresponding to the prediction of unconditional stability of the spatially uniform solution(21). Finally, we remark that µ c = O (1) is most apparent for l ∼ δ , a regime in which theinter-droplet coupling is dominated by nearest-neighbor interactions, as might be anticipatedfrom our local approximation of convolutions (see § l is much smaller than δ , eachdroplet only interacts weakly with all other droplets (including its nearest neighbors), andspatial variations are less propitious. We proceed to explore the nonlinear dynamics predicted by the amplitude equations (9)beyond the onset of spatial variations, µ < µ c . The amplitude equations are evolved usinga spectral method over the droplet number, n , and a fourth-order Runge-Kutta method intime, for which the linear terms are transformed using an integrating factor (see Appendix Cfor details) [41]. Initially considering µ just below the instability threshold, 0 < µ c − µ (cid:28) A n = ρ + ζ sin( nα ) and B n = ζ sin( nα ) ( ζ (cid:28)
1) until a periodic state is attained. Thereafter, we decrement µ by 0.02 and initialize the following simulation at the final values obtained in the precedingsimulation. The MATLAB code used to simulate these dynamical states is provided in thesupplementary material.As suggested by Figure 4, there is a vast parameter space we could explore with equations(9) by varying the parameters l , δ , and N . To fix ideas, we focus on the case of two adjacentdroplet numbers, N = 40 and N = 41, and set l = δ = 2 .
6, which serves to elucidatethe key phenomenology exhibited by the amplitude equations (9). A deeper exploration ofthe parameter space is reserved for future work. Before we present the solutions, we note23igure 4: The onset of spatial modulations for N = 40 droplets, as predicted by the linearstability analysis of the amplitude equations (9) for the wave kernel defined in equation (25)and p = 1 nearest neighbors. ( a ) When the initial instability of the equidistant lattice arises via a supercritical Hopf bifurcation, we color-code each value of the spacing parameter, δ ,and the decay length, l , by min( µ c , µ < µ c . Our theoryis valid when µ c = O (1). When µ c = 0, the periodic state is predicted to be unconditionallystable. Large values of µ c (i.e. those exceeding the threshold of 2) arise near the boundarybetween super- and sub-critical Hopf bifurcations. ( b ) The corresponding value of k c , aspredicted by the linear stability analysis, demonstrating the correlation between large µ c and k c < N/ µ , since the coefficientsare fixed for a particular choice of wave kernel (25) and constituent parameters l , δ , and N . Thus, varying µ corresponds to traversing a particular path through parameter space , incontrast to varying each coefficient in the amplitude equations independently. As we shallsee, this variation in µ gives rise to a series of bifurcations between qualitatively differentdynamical behaviors. 24n the case of 40 droplets (see Figure 5), close to the onset of spatial modulations at µ = µ c , we first observe a time-independent solution for | A n | (panel (i)), which destabilizesinto a periodic, breather-like state (panel (ii)) with dips apparent in | A n | . These dips persistas µ c − µ is increased (panels (iii) and (iv)), but | A n | is instead constant in time. When µ c − µ is sufficiently large (panels (v) and beyond), a robust stationary bright soliton emergesover a large range of µ , before destabilizing far from the instability threshold (panel (vii)),leading to chaotic fluctuations in the soliton form. We note that when N = 40, the groupspeed, c g , is identically zero. For N = 41 droplets and non-zero group speed (Figure 6),similar dynamical transitions occur, although the asymmetry of the system instead yieldstraveling waves and propagating solitons (panels (i) and (vii)). Notably, we also observeparameter regimes for which dark breathers (panel (ii)) and dark solitons (panels (iii)–(v))arise, characterized by the sharp dips in the amplitude, | A n ( T ) | , towards zero, which are notpresent when N = 40. All of the aforementioned features arise over a large range of N , andthus appear to be canonical features of the discrete amplitude equations (9). We also notethat the jumps in the bounds of | A n | and B n are indicative of hysteresis between dynamicalstates, an effect to be explore in greater detail elsewhere. In this paper, we have presented a rigorous mathematical framework to derive a discrete setof amplitude equations from a driven and dissipative oscillator model, inspired by the physicsof droplet lattices bouncing on a vibrating fluid bath. Our systematic derivation providesa direct link between the constitutive properties of the lattice model (2) (specifically, thewave kernel, H ) and the coefficients arising in the amplitude equations (9). A linear stabilityof the amplitude equations (9) reveals the importance of the coupling to the discrete meanequation (9b) in destabilising the system from a spatially uniform state, leading to spatialmodulations in the droplet amplitude following a second bifurcation, similar in spirit tothe BFN instability. Beyond this second bifurcation, numerical solutions of the amplitudeequations (9) reveal a fascinating family of dynamical behaviours including bright and darksolitons, breather states, and traveling waves. As computational models of the droplet25igure 5: The nonlinear dynamics predicted by the amplitude equations (9) for N = 40droplets and µ < µ c = 0 . a ) The upper and lower bounds of | A n ( T ) | attained overthe entire simulation. We identify seven dynamical regimes, which are denoted by Romannumerals and divided by the dotted vertical lines. ( b ) The upper and lower bounds of δB n ( T ) = B n ( T ) − (cid:104) B n ( T ) (cid:105) attained over the entire simulation, where angled bracketsdenote the average over n . ( c ) Snapshots of | A n | ( T fixed) for each dynamical regime, where µ c − µ takes values (i) 0.03, (ii) 0.08, (iii) 0.12, (iv) 0.14, (v) 0.2, (vi) 0.4, and (vii) 0.48. ( d )Corresponding spacetime plots of | A n ( T ) | for the parameter values in ( c ). The dashed linescorrespond to the snapshot time, T = 1200. The system parameters here δ = 2 . l = 2 . p = 1 nearest neighbors, yielding k c = 20 and ν c = 1 . N = 41droplets and µ < µ c = 0 . a ) The upper and lower bounds of | A n ( T ) | attained overthe entire simulation. We identify seven dynamical regimes, which are denoted by Romannumerals and divided by the dotted vertical lines. ( b ) The upper and lower bounds of δB n ( T ) = B n ( T ) − (cid:104) B n ( T ) (cid:105) attained over the entire simulation, where angled bracketsdenote the average over n . ( c ) Snapshots of | A n | ( T fixed) for each dynamical regime, where µ c − µ takes values (i) 0.014, (ii) 0.04, (iii) 0.06, (iv) 0.08, (v) 0.1, (vi) 0.15, and (vii) 0.3. ( d )Corresponding spacetime plots of | A n ( T ) | for the parameter values in ( c ). The dashed linescorrespond to the snapshot time, T = 1200. The system parameters here δ = 2 . l = 2 . p = 1 nearest neighbors, yielding k c = 20 and ν c = 1 . x n + ˙ x n = − ∂h∂x ( x n , t ) , (26a) P h = N (cid:88) m =1 H ( x − x m ) , (26b)with the linear operator P = ∂/∂t + ν in equation (2b) serving as a particular example.The novelty of the model (26) is that the inter-particle coupling potential, h , is dynamic ,continuously evolving with the particle motion, rather than being fixed in space or withrespect to the particles. Other potential choices of P are numerous, and could lead to aneven richer family of dynamics. For a particular choice of P , if the bifurcation leading toinstability of the oscillator is of supercritical Hopf type, then one should expect a complexGinzburg-Landau equation in the vicinity of the bifurcation point. However, the preciseform of the amplitude equations will change depending on the type of primary bifurcationthat arises and the inherent symmetries of the system [1, 3]. Such an investigation may leadto further insights into the dynamics and pattern-forming behaviour of active particles incomplex environments [26]. Appendix A Derivation of the Ginzburg-Landau andmean equations
Further details are provided of the multiple-scales expansion leading to the complex Ginzburg-Landau and mean equations (9). The general procedure toward obtaining equations (9) is tosubstitute the asymptotic expansions (10) into (2) and gather successive powers of ε . At eachsuccessive order, we suppress resonant terms proportional to e i φ n ( t ) (where φ n ( t ) = k c nα + ω c t )or those constant in t . To extract all relevant terms at each order, we must introduce aux-iliary variables to solve for the free surface, h , and also expand convolutions in the mannersummarized in §
2. For notational efficiency, we denote H m = H ( mδ ), H (cid:48) m = H (cid:48) ( mδ ), and soforth.At leading order we obtain ∂h (0) ∂x (cid:12)(cid:12)(cid:12)(cid:12) x = nδ = 0 , h (0) ( x ) = 1 ν c N (cid:88) m =1 H ( x − mδ ) , (27)28eflecting the fact that the free-surface gradient beneath each droplet vanishes in the steady-state. Consequently, all odd-derivatives of h (0) vanish beneath each droplet at equilibrium,a fact we will make repeated use of in simplifying the forthcoming terms in arising ourexpansion.The problem at O ( ε ) has already been discussed in § O ( ε ), remembering that terms from (20) are promoted to O ( ε ) after setting α = µε . Theequations at O ( ε ) are ∂ x (2) n ∂t + ∂x (2) n ∂t = − (cid:26) ∂h (2) ∂x + x (2) n ∂ h (0) ∂x + x (1) n ∂ h (1) ∂x (cid:27) (cid:12)(cid:12)(cid:12) x = nδ , (28) ∂h (2) ∂t + ν c h (2) = h (0) + N (cid:88) m =1 (cid:26) x (1)2 m H (cid:48)(cid:48) ( x − mδ ) − x (2) m H (cid:48) ( x − mδ ) (cid:27) . (29)Following the procedure outlined in § h (2) by introducing two further auxiliaryvariables, Y n and Z n , satisfying ∂Y n ∂t + ν c Y n = 12 x (1)2 n , ∂Z n ∂t + ν c Z n = x (2) n , (30)chosen to match the coefficients of the wave kernel on the right-hand side of (29). Hence, aparticular solution of (29) is h (2) = 1 ν c h (0) + N (cid:88) m =1 { Y m H (cid:48)(cid:48) ( x − mδ ) − Z m H (cid:48) ( x − mδ ) } . (31)By substituting the form of x (1) n , given by equation (15), into the first equation of (30), wefind Y n = 1 ν c (cid:20) | A n | + 12 B n (cid:21) + B n (cid:20) A n ν c + i ω c e i φ n + c.c. (cid:21) + 12 (cid:20) A n ν c + 2i ω c e φ n + c.c. (cid:21) . (32)Substituting (31)–(32) into (28), and then using (27), yields L n x (2) = (cid:104) B n + { A n e i φ n + c . c . } (cid:105)(cid:20) ν c N (cid:88) m =1 B n − m H (cid:48)(cid:48)(cid:48) m + (cid:26) e i φ n N (cid:88) m =1 A n − m e − i k c mα ν c + i ω c H (cid:48)(cid:48)(cid:48) m + c . c . (cid:27)(cid:21) − ν c N (cid:88) m =1 (cid:20) B n − m + | A n − m | (cid:21) H (cid:48)(cid:48)(cid:48) m − (cid:20) e i φ n N (cid:88) m =1 A n − m B n − m ν c + i ω c e − i k c mα H (cid:48)(cid:48)(cid:48) m + c . c . (cid:21) − (cid:20) e φ n N (cid:88) m =1 A n − m e − k c mα H (cid:48)(cid:48)(cid:48) m ν c + 2i ω c + c . c . (cid:21) . (33)29nalogous to § B n − m onthe right-hand side of (33), and similarly for all other terms of this form (such as A n − m , B n − m , etc.). After some arduous algebra, we reduce (33) to the highly simplified form L n x (2) = 2Re (cid:20) N (cid:88) m =1 e − i k c mα H (cid:48)(cid:48)(cid:48) m ν c + i ω c (cid:21) | A n | + α (cid:2) ˆ γ e i φ n A n ∇ B n + c . c . (cid:3) + 2 α Re (cid:2) ˆ γ A n ∇ A ∗ n (cid:3) + (cid:2) ˆ c A n e φ n + c . c . (cid:3) + α (cid:2) ˆ c e φ n A n ∇ A n + c . c . (cid:3) + O ( α ) , (34)where the complex coefficients are defined asˆ γ = N (cid:88) m =1 (cid:18) e − i k c mα ν c + i ω c − ν c (cid:19) a m H (cid:48)(cid:48)(cid:48) m , ˆ γ = N (cid:88) m =1 (cid:18) ν c − e i k c mα ν c − i ω c (cid:19) a m H (cid:48)(cid:48)(cid:48) m , ˆ c = N (cid:88) m =1 (cid:20) e − i k c mα ν c + i ω c − e − k c mα ν c + 2i ω c ) (cid:21) H (cid:48)(cid:48)(cid:48) m , ˆ c = N (cid:88) m =1 (cid:18) e − k c mα ν c + 2i ω c − e − i k c mα ν c + i ω c (cid:19) a m H (cid:48)(cid:48)(cid:48) m . (35)To apply the weak-asymmetry approximation developed in § | A n | ,we first note that2Re (cid:20) N (cid:88) m =1 e − i k c mα H (cid:48)(cid:48)(cid:48) m ν c + i ω c (cid:21) = 2Re (cid:20) i ν c + i ω c (cid:21) N (cid:88) m =1 ( − m sin( mαχ ) H (cid:48)(cid:48)(cid:48) m , where χ = N/ − k c . By a similar argument to which the group velocity was promotedto O ( ε ) in § | A n | has size O ( α ) for α (cid:28) χ = O (1), so thecorresponding term in (34) should likewise appear at O ( ε ). We thus write2Re (cid:20) N (cid:88) m =1 e − i k c mα H (cid:48)(cid:48)(cid:48) m ν c + i ω c (cid:21) | A n | = α ˆ γ | A n | , where the O (1) real coefficient ˆ γ is defined asˆ γ = 2 α Re (cid:20) N (cid:88) m =1 e − i k c mα H (cid:48)(cid:48)(cid:48) m ν c + i ω c (cid:21) . (36)In a similar spirit, we deduce that the coefficient ˆ c is size O ( α ) when χ (cid:54) = 0 and zerootherwise. Therefore the non-secular terms in (34) (those that are proportional to e φ n ) areboth of size O ( α ) and so should actually appear at O ( ε ). However, as these terms will stillbe non-secular at that order, they play no role in the derived amplitude equations for A n and B n . We conclude that all the inhomogeneities in (34) (which are of size O ( α )) should30nstead appear at O ( ε ). Hence, at O ( ε ), we have L n x (2) = 0, which is identical to theproblem for x (1) . Akin to the solution ansatz at O ( ε ), we therefore pose x (2) n = E n ( T ) + (cid:2) C n ( T )e i φ n + c.c. (cid:3) , Z n = 1 ν c E n ( T ) + (cid:20) C n ( T ) ν c + i ω c e i φ n + c.c. (cid:21) , which, when applying the approximation to the convolution (see § O ( α ) terms that arise out of the approximation tothis convolution appear at O ( ε ), which is beyond the order presented in this calculation.At O ( ε ), we have a system for x (3) n and h (3) , namely ∂ x (3) n ∂t + ∂x (3) n ∂t + x (3) n ∂ h (0) ∂x (cid:12)(cid:12)(cid:12)(cid:12) x = nδ = − (cid:34) ∂ x (1) n ∂t∂T + ∂x (1) n ∂T (cid:35) − (cid:20) ∂h (3) ∂x + x (1) n ∂ h (2) ∂x + 12 x (1)2 n ∂ h (1) ∂x + x (2) n ∂ h (1) ∂x + 16 x (1)3 n ∂ h (0) ∂x (cid:21) (cid:12)(cid:12)(cid:12) x = nδ (37)and ∂h (3) ∂t + ν c h (3) = − (cid:20) ∂h (1) ∂T − h (1) (cid:21) + N (cid:88) m =1 (cid:26) x (1) m x (2) m H (cid:48)(cid:48) ( x − mδ ) − x (3) m H (cid:48) ( x − mδ ) − x (1)3 m H (cid:48)(cid:48)(cid:48) ( x − mδ ) (cid:27) . (38)Appended to the right-hand side of (37) will be the terms promoted from both O ( ε ) and O ( ε ). We follow an identical procedure to our analysis at O ( ε ) and O ( ε ). First we in-troduce three auxiliary variables ( § h (3) . We then substitute this solution into (37), along with the dropletpositions and wave field terms computed from lower orders, and then apply the slowly vary-ing approximation ( § L n x (3) = RHS, where the right-hand-side (RHS) iscomposed of terms that are constant in t , terms with coefficients e ± i φ n ( t ) , and non-secularterms (whose form can be ignored at this stage). For a bounded solution, we require thatthe constant and e i φ n secular terms have vanishing coefficients, which yields the followingevolution equations for the complex amplitude, A n , and the real drift, B n :ˆ σ d A n d T + µ ˆ c g ∇ A n = ˆ σ A n − ˆ σ | A n | A n + µ ˆ γ A n ∇ B n + µ ˆ D ∆ A n , (39a)ˆ b d B n d T = µ ˆ D ∆ B n + 2 µ Re (cid:2) ˆ γ A n ∇ A ∗ n (cid:3) + µ ˆ γ | A n | , (39b)31here µ = α/ε and we have neglected terms whose coefficients are of size O ( α ) at O ( ε )( § C n and E n ) with-out proceeding to O ( ε ) and higher. However, a satisfactory approximation is obtained byconsidering A n and B n alone, which form a closed system.Recalling that D k ( λ ; ν ) is the dispersion relation (4), the coefficients (other than the ˆ γ i defined in equations (35) and (36)) appearing in (39a) are as follows:ˆ σ = ∂ D k c ∂λ (i ω c ; ν c ) , ˆ σ = ∂ D k c ∂ν (i ω c ; ν c ) , ˆ σ = 32 ν c N (cid:88) m =1 H (cid:48)(cid:48)(cid:48)(cid:48) m − N (cid:88) m =1 H (cid:48)(cid:48)(cid:48)(cid:48) m Re (cid:20) e − i k c mα ν c + i ω c (cid:21) + 12 N (cid:88) m =1 e − k c mα H (cid:48)(cid:48)(cid:48)(cid:48) m ν c + 2i ω c − N (cid:88) m =1 H (cid:48)(cid:48)(cid:48)(cid:48) m e − i k c mα ν c + i ω c , ˆ D = 1 ν c + i ω c N (cid:88) m =1 b m e − i k c mα H (cid:48)(cid:48) m , while those that appear in (39b) areˆ b = ∂ D ∂λ (0; ν c ) = 1 + 1 ν c N (cid:88) m =1 H (cid:48)(cid:48) m , ˆ D = 1 ν c N (cid:88) m =1 b m H (cid:48)(cid:48) m . Upon dividing (39a) by ˆ a and (39b) by ˆ b we arrive at equations (9) in the main text, where( c g , σ , σ , γ , D ) = (ˆ c g , ˆ σ , ˆ σ , ˆ γ , ˆ D ) / ˆ σ and ( D , γ , γ ) = ( ˆ D , ˆ γ , ˆ γ ) / ˆ b . Appendix B Approximation of convolutions
We introduce an interpolating polynomial of degree 2 p passing through the points F n − p , . . . , F n + p ,where F m = F ( mα ), so that F n − m = F n − m D F n + m D F n − . . . + m p (2 p )! D p F n (40)for | m | ≤ p . The operators D j /α j are the symmetric finite difference operators approximatingthe j th derivative of F ( θ ) over 2 p + 1 grid points spaced α apart. Hence, as F ( j ) ( αn ) = O (1)(by assumption) and D j F n α j = d j F d θ j ( nα ) + O ( α q )32or some integer q >
0, we conclude that D j F n /α j = O (1). Therefore, for m = − p, . . . , p , wemay recast (40) in the following form: F n − m = F n − αm ∇ F n + α m F n + O ( α ) , where ∇ and ∆ are the central finite difference operators approximating the first and secondderivatives of F ( θ ) using 2 p + 1 points spaced α apart.For example, for p = 1, the symmetric finite difference stencils are defined as ∇ F n = 12 α (cid:0) F n − + F n +1 (cid:1) and ∆ F n = 1 α (cid:0) F n − − F n + F n +1 (cid:1) , and the O ( α ) term in (18) vanishes. For p = 2, we instead define ∇ F n = 1 α (cid:18) F n − − F n − + 23 F n +1 − F n +2 (cid:19) and ∆ F n = 1 α (cid:18) − F n − + 43 F n − − F n + 43 F n +1 − F n +2 (cid:19) . Appendix C Numerical implementation
To evolve the amplitude equations (9), we apply a discrete Fourier transform to A n and B n ,and then introduce an integrating factor to integrate the linear components exactly [41].Specifically, we denote the discrete Fourier transform of A n as ˆ A k = F k [ A n ], and likewisefor B n . By applying the Fourier transform to the amplitude equations (9), we obtaind ˆ A k d T + M k ˆ A k = F k (cid:2) µγ A n ∇ B n − σ | A n | A n (cid:3) , (41a)d ˆ B k d T + N k ˆ B k = F k (cid:2) µ Re [ γ A n ∇ A ∗ n ] + µγ | A n | (cid:3) , (41b)where M k = µ c g ∇ k − σ − µ D ∆ k and N k = − µ D ∆ k . We recall that ∇ k and ∆ k are the Fourier multipliers of the difference operators ∇ and ∆,respectively. To account for the stiffness manifest in the operators M k and N k , we introduce33n integrating factor. When evolving from time T = T n and T = T n +1 , we therefore recast(41) as dd T (cid:0) ˆ A k e M k ( T − T n ) (cid:1) = F k (cid:2) µγ A n ∇ B n − σ | A n | A n (cid:3) e M k ( T − T n ) , (42a)dd T (cid:0) ˆ B k e N k ( T − T n ) (cid:1) = F k (cid:2) µ Re [ γ A n ∇ A ∗ n ] + µγ | A n | (cid:3) e N k ( T − T n ) , (42b)and then evolve the new variables, ˆ A k e M k ( T − T n ) and ˆ B k e N k ( T − T n ) , using a fourth-order Runge-Kutta method. For the numerical results presented in §
4, we use a time step of 0.01. MAT-LAB code implementing this numerical scheme is provided in the supplementary material.
References [1] I. S. Aranson and L. Kramer. The world of the complex Ginzburg–Landau equation.
Reviews of Modern Physics , 74(1):99, 2002.[2] V. Garc´ıa-Morales and K. Krischer. The complex Ginzburg–Landau equation: an in-troduction.
Contemporary Physics , 53(2):79–95, 2012.[3] M. C. Cross and P. C. Hohenberg. Pattern formation outside of equilibrium.
Reviewsof Modern Physics , 65(3):851, 1993.[4] A. C. Newell, T. Passot, and J. Lega. Order parameter equations for patterns.
AnnualReview of Fluid Mechanics , 25(1):399–453, 1993.[5] P. Coullet, L. Gil, and J. Lega. Defect-mediated turbulence.
Physical Review Letters ,62(14):1619, 1989.[6] D. Tanaka and Y. Kuramoto. Complex Ginzburg-Landau equation with nonlocal cou-pling.
Physical Review E , 68(2):026219, 2003.[7] V. Garc´ıa-Morales and K. Krischer. Nonlocal complex Ginzburg-Landau equation forelectrochemical systems.
Physical Review Letters , 100(5):054101, 2008.[8] J. Denk, L. Huber, E. Reithmann, and E. Frey. Active curved polymers form vortexpatterns on membranes.
Physical Review Letters , 116(17):178301, 2016.349] T. H. Tan, J. Liu, P. W. Miller, M. Tekant, J. Dunkel, and N. Fakhri. Topologicalturbulence in the membrane of a living cell.
Nature Physics , 16(6):657–662, 2020.[10] A. C. Newell and J. A. Whitehead. Finite bandwidth, finite amplitude convection.
Journal of Fluid Mechanics , 38(2):279–303, 1969.[11] L. A. Segel. Distant side-walls cause slow amplitude modulation of cellular convection.
Journal of Fluid Mechanics , 38(1):203–224, 1969.[12] M. C. Cross. Derivation of the amplitude equation at the Rayleigh–B´enard instability.
Physics of Fluids , 23(9):1727–1731, 1980.[13] N. Stoop, R. Lagrange, D. Terwagne, P. M. Reis, and J. Dunkel. Curvature-inducedsymmetry breaking determines elastic surface patterns.
Nature Materials , 14(3):337–342, 2015.[14] V. Hakim and W-J. Rappel. Dynamics of the globally coupled complex Ginzburg-Landau equation.
Physical Review A , 46(12):R7347, 1992.[15] J. F. Ravoux, S. Le Dizes, and P. Le Gal. Stability analysis of plane wave solutions ofthe discrete Ginzburg-Landau equation.
Physical Review E , 61(1):390, 2000.[16] K-I. Maruno, A. Ankiewicz, and N. Akhmediev. Exact localized and periodic solutionsof the discrete complex Ginzburg–Landau equation.
Optics Communications , 221(1-3):199–209, 2003.[17] S. J. Thomson, M. M. P. Couchman, and J. W. M. Bush. Collective vibrations ofconfined levitating droplets.
Physical Review Fluids , 5:083601, Aug 2020.[18] S. J. Thomson, M. Durey, and R. R. Rosales. Collective vibrations of a hydrodynamicactive lattice.
Proceedings of the Royal Society A , 476(2239), 2020.[19] G. C Sethia and A. Sen. Chimera states: the existence criteria revisited.
PhysicalReview Letters , 112(14):144101, 2014. 3520] B. Thakur and A. Sen. Collective dynamics of globally delay-coupled complexGinzburg-Landau oscillators.
Chaos: an Interdisciplinary Journal of Nonlinear Sci-ence , 29(5):053104, 2019.[21] Y. Couder, S. Protiere, E. Fort, and A. Boudaoud. Walking and orbiting droplets.
Nature , 437(7056):208–208, 2005.[22] J. W. M. Bush. Pilot-wave hydrodynamics.
Annual Review of Fluid Mechanics , 47:269–292, 2015.[23] J. W. M. Bush, Y. Couder, T. Gilet, P. A. Milewski, and A. Nachbin. Introduction tofocus issue on hydrodynamic quantum analogs.
Chaos: an Interdisciplinary Journal ofNonlinear Science , 28(9):096001, 2018.[24] A. Eddi, E. Sultan, J. Moukhtar, E. Fort, M. Rossi, and Y. Couder. Information storedin Faraday waves: the origin of a path memory.
Journal of Fluid Mechanics , 674:433–463, 2011.[25] J. Mol´aˇcek and J. W. M. Bush. Drops walking on a vibrating bath: towards a hydro-dynamic pilot-wave theory.
Journal of Fluid Mechanics , 727:612–647, 2013.[26] C. Bechinger, R. Di Leonardo, H. L¨owen, C. Reichhardt, G. Volpe, and G. Volpe.Active particles in complex and crowded environments.
Reviews of Modern Physics ,88(4):045006, 2016.[27] C. Scholz, S. Jahanshahi, A. Ldov, and H. L¨owen. Inertial delay of self-propelled par-ticles.
Nature Communications , 9(1):1–9, 2018.[28] M. X. Lim, A. Souslov, V. Vitelli, and H. M. Jaeger. Cluster formation by acousticforces and active fluctuations in levitated granular matter.
Nature Physics , 15(5):460–464, 2019.[29] H. L¨owen. Inertial effects of self-propelled particles: From active Brownian to activeLangevin motion.
The Journal of Chemical Physics , 152(4):040901, 2020.3630] A. U. Oza, R. R. Rosales, and J. W. M. Bush. A trajectory equation for walkingdroplets: hydrodynamic pilot-wave theory.
Journal of Fluid Mechanics , 737:552–570,2013.[31] S. H. Strogatz.
Nonlinear dynamics and chaos: with applications to physics, biology,chemistry, and engineering . CRC press, 2018.[32] J. K. Kevorkian and J. D. Cole.
Multiple scale and singular perturbation methods ,volume 114. Springer Science & Business Media, 2012.[33] C. R. Doering, J. D. Gibbon, D. D. Holm, and B. Nicolaenko. Low-dimensional be-haviour in the complex Ginzburg-Landau equation.
Nonlinearity , 1(2):279, 1988.[34] M. Bartuccelli, P. Constantin, C. R. Doering, J. D. Gibbon, and M. Gisself¨alt. Onthe possibility of soft and hard turbulence in the complex Ginzburg-Landau equation.
Physica D: Nonlinear Phenomena , 44(3):421–444, 1990.[35] P. Coullet and S. Fauve. Propagative phase dynamics for systems with Galilean invari-ance.
Physical Review Letters , 55(26):2857, 1985.[36] P. C. Matthews and S. M. Cox. Pattern formation with a conservation law.
Nonlinearity ,13(4):1293, 2000.[37] N. L. Komarova and A. C. Newell. Nonlinear dynamics of sand banks and sand waves.
Journal of Fluid Mechanics , 415:285–321, 2000.[38] S. M. Cox and P. C. Matthews. New instabilities in two-dimensional rotating convectionand magnetoconvection.
Physica D: Nonlinear Phenomena , 149(3):210–229, 2001.[39] T. B. Benjamin and J. E. Feir. The disintegration of wave trains on deep water.
Journalof Fluid Mechanics , 27(3):417–430, 1967.[40] J. T. Stuart and R. C. DiPrima. The Eckhaus and Benjamin-Feir resonance mechanisms.
Proceedings of the Royal Society of London. A: Mathematical and Physical Sciences ,362(1708):27–41, 1978. 3741] P. A. Milewski and E. G. Tabak. A pseudospectral procedure for the solution of nonlinearwave equations with examples from free-surface flows.
SIAM Journal on ScientificComputing , 21(3):1102–1114, 1999.[42] M. Durey, P. A. Milewski, and Z. Wang. Faraday pilot-wave dynamics in a circularcorral.
Journal of Fluid Mechanics , 891:A3, 2020.[43] D. M. Abrams and S. H. Strogatz. Chimera states for coupled oscillators.
PhysicalReview Letters , 93(17):174102, 2004.[44] G. C. Sethia, A. Sen, and F. M. Atay. Clustered chimera states in delay-coupled oscil-lator systems.
Physical Review Letters , 100(14):144102, 2008.[45] S. Nkomo, M. R. Tinsley, and K. Showalter. Chimera states in populations of nonlocallycoupled chemical oscillators.
Physical Review Letters , 110(24):244102, 2013.[46] E. A. Martens, S. Thutupalli, A. Fourri`ere, and O. Hallatschek. Chimera statesin mechanical oscillator networks.
Proceedings of the National Academy of Sciences ,110(26):10563–10567, 2013.[47] J. Wojewoda, K. Czolczynski, Y. Maistrenko, and T. Kapitaniak. The smallest chimerastate for coupled pendula.
Scientific Reports , 6(1):1–5, 2016.[48] J. F. Totz, J. Rode, M. R. Tinsley, K. Showalter, and H. Engel. Spiral wave chimerastates in large populations of coupled chemical oscillators.