Slow passage through the Busse balloon – predicting steps on the Eckhaus staircase

Abstract

Motivated by the impact of worsening climate conditions on vegetation patches, we study dynamic instabilities in an idealised Ginzburg–Landau model. Our main results predict time instances of sudden drops in wavenumber and the resulting target states. The changes in wavenumber correspond to the annihilation of individual vegetation patches when resources are scarce and cannot support the original number of patches. Drops happen well after the primary pattern has destabilised at the Eckhaus boundary and key to distinguishing between the disappearance of 1,2 or more patches during the drop are complex spatio-temporal resonances in the linearisation at the unstable pattern. We support our results with numerical simulations and expect our results to be conceptually applicable universally near the Eckhaus boundary, in particular in more realistic models.

1. Introduction

We study the evolution of spatially periodic patterns as system parameters vary slowly. Our motivation stems from ecological problems in which slowly varying parameters model worsening environmental conditions due to climate change. Of particular interest are dryland ecosystem models describing the interaction between vegetation density and available resources such as water. An important goal is to understand the process of desertification, predicting how much vegetation will remain as, for instance, average yearly rainfall decreases due to changing climate conditions. In ODE models that do not account for spatial variation of vegetation density or available resources, one often predicts tipping, in which a sudden collapse from a vegetated to a desert state occurs once the yearly rainfall decreases past some critical value. One also often observes hysteresis in such models, so that it is difficult to reverse desertification once it has occurred even if average rainfall begins to recover [49].

There has been great recent interest [2, 3, 9, 20, 40, 49, 50] in understanding the role of spatial dependence of the distribution of vegetation and resources on the ecosystem dynamics and eventual collapse. In particular, it was recently observed that in spatially extended models, a Turing instability of the uniformly vegetated state may occur prior to reaching the tipping point in the spatially independent system; see [49] and references therein. Uniform vegetation then gives way to periodic vegetation patches, and the pattern now evolves through the space of stable periodic patterns, called the Busse balloon; see [7, 8] and [14, Ch.IIIB2b]. In some cases, the system may avoid tipping completely, with stable periodic vegetation patterns persisting past the former tipping point. This ‘Turing before tipping’ mechanism then increases the resilience of the ecosystem as it avoids total collapse and the system may in principle recover by evolving back through the Busse balloon of stable patterns if environmental conditions begin to improve.

The purpose of this work is to better understand the evolution of patterns in and beyond the Busse balloon as system parameters vary slowly. One expects that the pattern will vary continuously, slightly adjusting the amplitude and shape of vegetation patches, until it reaches the boundary of the Busse balloon, at which point the pattern becomes unstable in the system with frozen parameters. In ecological terms, the available resources are no longer sufficient to support the present number of vegetation patches, so the system ‘sacrifices’ some number of vegetation patches, giving way to a new stable pattern with fewer vegetation patches. This self-organised sacrifice intuitively requires coordination between the different vegetation patches, all of which experience the same scarcity of resources. In parameter space, this crisis is represented by a sudden jump from the boundary to the interior of the Busse balloon, to a pattern with a lower wavenumber; see Figure 1 for an illustration.

Figure 1. Top: schematic depiction of vegetation patches in one- and two-dimensional domains and drops of density as parameters evolve slowly; bottom: schematic depiction of the Busse balloon near a Turing instability with parameter $\mu$ vertical, wavenumber horizontal. Wavenumbers of equilibria are quantised due to imposed spatially periodic boundary conditions with small (left) and large (right) period; see for example the discussion in (1.3). Patterns exist above $\mu _{\mathrm{ex}}$ but are stable only above $\mu _{\mathrm{eck}}$ . Red curves are schematic sample paths of observed wavenumbers when the parameter $\mu$ is slowly decreased, evolving in a staircase along the Eckhaus boundary. See also Figure 3 for numerical simulations.

While this picture of ‘bouncing off the boundary of the Busse balloon’ is somewhat intuitive and is born out in numerical simulations (see for instance Figure 5), there do not appear to be theoretical predictions of these dynamics of evolution through the Busse balloon in simple pattern-forming models.1 For instance, while it is intuitive that one must jump to another stable pattern once the original underlying pattern becomes unstable past the boundary of the Busse balloon, it is not clear how many vegetation patches will be lost in such a transition. A natural goal is then to answer the following question.

Given an initial periodic pattern with small fluctuations, can we predict the evolution of the number of patches as parameters vary slowly? In particular, can we predict the time at which the number of patches will change, and how much it will drop by?

We formulated our motivation in terms of dryland vegetation patches, but we believe that there are many settings in which this evolution of patterns with varying parameters is at the centre of the description of dynamics; see for instance studies of Turing patterns in growing domains [37] or the analysis of melting-boundary convection in [59].

The question embeds naturally into the larger set of problems known as delay of instability; see for instance [15, 27, 44, 54], in particular the early work in [43], and see [10, 13, 23, 26] for dynamic bifurcations in the context of pattern-forming instabilities. Such delays manifest themselves in contexts where a branch of equilibria undergoes a bifurcation and destabilises, or vanishes entirely. ‘Dynamically’ passing through this bifurcation diagram by slowly changing the parameter with rate $\varepsilon$ , one expects to follow the ‘static’ picture, that is, one expects to track the stable branch and, in the absence of the stable branch, a strong unstable manifold of the equilibrium at criticality; see Figure 2. In the simplest case of a saddle-node bifurcation, the solution stays near the remnant of the equilibria for a time $\varepsilon ^{-1/3}$ [38]. Similar results are available in the case of pitchfork and transcritical bifurcations when the parameter change includes a generic disturbance of the trivial branch [39]. If the trivial branch is undisturbed, the delay is of order $\varepsilon ^{-1}$ , that is, of order one in the parameter, with take-off point determined by the initial parameter value. Of course such a scenario is global in the sense that it requires assumptions on the system far away from the bifurcation point, at which point also other eigenvalues might be relevant; see for instance [33] for an analysis of two such eigenvalues interacting and [32] for a case with continuous spectrum crossing. More robust delays occur in Hopf bifurcations [27, 45, 46], where recent work, closer in spirit to our work here, also explores a PDE setting where many modes destabilise instantaneously [25, 34].

Figure 2. Delay of bifurcation in pitchfork (left) and saddle-node (right) bifurcation when the parameter is slowly varied. Left, top: slow passage through a subcritical pitchfork bifurcation, $a'=\mu a + a^3,\ \mu '=\varepsilon$ , yields $\mathrm{O}(1)$ delay in parameter space of the departure from the unstable state. Left, bottom: the reason for the delay in the pitchfork bifurcation is an accumulated exponential smallness from the dynamics in the stable regime, here shown in a schematic log plot of the amplitude. Right: in a slow passage near a fold, $a'=-\mu -a^2,\ \mu '=\varepsilon$ , the delay is small, $\mathrm{O}(\varepsilon ^{2/3})$ in the parameter. See text for details on how our results relate to the delay in the pitchfork bifurcation.

Our setting relates to the global aspect of the delay of instability in pitchfork bifurcations. In fact, the Eckhaus instability, which generically describes the boundary of the Busse balloon near a Turing bifurcation, manifests itself in bounded domains as a subcritical pitchfork bifurcation. With this analogy, we predict delays of order one in the parameter before leaving a fixed pattern. A simple linear calculation indeed correctly predicts the take-off point through an integral criterion: take-off points are found by requiring that the average of the marginally stable eigenvalue along part of the branch followed by the solution is zero. Intuitively, the cumulative exponential decay in the stable regime is ‘spent’ in the unstable regime where the solution grows exponentially, until both cancel and the solution reaches finite distance from the equilibrium branch again.

Beyond this simple linear calculation that we corroborate and supplement with global information, our key insight here is that the global nature of the delay makes it necessary to take into account other potentially unstable modes. Our results in this direction answer the main question posed above through linear predictions based on the presence of nonlinear spatio-temporal resonances. We predict, depending on the initial parameter value, the time when the solution leaves a neighbourhood of the equilibrium and the new equilibrium that the solution approaches. Since in fact the number of maxima (or vegetation patches, in this particular scenario) drops, we refer to this critical time as the drop time, to the associated parameter value as the drop parameter value, and to the decrease in the number of maxima as the drop number. Our results predict drop numbers in a range from 1 to 4 and associated drop times. The predictions are based on spatio-temporal resonances in the complex plane and are reminiscent of criteria for front invasion speeds and transitions from convective to absolute instability in large bounded domains as explored in [1, 19] and point to a potentially more comprehensive theory of spatio-temporal instabilities in extended systems. We emphasise that the results here address one particular, albeit ubiquitous example of instability, the long-wavelength Eckhaus destabilisation. Boundaries of the Busse balloon in systems far from equilibrium often involve different instability modes, in particular period-doublings; see [16, 17] for examples and see [48] for a conceptual view on instabilities of Turing patterns in one-dimensional systems.

To be specific, we study this question here in the Ginzburg–Landau equation

(1.1) \begin{align} A_t = A_{xx} + \mu (\varepsilon t) A - A |A|^2, \quad A = A(x,t) \in \mathbb{C}, \quad x\in \mathbb{R}, \quad t \gt 0. \end{align}

The Ginzburg–Landau equation is one of the simplest models of pattern-forming systems, and a universal amplitude equation that describes spatio-temporal dynamics near a Turing instability [14, 47]. It is therefore a natural starting point to understand the evolution of periodic patterns through the Busse balloon. In the Ginzburg–Landau equation, the boundary of the Busse balloon is determined by the Eckhaus instability; see [14, iVA1a(ii)] for background, [36, 55] for a study of the dynamics of the instability, [29, 56] for the effect of noise, [51] for finite-size effects, and Section 2 below for a basic review. We allow the linear coefficient $\mu$ to vary slowly in time, with time scale $\varepsilon^{-1}$ , and consider the problem on bounded domains with periodic boundary conditions. We will often make comparisons to the static problem, where $\mu$ does not vary in time. We emphasise here that when interpreting results obtained for the Ginzburg–Landau equation in the context of a Turing instability, one needs to account for the fact that $A(x,t)$ is an amplitude-phase modulation, so that actual solutions are of the form $A(\varepsilon x,\varepsilon ^2 t)\mathrm{e}^{\mathrm{i} x}$ , say. In particular, the preferred constant state $A\equiv const$ in Ginzburg–Landau corresponds to a preferred patterned state in the original model.

The equation (1.1) with $\varepsilon =0$ , $\mu (\varepsilon t)=\mu \gt 0$ fixed possesses equilibria $A_{*}(x;j)=\sqrt{\mu -j^2}\mathrm{e}^{\mathrm{i} jx}$ for all $|j|\leqslant \sqrt{\mu }$ . We are interested in the fate of these ‘patterned’ equilibria for a fixed $j$ when the parameter $\mu$ is slowly decreased and, as it turns out, the equilibrium is destabilised. We wish to consider instabilities on a fixed-size domain which we assume, after possibly rescaling $x,t$ , and $|A|$ , to be $x\in [0,2\pi ]$ . To accommodate an initially fixed equilibrium with arbitrary $j$ we impose $2\pi$ -twist-periodic boundary conditions, relative to this equilibrium,

(1.2) \begin{align} A(2 \pi, t) = \mathrm{e}^{2\pi \mathrm{i} j} A(0, t), \quad A_x(2 \pi, t) = \mathrm{e}^{2 \pi \mathrm{i} j} A_x (0, t). \end{align}

This allows us to treat $j$ as a continuous rather than a quantised parameter.

Next, defining a new variable $B$ through $A(x,t) = \mathrm{e}^{\mathrm{i} j x} B(x,t)$ , we find that $B$ solves the conjugated equation

(1.3) \begin{align} B_t = (\partial _x + \mathrm{i} j)^2 B + \mu (\varepsilon t) B - B |B|^2, \end{align}

with periodic boundary conditions

(1.4) \begin{align} B(2 \pi, t) = B(0, t), \quad B_x (2 \pi, t) = B_x (0, t). \end{align}

This equation has a family of constant solutions

(1.5) \begin{equation} E_j=\{\mathrm{e}^{\mathrm{i}\varphi }B_*(j),\varphi \in [0,2\pi )\},\quad B_*(j) =\sqrt{\mu - j^2}, \end{equation}

modulated solutions

(1.6) \begin{equation} E_{j'}=\{\mathrm{e}^{\mathrm{i}\varphi }B_*(j') \mathrm{e}^{i(j'-j)x},\varphi \in [0,2\pi )\},\quad B_*(j') =\sqrt{\mu - j'^2}, \end{equation}

and the trivial solution $B_*=0$ . Drop numbers as described above, when, say, $j\gt 0$ , are now simply the differences $j-j'$ for a trajectory that hovers subsequently first near an equilibrium $E_j$ and later near an equilibrium $E_{j'}$ with $j\gt |j'|$ . Together with the drop number, we also wish to predict the drop time $t_*$ , when a solution leaves a small fixed neighbourhood of $E_{j}$ , or, more naturally, the associated drop parameter value $\mu (\varepsilon t_*)$ . The main insight of the analysis presented in this paper can be roughly described as follows.

Main results. Consider solutions to (1.3) with periodic boundary conditions (1.4), and with initial conditions close to $E_j$ for some $j$ , and assume that $E_j$ is initially stable for $\mu = \mu _0$ . Under conceptual assumptions on the existence of connecting orbits in the static problem, drop times and drop numbers as $\mu$ varies are determined by resonance criteria for eigenvalues, and associated integral conditions summarised in Section 3.1.

As a consequence, our results predict the drops in the intricate web of winding numbers $j$ observed over time, when only the initial parameter value is varied; see Figure 3 for this increase of the drop with distance from criticality and for the ensuing sequence of drop events.

Figure 3. Simulations of (1.1) with initial condition $A_*(x;j)$ , $j=8$ and $\varepsilon =0.05$ , with initial $\mu$ -values $\mu _{\mathrm{in}}=\mu _{\mathrm{c}}+\hat{\mu }$ , $\hat{\mu }=6,16,\ldots,56$ . Plotted are winding numbers of $A$ as time progresses and $\mu$ decreases, over the parameter value $\mu$ at a given time instance. Trajectories crisscross the Eckhaus boundary in a staircase pattern that depends on the initial parameter value. The value $\hat{\mu }=0$ corresponds to the onset of the Eckhaus instability at $j=8$ . Initial conditions further away from the instability lead to later drops and higher drop numbers. Fractions of 1 are added to show all itineraries simultaneously, so that the actual current winding number is the largest integer below the plotted curve. Also shown are the boundary of instability (red) and the boundary of existence of equilibria (black) with a given $j$ ; see Section 2 for details. Note that the figure is reflected along the diagonal when compared to Figure 1, using the traditional bifurcation-theoretic depiction of phase-space versus parameter also used in Figure 2.

2. The Eckhaus instability in the static problem

In this section, we recall stability of periodic patterns in the static problem

(2.1) \begin{align} B_t = (\partial _x + \mathrm{i} j)^2 B + \mu B - B |B|^2, \end{align}

with $\mu$ fixed and with $2\pi$ -periodic boundary conditions. Dynamics of (2.1) are influenced by its symmetries, which we now list for reference:

• • (spatial translation) $T_{x_0}\,:\, B(x,t) \mapsto B(x+x_0, t)$ , viewing $x$ on the circle ${\mathbb{R}}/ 2 \pi \mathbb{Z}$ ;

• • (reflection) $T_-\,:\, B(x,t) \mapsto \mathrm{e}^{-2\mathrm{i} jx} B({-}x,t)$ ;

• • (gauge rotation) $R_{\phi }\,:\, B(x,t) \mapsto \mathrm{e}^{i\phi } B(x,t)$ for $\phi \in \mathbb{R}/2\pi \mathbb{Z}$ ;

• • (complex reflection) $S_-\,:\,B(x,t)\mapsto \mathrm{e}^{-2\mathrm{i} jx}\bar{B}(x,t)$ .

From the discussion above, we recall that we are interested in the stability properties of the spatially constant equilibria $E_{j}$ , defined in (1.5). We remark that one commonly studies stability properties on the unbounded real line rather than on a fixed, possibly large finite domain. Separating into phase and amplitude dynamics and linearising, one there finds that the $E_{j}$ are linearly stable provided $\mu \gt 3 j^2$ , and the long wavelength instability at $\mu = 3 j^2$ is referred to as the Eckhaus instability; see for instance [31, Chapter 8.2] for a review. Our focus on finite domains here gives well-known corrections to this threshold. The restriction to finite domains here appears essential to the techniques used below.

2.1 Local stability analysis

The Eckhaus instability in finite domains was analysed in [57], and we briefly review the basic linear analysis here. Separating (2.1) into real and imaginary parts $B = u + iv$ and linearising about $(u,v) = (B_*(j), 0)$ , we find

\begin{align*} u_t &= u_{xx} - 2 j v_x - 2 (\mu - j^2) u \\[2pt] v_t &= v_{xx} + 2 j u_x. \end{align*}

The linearisation is block diagonal in Fourier modes $ u(x) = \sum _{\ell \in \mathbb{Z}} u_\ell \mathrm{e}^{\mathrm{i} \ell x},\ v(x) = \sum _{\ell \in \mathbb{Z}} v_\ell \mathrm{e}^{\mathrm{i} \ell x}$ ,

(2.2) \begin{align} \begin{pmatrix} \dot{u_\ell } \\[5pt] \dot{v_\ell } \end{pmatrix} = \begin{pmatrix} -\ell ^2 - 2( \mu -j^2) & - 2 \mathrm{i} j \ell \\[5pt] 2 \mathrm{i} j \ell & -\ell ^2 \end{pmatrix} \begin{pmatrix} u_\ell \\[5pt] v_\ell \end{pmatrix}. \end{align}

This matrix is self-adjoint, a reminder that the Ginzburg–Landau equation is a gradient flow. It has a negative trace, so that it always has at least one negative real eigenvalue. From a short calculation, one finds that the other eigenvalue is positive when

(2.3) \begin{align} \mu \lt \mu _{\ell,\mathrm{c}} \,:\!=\, 3j^2 - \frac{\ell ^2}{2} \end{align}

for $\ell \neq 0$ , while the $\ell = 0$ block always has one zero and one negative eigenvalue. As $\mu$ decreases, the instability is therefore always first seen in the $\ell = \pm 1$ mode. The eigenvalue that becomes unstable past $\mu _{\ell,\mathrm{c}}$ is given by

(2.4) \begin{align} \lambda _{\ell,+}(\mu ) = - \ell ^2 - (\mu -j^2) + \sqrt{4\ell ^2 j^2 + (\mu -j^2)^2}. \end{align}

For notational simplicity, we usually suppress the dependence of $\lambda _{\ell,+}(\mu )$ and $\mu _{\ell,\mathrm{c}}$ on the base wavenumber $j$ . To further simplify notation, we write $\lambda _\ell$ instead of $\lambda _{\ell,+}$ since the eigenvalue $\lambda _{\ell,-}$ will be irrelevant for the discussion in this work. We also write simply $\mu _{\mathrm{c}}\,:\!=\,\mu _{\mathrm{1,c}}$ for the boundary of stability, noting that modes with $|\ell | \neq 1$ destabilise for smaller values of $\mu$ , reflecting the side-band nature of the Eckhaus instability. Eigenvectors are

\begin{equation*} (u_{\ell },v_{\ell })(\mu )=\left (j^2-\mu +\sqrt {4\ell ^2 j+(\mu -j^2)^2},2\mathrm {i}\ell j\right ). \end{equation*}

Note that, since the equilibrium is fixed under translations $T_{x_0}$ and the reflection $T_-S_-: B(x)\mapsto \bar{B}({-}x)$ , the linearisation commutes with these symmetries which henceforth act on eigenspaces (and on the flows on invariant manifolds tangent to those eigenspaces). Since the algebraic expression for eigenvalues and eigenvectors are somewhat cumbersome, we will later verify assumptions explicitly in the limit of large $j$ , where we set $\mu =\mu _{\mathrm{c}}+\hat{\mu }$ , and find

(2.5) \begin{align} \lambda _{\ell }(\mu )=\frac{-\ell ^2}{2}\left (\frac{\ell ^2}{2}-\frac{1}{2}+\hat{\mu }\right )j^{-2}+\mathrm{O}\left (j^{-4}\right ),\qquad (u_{\ell },v_{\ell })(\mu )=\left (\frac{\ell }{2}+\mathrm{O}\left (j^{-2}\right ),\mathrm{i} j\right ), \end{align}

also writing $e_{\ell }(\mu )$ short for the normalised eigenvector; see Figure 4 for an illustration of the dependence of eigenvalues on $\mu$ . The stability threshold $\hat{\mu }=0$ corresponds to the Eckhaus boundary $\mu = 3 j^2$ in the infinite domain, but with the finite size correction $-\frac{1}{2}$ . Note that since $u$ and $v$ are real, we have $u_\ell = \overline{u_{-\ell }}, v_\ell = \overline{v_{-\ell }}$ , and in particular $\lambda _{\ell }(\mu ) = \lambda _{-\ell }(\mu ).$

Figure 4. Eigenvalues of the linearisation at $E_j$ with $j=8$ (top) and $j=4$ (bottom), for parameter values $\mu =j^2$ at onset to $\mu _{\mathrm{c}}+10$ . Enlargement near criticality (right) reveals the intricate crossovers that are explicit in the limit of large $j$ . Eigenvalues are shown from $\ell =1$ in dark blue to $\ell =7$ (top) and $\ell =5$ (bottom), respectively, in light blue.

2.2 Connecting orbits

To understand the evolution of periodic patterns in the dynamic problem, we are interested not only in stability boundaries in the static problem but also in connecting orbits. That is, we ask the following question: if we start with an initial condition $B_0$ near the equilibrium $E_j$ which is Eckhaus unstable, where does the solution $B(t)$ go as $t \to \infty$ ?

The complex Ginzburg–Landau equation is in fact a gradient flow in $L^2([0,2\pi ),\mathbb{C})$ ,

(2.6) \begin{equation} B_t=-\nabla _{L^2}\mathcal{E}[B],\qquad \text{ with energy }\quad \mathcal{E}[B] = \frac{1}{2 \pi }\int _0^{2 \pi } \left (\frac{1}{2} |B_x+\mathrm{i} jB|^2 - \frac{\mu }{2} |B|^2 + \frac{1}{4} |B|^4 \right )\, \mathrm{d} x. \end{equation}

As a consequence, the solution $B(t)$ necessarily converges to an equilibrium with energy lower than $E_j$ . A direct calculation shows that from this perspective, out of all equilibria $E_k$ , only wavenumbers $|k| \lt j$ are eligible, and one may suspect that solutions $B(t)$ originating near $E_j$ converge to a specific $E_k$ . Unfortunately, it seems difficult to establish analytically which equilibrium $E_k$ is selected in this sense. In fact, there are also spatially patterned equilibria other than the $E_k$ , believed to all be unstable as they limit, in a large-domain limit, onto unstable phase defect solutions such as

(2.7) \begin{equation} A_{\mathrm{d}}(x;\,k)=\left (\sqrt{2}k+\mathrm{i} \sqrt{1-3k^2}\tanh (\sqrt{1-3k^2}x/\sqrt{2})\right )\mathrm{e}^{\mathrm{i} k x}; \end{equation}

see [42]. Some, for our purposes rather restrictive, results have been obtained in [18], where heteroclinic orbits between two distinct equilibrium sets $E_j$ and $E_{j-k}$ have been constructed. The arguments there are based on the variational structure (2.6), roughly demonstrating that there is in fact a unique candidate for connecting orbits with given co-periodicity. In our contexts, the results only give conclusive information in the case $\frac{1}{2} \lt j \lt 1$ , in which case a generic perturbation of $E_j$ will converge to $E_{j-1}$ .

On the other hand, the fate of perturbations is rather straightforward to establish numerically. We distill our findings in the following conceptual assumption about the static problem, which we will later rely on to draw conclusions in the problem with slowly varying parameters.

First note that for a generic $\mu$ , the eigenvalues $\lambda _{\ell } (\mu )$ will be distinct up to the symmetry $\lambda _{\ell } (\mu ) = \lambda _{-\ell } (\mu )$ . The special crossover points at which this condition is violated, that is, at which eigenvalues for different wavenumbers $\ell$ collide, will play a role later but are excluded in the following discussion. We denote the crossover point, at which $\lambda _{\ell }(\mu ) = \lambda _{k}(\mu )$ , by $\mu = \mu _{\ell \,:\, k}$ . Away from these crossover points, we let $\mathrm{e}^{\mathrm{u}}_\ell$ denote the real two-dimensional eigenspace associated with the unstable eigenvalues $\lambda _{\pm \ell }(\mu )$ .

The following hypothesis is sometimes colloquially referred to as node conservation or linear mode selection.

Hypothesis 2.1 (Linear mode selection). For $0\lt \delta _1\ll \delta _0\ll 1$ sufficiently small fixed, trajectories with initial conditions given as a perturbation of $E_{j}$ of size $\delta _0$ will limit on $E_{j-\ell _*}$ if their amplitude is maximal in the unstable eigenspace with index $\ell _*$ . To be precise, writing the projection of the initial condition on the unstable subspace of $E_j$ as a linear combination $\sum \alpha _\ell e_{\ell } \mathrm{e}^{\mathrm{i}\ell x}$ , we have $\delta _0=|\alpha _{\ell _*}|\gg{\delta _1}=\max _{\ell \neq \ell _*}|\alpha _\ell |$ . Here, $e_{\ell }$ is the normalised eigenvector to the eigenvalue $\lambda _{\ell }$ from (2.5); see Figure 6 for a schematic.

To explain the basic intuition, notice that perturbations of the constant equilibrium in the direction of the eigenvector $\lambda _{\ell }$ , that is perturbations where $\delta _1=0$ , have a sinusoidal shape and winding number $\ell$ relative to $E_j$ . The hypothesis states that as the amplitude grows, this ‘modal’ shape is preserved and terminates in the equilibrium $E_{j-\ell }$ , with winding number $-\ell$ relative to the origin.

Figure 5. Left: existence and stability boundaries, as well as spatio-temporal resonances responsible for cross-overs in the drop number changes, together with mode number of a sample trajectory. Resonance curves $\mu _{1^2:2}$ , $\mu _{1,2:3}$ , etc. correspond to parameter values where $\lambda _1+\lambda _1=\lambda _2$ , $\lambda _1+\lambda _2=\lambda _3$ , and so forth. In the region between solid and dashed part of the curves, $\lambda _{j+1}\gt \lambda _1+\lambda _j$ ; see Section 3.1 for details on resonances and their relevance. Right: space-time plot of $\mathrm{Re}(A)$ of the same sample trajectory, with time axis represented in terms of the slowly varying $\mu$ , varying up until $\mu \lt 0$ and no nontrivial equilibria exist.

Figure 6. Schematic of the picture implied by Hypothesis 2.1: the unstable manifold of $E_j$ , here 2-dimensional, contains open sets of orbits that connect to $E_{j-1}$ and $E_{j-2}$ , respectively. The hypothesis guarantees that the connecting orbits to $E_{j-1}$ and $E_{j-2}$ , respectively, contain the caps of conical regions around the respective eigenspaces. In particular, trajectories in the boundary between connecting orbits to $E_{j-1}$ and $E_{j-2}$ are not arbitrarily close to the eigenspaces. It seems plausible that the boundary of validity of the hypothesis is related to trajectories on a codimension-one stable manifold of saddle equilibria, here referred to as ‘defect’, since equilibria of the form (2.7) are natural candidates for such a role in the setting of an unbounded domain.

Figure 7. Testing the linear prediction hypothesis in the fixed- $\mu$ problem, we computed trajectories with small perturbations of an unstable equilibrium $E_j$ in the modes $\ell =1$ (left) and $\ell =2$ (right). We found that the trajectories converged to the equilibrium $E_{j-\ell }$ in a range (shaded grey) far exceeding the possible drop ranges investigated below. In the left figure, $\ell =1$ , drops to $E_{j-1}$ were confirmed up for all $\mu$ larger than and up to a critical value $\mu _{1,\mathrm{bdy}}$ , which includes a region where $\lambda _2\gt \lambda _1$ (shown as $\mu \gt \mu _{1:2}$ ) and a region where $\lambda _2\gt 2\lambda _1$ (shown as $\mu \gt \mu _{1^2:2}$ ). However, it does not encompass the entire parameter region where $E_{j-\ell }$ is stable (shown as $\mu \gt \mu _{\mathrm{stab}}^{-1}$ ). Remarkably, higher drops appear to occur past a somewhat fixed distance to the instability boundary $\mu _{\mathrm{stab}}$ . On the right, the region with consistent drop-by-2 (shaded grey) also encompasses well the region where we predict a local drop-by-two. We refer to (3.4) for precise definitions of the various resonance curves shown and the main prediction in Section 3.1 relying on these resonances to predict drop numbers and drop times.

Hypothesis 2.1 is reasonably well supported by numerical experiments; see Figure 7 for experiments in the cases $\ell =1$ and $\ell =2$ . In the numerical experiments, we perturbed the unstable equilibrium in the eigenvector to $\lambda _{\ell }$ , only, thus checking the milder version of the hypothesis with $\delta _1=0$ . We found consistent drops by $\ell$ in a large region (shaded grey in the figure). We superimposed this region with eigenvalue resonance curves that we computed using numerical continuation. In the left panel, we show the transition from drop-by-one to drop-by-two as predicted in the next section at $\mu _{1^2:2}$ . This resonance occurs well inside the shaded grey region. Note that for smaller $\mu$ the resonance is not present, the numerical continuation of the resonance condition undergoes a saddle-node with a second resonance curve $\mu ^*_{1^2:2}$ . In the right panel, we show in particular the resonance $\mu _{1,2:3}$ , which indicates the transition to a drop-by-three, again well inside the region where the hypothesis holds for $\ell =2$ ; for details, compare the definition of the log amplitudes $b_j$ (3.5)–(3.8) and the main prediction (3.10).

In the spirit of this assumption, we fix $\delta \gt 0$ sufficiently small and consider a trajectory that starts in a neighbourhood of the unstable equilibrium $E_j$ . We define the local drop number (depending on $\delta$ ) at the time when the distance to the $E_j$ reaches $\delta$ by finding the index of the unstable eigenspace with maximal amplitude in the sense of Hypothesis 2.1. In the results presented in the next section, it then turns out that the distance to the eigenspace is in fact very small, that is, $\delta _1\ll \delta _0$ .

We say that the global drop number of such a trajectory is $\ell$ if $E_{j-\ell }$ is the first equilibrium so that the trajectory reaches its neighbourhood of size $\delta$ . Hypothesis 2.1 is geared to imply that local drop numbers and global drop numbers are equal in the static ( $\mu$ fixed) problem. In fact, local and global drop numbers agree well also in the dynamical problem, where $\mu$ varies slowly; compare Figure 11.

We note that, thinking about proving a version of Hypothesis 2.1, detailed information on connecting orbits usually requires information in addition to energy levels of equilibria or more general topological information [41]. One such property are comparison principles with the associated refined Sturm oscillation properties in the real scalar version of (1.1), known as the Allen-Cahn equation, $u_t = u_{xx} + \mu u - u^3$ . Here, in fact the structure of connecting orbits is completely determined; see for instance [28, Chapter 5.3], the review [22], or [21] for a more recent perspective. Some extensions towards systems of the type considered here were investigated in [6].

3. Slow passage through the Eckhaus instability

Based on our understanding of the static problem, in particular Hypothesis 2.1, we now turn to predictions for the dynamic problem

(3.1) \begin{align} B_t &= (\partial _x+\mathrm{i} j)^2 B + \mu B - B|B|^2, \end{align}

(3.2) \begin{align} \dot{\mu } &= - \varepsilon. \end{align}

Throughout, we envision an initial condition given as a small generic perturbation of the $j$ -modal equilibrium $E_j$ , for a $\mu$ -value initially large enough such that $E_j$ is stable. As time evolves, the solution will follow the $\mu$ -dependent equilibrium very closely as long as $\mu \gt \mu _{\mathrm{c}}$ . Our goal is to predict when (i.e., for which $\mu$ -value less than $\mu _{\mathrm{c}}$ ) the solution will leave a small neighbourhood of $E_j$ and what the dominant linear mode $\ell$ is at that time instance, hence predicting a drop-by- $\ell$ based on our linear mode selection assumption, Hypothesis 2.1.

3.1 Main result – spatio-temporal resonances and drop criteria

We summarise here the results of the analysis in the subsequent sections in the form of an algorithm that predicts, at leading order, local drop parameter values $\mu _{\mathrm{out}}$ and local drop number $m$ given an initial parameter value $\mu _{\mathrm{in}}$ and initial wavenumber $j$ in (3.2). Recall the eigenvalues of the linearisation $\lambda _{\ell }(\mu )$ and define resonant parameter values $\mu =\mu _{k_1,\ldots,k_s\,:\,k_0}$ as the values of $\mu$ where a resonance condition in the eigenvalues holds:

(3.3) \begin{align} \lambda _{k_0}(\mu )&=\sum _{\sigma =1}^s \lambda _{k_\sigma }(\mu ) \quad \text{ at } \quad \mu =\mu _{k_1,\ldots,k_s\,:\,k_0}, \end{align}

(3.4) \begin{align} k_0&=\sum _{\sigma =1}^s k_\sigma. \end{align}

For instance,

\begin{equation*} \lambda _{3}(\mu )=\lambda _{1}(\mu )+\lambda _{2}(\mu ),\quad \text { at } \quad \mu =\mu _{1,2:3}, \end{equation*}

and introducing short-hand notation,

\begin{equation*} \lambda _{3}(\mu )=3\lambda _{1}(\mu )\quad \text { at } \quad \mu =\mu _{1,1,1:3}\,=\!:\,\mu _{1^3:3}. \end{equation*}

The parameter values $\mu _{k_1,\ldots,k_s\,:\,k_0}$ represent spatio-temporal resonances in the sense that at these parameter values linear solutions $\mathrm{e}^{\mathrm{i} kx+\lambda t}$ exist that are resonant both in time (3.3) and space (3.4).

Next, we recursively define $\log$ -amplitudes of modes $b_\ell$ given $\mu _{\mathrm{in}}$ as follows:

(3.5) \begin{align} b_1(\mu )&=\int _{\mu }^{\mu _{\mathrm{in}}}\lambda _1(\hat \mu )\mathrm{d}\hat \mu, \end{align}

(3.6) \begin{align} b_2(\mu )&=\int _{\mu }^{\mu _{1^2:2}}\lambda _2(\hat \mu )\mathrm{d}\hat \mu + 2 b_1(\mu _{1^2:2}), \qquad \mu \lt \mu _{1^2:2} \end{align}

(3.7) \begin{align} b_3(\mu )&=\int _{\mu }^{\mu _{1,2:3}}\lambda _3(\hat \mu )\mathrm{d}\hat \mu + b_1(\mu _{1,2:3})+b_2(\mu _{1,2:3}), \qquad \mu \lt \mu _{1,2:3} \end{align}

(3.8) \begin{align} b_4(\mu )&=\int _{\mu }^{\mu _{2,2:4}}\lambda _4(\hat \mu )\mathrm{d}\hat \mu + 2b_2(\mu _{2^2:4}), \qquad \mu \lt \mu _{2,2:4}. \end{align}

Main prediction – drops $m\leqslant 4$ . Fix $\ell$ and an initial parameter value $\mu _{\mathrm{in}}\gt \mu _{\mathrm{c}}$ and consider dynamics with slowly decreasing $\mu =\mu _{\mathrm{in}}-\varepsilon t$ . The drop parameter value $\mu _{\mathrm{out}}$ is given to zeroth order in $\varepsilon$ by

(3.9) \begin{equation} \mu _{\mathrm{out}}=\mathrm{argmax}_{\mu \lt \mu _{\mathrm{c}}} \mathrm{max}_{1\leqslant \ell \leqslant 4 } b_\ell (\mu )=0, \end{equation}

that is, the largest (dynamically the first) $\mu$ -value at which any of the log amplitudes return to 0. The local drop number is

(3.10) \begin{equation} m= \mathrm{argmax}_{1\leqslant \ell }\, b_\ell (\mu _{\mathrm{out}}), \end{equation}

that is, the index $\ell$ that achieves $b_\ell (\mu _{\mathrm{out}})=0.$ The maxima are taken only over log amplitudes that are defined for the given parameter value $\mu$ . As a function of $\mu _{\mathrm{in}}$ , $\mu _{\mathrm{out}}\lt \mu _{\mathrm{c}}$ is strictly decreasing and $m$ is strictly increasing.

The $\log$ -amplitudes before the resonances can be defined as well, but they are always much smaller than the $\log$ -amplitudes with smaller index.

The reader will notice that only specific resonances occur in the definition of the $b_\ell$ . We found these resonances to be relevant in the sense that their effects dominate the effects of all other resonances in the examples we studied, in particular in the limit of large $j$ . We do not know of a general rule to predict which resonances dominate for larger potential drop numbers, so the algorithm for drops by $5$ or higher is more complex. For the mode $b_5$ , one needs to find $\log$ -amplitudes depending on different spatio-temporal resonances and maximise over all of those, for instance defining

\begin{align*} b_5^{2,3:5}(\mu )&=\int _{\mu }^{\mu _{2,3:5}}\lambda _5(\hat \mu )\mathrm{d}\hat \mu + b_2(\mu _{2,3:5}) + b_3(\mu _{2,3:5}) \qquad \mu \lt \mu _{2,3:5},\\[5pt] b_5^{1,4:5}(\mu )&=\int _{\mu }^{\mu _{1,4:5}}\lambda _5(\hat \mu )\mathrm{d}\hat \mu + b_1(\mu _{1,4:5}) + b_4(\mu _{1,4:5}) \qquad \mu \lt \mu _{1,4:5},\\[5pt] b_5^{1,2^2:5}(\mu )&=\int _{\mu }^{\mu _{1,2^2:5}}\lambda _5(\hat \mu )\mathrm{d}\hat \mu + b_1(\mu _{1,2^2:5}) + 2b_2(\mu _{1,2^2:5}) \qquad \mu \lt \mu _{1,2^2:5},\qquad \text{etc.} \end{align*}

With this information, we can define critical parameter values where changes in drop numbers occur. At these parameter values, the argmax in (3.10) is realised simultaneously at two different values of $\ell$ . In all situations we encountered, those two values of $\ell$ were adjacent, $\ell \in \{m,m+1\}$ and the drop number increases by one as $\mu _{\mathrm{in}}$ increases through this critical parameter value. We then denote this critical parameter value by $\mu _{\mathrm{in},m\to m+1}$ , with associated local drop parameter values $\mu _{\mathrm{out},m\to m+1}$ .

It turns out that $\mu _{\mathrm{out},1\to 2}=\mu _{1^2:2}$ , but for $m\geqslant 2$ , $\mu _{\mathrm{out},m\to m+1}$ is strictly less than the resonances associated with the mode $m+1$ , for instance $\mu _{1,m:m+1}$ , $\mu _{1^2,m-1:m+1}$ , etc. We will see that, as a consequence, corrections to the transition points are $\mathrm{O}(\sqrt{\varepsilon })$ for $\mu _{\mathrm{out},1\to 2}$ but $\mathrm{O}(\varepsilon )$ for other transitions.

We refer to Figures 8, 10, and 11 for numerical demonstration of the validity of the results presented here. The remainder of this section presents the rationale behind these predictions, explaining to what extent one should expect rigorous results, and some details on $\varepsilon$ -corrections as well as asymptotics for $j\to \infty$ . The numerical results also show that global drop numbers agree well with local drop number predictions. This is also confirmed by the direct testing of Hypothesis 2.1, Figure 7, which shows that the hypothesis holds well beyond the region where drop number changes are expected, at least for drop numbers of 1 and 2.

3.2 Drop-by-1 transitions

First, we consider initial conditions for which $B_0$ is close to an equilibrium $E_j$ , and the initial parameter value $\mu _{\mathrm{in}}$ is just above the stability boundary $\mu _{\mathrm{c}}=\mu _{1,\mathrm{c}}$ at which $\lambda _{1}(\mu )$ becomes positive. Starting very close to the stability boundary, we expect the first drop in the winding number to occur while $\lambda _{1}(\mu )$ is the only unstable eigenvalue, and so we study the reduced problem on the associated centre manifold. In this case, no further assumptions are necessary to make precise and rigorous statements about the local drop time.

When $\varepsilon = 0$ and $\mu =\mu _{1,\mathrm{c}}$ , the static problem has a complex one-dimensional, real two-dimensional centre manifold. A parameter-dependent centre manifold reduction [57] gives the leading order equation for the complex amplitude $a_1$ near $\mu = \mu _{1,\mathrm{c}}$ ,

(3.11) \begin{align} \dot{a}_1 = \lambda _{1} (\mu ) a_1 + c_1(\mu ) a_1 |a_1|^2, \quad c_1(\mu ) \gt 0, \end{align}

which is the normal form for a subcritical pitchfork bifurcation with rotational symmetry. The cubic coefficient, obtained through a centre-manifold expansion was shown to be positive in [57]. Solutions to (2.1) on this centre manifold are then recovered using the eigenvectors $e_1(\mu )$ (2.5) through

(3.12) \begin{align} B(x,t) = B_*(j) + u + \mathrm{i} v,\qquad u + \mathrm{i} v = a_1(t) e_{1}(\mu ) \mathrm{e}^{\mathrm{i} x} + \bar{a}_1(t) \bar{e}_{1}(\mu ) \mathrm{e}^{-\mathrm{i} x} + \mathrm{O}\left (|a_1|^2\right ). \end{align}

Considering these centre manifolds in the system (3.1)–(3.2) for $\varepsilon = 0$ , we find an invariant manifold, of real dimension three, given by the union over $\mu$ of the centre manifolds constructed above. This manifold is normally hyperbolic from the explicit splitting in the linearisation at the equilibria and hence persists for $\varepsilon$ small [4, 28], with dynamics given by

(3.13) \begin{align} \dot{a}_1 &= \lambda _{1} (\mu ) a_1 + c_1 (\mu ) a_1 |a_1|^2+\mathrm{O}\left (|a_1|^5\right ), \end{align}

(3.14) \begin{align} \dot{\mu } &= - \varepsilon. \end{align}

As noted above, the flow commutes with the action of translations $a_1\mapsto \mathrm{e}^{\mathrm{i}\varphi }a_1$ and reflections $a_1\mapsto \bar{a}_1$ , so that the real subspace is invariant and dynamics for arbitrary initial conditions are conjugate to the dynamics in the real space by complex rotation. Alternatively, one could also use a time-dependent centre-manifold reduction [11] after an appropriate modification of the $\mu$ -dynamics outside of the parameter regime considered and arrive at the same reduced equation.

We then consider a boundary-value problem where we prescribe the initial amplitude $a_{1,\mathrm{in}}=\delta \ll 1$ and initial parameter value $\mu _{\mathrm{in}}$ at a time $-T_{\mathrm{in}}$ , and seek time $T_{\mathrm{out}}$ and associated parameter value $\mu _{\mathrm{out}}$ where the amplitude $a_1$ is of size $\delta$ again,

(3.15) \begin{align} a_1({-}T_{\mathrm{in}}) &= a_{1,\mathrm{in}}=\delta,\qquad \mu ({-}T_{\mathrm{in}}) = \mu _{\mathrm{in}} \nonumber \\[5pt] a_1(T_{\mathrm{out}}) &= a_{1,\mathrm{out}} = \delta,\qquad \mu (T_{\mathrm{out}}) = \mu _{\mathrm{out}}. \end{align}

For convenience, we set

(3.16) \begin{equation} \mu _{\mathrm{in}}-\mu _{1,\mathrm{c}} = \varepsilon T_{\mathrm{in}}, \text{ so that } \mu _{1,\mathrm{c}}-\mu _{\mathrm{out}} = \varepsilon T_{\mathrm{out}},\quad \mu (0) =\mu _{1,\mathrm{c}}. \end{equation}

Because $\lambda _{1}(\mu _{\mathrm{in}})$ is stable, we expect that $a_1$ will initially decay exponentially, and then grow once $\mu$ decreases past $\mu _{1,\mathrm{c}}$ , so that the time $T_{\mathrm{out}}$ and thereby the value of $\mu _{\mathrm{out}}$ at which $a_1$ returns to its initial size, $a_1(T_{\mathrm{out}}) = \delta$ , are well defined.

Since we are only considering a time interval in which $a_1$ remains small, as a first approximation we analyse (3.13)–(3.14) by neglecting the nonlinear terms and solving the resulting linear equation explicitly for initial conditions $(a_{1,\mathrm{in}},\mu _{\mathrm{in}})$ , with solution

(3.17) \begin{align} a_1(t) = \exp \left ( \int _{-T_{\mathrm{in}}}^t \lambda _{1}(\mu (\varepsilon s)) \, \mathrm{d} s \right ) a_{1,\mathrm{in}}. \end{align}

Since $\lambda _{1}(\mu _{\mathrm{in}}) \lt 0$ , the integral is initially negative and begins increasing as $\mu$ decreases past $\mu _{1,\mathrm{c}}$ . Then in this linear prediction, the time $T_{\mathrm{out}}$ is the first time after $T_{\mathrm{in}}$ at which the equal area condition

(3.18) \begin{equation} \int _{-T_{\mathrm{in}}}^{T_{\mathrm{out}}} \lambda _{1}(\mu (\varepsilon s)) \, \mathrm{d} s = 0, \qquad \text{ or, equivalently, }\quad \int _{\mu _{\mathrm{in}}}^{\mu _{\mathrm{out}}} \lambda _{1}(\mu ) \, \mathrm{d}\mu = 0, \end{equation}

holds. Note that the latter integral describes precisely the vanishing $\log$ -amplitude $b_1=0$ from Section 3.1.

Lemma 3.1 (Equal-area prediction). The system (3.12) together with the boundary conditions (3.16) has a unique solution for $0\lt \delta, \mu _{1,\mathrm{in}}-\mu _{1,\mathrm{c}} \ll 1$ , sufficiently small, and we have the expansion

(3.19) \begin{equation} \mu _{\mathrm{out}}=\mu _{\mathrm{out}}^0+\mathrm{O}(\varepsilon ), \end{equation}

where $\mu _{\mathrm{out}}^0\lt \mu _{\mathrm{in}}$ is the largest solution to the equal-area condition (3.18). In particular, $\mu _{1,\mathrm{c}}-\mu _{\mathrm{out}}=\mu _{\mathrm{in}}-\mu _{1,\mathrm{c}}+\mathrm{O}(|\delta |+|\varepsilon |+|\mu _{\mathrm{in}}-\mu _{1,\mathrm{c}}|^2)$ .

Proof. The result states that at leading order, the linear prediction (3.18) is correct. For small $\delta$ , the eigenvalue $\lambda _{1}(\mu )$ can be approximated near $\mu _{\mathrm{1,c}}$ by $\lambda _{1}(\mu )=\lambda _{1}'(\mu _{1,\mathrm{c}})(\mu -\mu _{1,\mathrm{c}})+\mathrm{O}(2)$ , which implies the last statement.

In order to show that nonlinear terms contribute only at higher order, one follows the analysis of the slow passage near a pitchfork bifurcation in [39]. We omit the somewhat lengthy but straightforward details.

Corollary 3.2 (drop-by-1 result). Assume Hypothesis 2.1 and an initial condition in a sufficiently small neighbourhood to the $j$ -mode equilibrium for some $j\gt 0$ , a parameter-value $\mu _{\mathrm{in}}\gtrsim \mu _{\mathrm{1,c}}$ , and let $\varepsilon$ be sufficiently small. Then the modal number will drop by one at $\mu _{\mathrm{out}}=2\mu _{\mathrm{1,c}}-\mu _{\mathrm{in}}+\mathrm{O}(\varepsilon )$ .

Proof. The result follows immediately from Lemma 3.1 and Hypothesis 2.1.

From Corollary 3.2, we conclude that for $\mu _{\mathrm{in}}$ sufficiently close to $\mu _{1,\mathrm{c}}$ , we will always observe a drop by one and can easily characterise the drop time. Naturally, we wish to see how far this result extends. The key limiting factor in the result is the presence of a leading eigenvalue with $\lambda _{1}(\mu )\gt \lambda (\mu )$ for all other $\lambda (\mu )$ in the spectrum of the linearisation. Technically, centre manifold reductions in the literature also require a uniform splitting, $\lambda _{1}(\mu )\gt \eta _0\gt \lambda (\mu )$ with $\eta _0$ independent of $\mu$ , possibly with sufficiently large gaps to ensure smoothness. We believe that a simple splitting of the leading eigenvalue $\lambda _{1}(\mu )\gt \lambda (\mu )$ for all $\mu \in (\mu _{\mathrm{in}},\mu _{\mathrm{out}})$ should be sufficient to establish Lemma 3.1 but will not attempt to do so here. We observed excellent validity of our predictions even when the uniform gap condition does not hold.

The procedure thus far predicts how and when a first transition will lead to a drop by one in wavenumber or possibly a more severe drop. At the end of this drop, from say $E_j$ to $E_{j-1}$ near $\mu =\mu _{\mathrm{out}}$ , we have followed a global heteroclinic to a neighbourhood of $E_{j-1}$ , where we typically expect that upon entering the neighbourhood, we have a generic perturbation of this stable equilibrium and can now repeat the analysis with the new $\mu _{\mathrm{in}}$ given by the previous $\mu _{\mathrm{out}}$ . A short calculation using our main prediction (3.10) and expressions for eigenvalues in (2.5) shows that for $j\gg 1$ , the distance between $\mu _{1,\mathrm{c}}$ for modes $j$ and $j+1$ is $6j+1$ , while the distance between $\mu _{1,\mathrm{c}}$ and $\mu _{2,\mathrm{c}}$ is only $\frac{3}{2}$ . Hence, one expects that, for large $j$ , a drop-by-1 is followed by a large distance to criticality in parameter space and a potentially subsequent higher drop; compare also Figure 3. Therefore, we now analyse such higher drops, when more than one eigenvalue is unstable at the predicted drop time.

3.3 The drop-by-2 transition

We first consider the static problem with $\mu _{3,\mathrm{c}} \lt \mu \lt \mu _{\mathrm{2,c}}$ , so that the equilibrium $E_j$ has two unstable eigenvalues $\lambda _{1}(\mu )$ and $\lambda _{2}(\mu )$ . There is an associated complex two-dimensional, real four-dimensional centre-unstable manifold. The dynamics on this manifold are governed by equations for the complex amplitudes $a_1, a_2$ in the associated eigenspaces,

(3.20) \begin{align} \dot{a}_1 &= \lambda _{1}(\mu ) a_1 + c_{1,2} (\mu ) \bar{a}_1 a_2 +\mathrm{O}(|a_1|(|a_1|^2+|a_2|^2)),\nonumber \\[5pt] \dot{a}_2 &= \lambda _{2} (\mu ) a_2 + c_{1^2} (\mu ) a_1^2 +\mathrm{O}(|a_1|(|a_1|^2+|a_2|^2)), \end{align}

with coefficients $c_{1,2}(\mu ), c_{1^2}(\mu ) \in \mathbb{C}$ that can be readily computed by evaluating the nonlinearity on the associated eigenspace and projecting. Again, the vector field commutes with the isotropy of the equilibrium, that is, translations $(a_1,a_2)\mapsto (\mathrm{e}^{\mathrm{i}\varphi } a_1,\mathrm{e}^{2\mathrm{i}\varphi } a_2),$ and reflections $(a_1,a_2)\mapsto (\bar{a}_1,\bar{a}_2)$ . As consequence, $c_{1^2}$ and $c_{1,2}$ are real and the real subspace $(a_1,a_2)\in \mathbb{R}^2$ is invariant. Note that the nonlinear terms correspond to the spatial resonances referred to in (3.4): the $a_j$ are amplitudes of Fourier modes $\mathrm{e}^{\mathrm{i} j x}$ , and the complex rotation is induced by the spatial shift $x\mapsto x+\varphi$ in the full equation. The associated manifold is normally hyperbolic in parameter regimes where $\lambda _{1}(\mu )$ and $\lambda _{2}(\mu )$ are separated from other eigenvalues; its smoothness (making in particular the cubic terms meaningful) is determined by spectral gaps; see for instance [28, 30]. These gaps can be easily evaluated in the large- $j$ limit (2.5) where $|\lambda _{2}(\mu )|/|\lambda _{3}(\mu )|\gt 2$ for all $\hat{\mu }\gt 0$ , providing the desired smoothness of the manifold and validity of expansions.

For $\varepsilon \neq 0$ , this manifold persists as a slow manifold with dynamics at quadratic order given by

(3.21) \begin{align} \dot{a}_1 &= \lambda _{1}(\mu ) a_1 + c_{1,2} (\mu ) \bar{a}_1 a_2, \end{align}

(3.22) \begin{align} \dot{a}_2 &= \lambda _{2} (\mu ) a_2 + c_{1^2} (\mu ) a_1^2, \end{align}

(3.23) \begin{align} \dot{\mu } &= - \varepsilon. \end{align}

To predict drop times and drop numbers, we will try to determine when $(a_1,a_2)$ leave a small neighbourhood of the origin, and whether $|a_1|\gg |a_2|$ or $|a_2|\gg |a_1|$ at that time, thus concluding that the local drop number is $1$ or $2$ , respectively. Hypothesis 2.1 then allows us to conclude global drop numbers. We therefore consider the system (3.21)–(3.23) with initial conditions $a_1({-}T_{\mathrm{in}}) = a_2 ({-}T_{\mathrm{in}}) = \delta$ , $\mu ({-}T_{\mathrm{in}}) = \mu _{\mathrm{in}}$ , with $\mu _{\mathrm{in}} \gt \mu _{1,\mathrm{c}} \gt \mu _{2,\mathrm{c}}$ so that both eigenvalues are initially stable. As we shall quickly see, the initial value of $a_2$ is irrelevant so that this particular choice is not restrictive. The solution will therefore initially contract, and we look for the first time $T_{\mathrm{out}}$ at which

(3.24) \begin{align} \max \left \{ | a_1(T_{\mathrm{out}} )|, |a_2 (T_{\mathrm{out}} )| \right \} = \delta. \end{align}

Asymptotics from variation-of-constant formula. We neglect higher order terms and study (3.21)–(3.23) setting $c_{1,2}\equiv 0$ , noticing that $a_2\ll 1$ so that $|c_{1,2}\bar{a_1}a_2|\ll |a_1|$ ; see also (3.35) and the discussion there. We also use $\mu$ as an equivalent time variable and find amplitudes

(3.25) \begin{align} a_1(\mu )&=\mathrm{e}^{\frac{1}{\varepsilon }\int _{\mu }^{\mu _{\mathrm{in}}}\lambda _1(\nu )\mathrm{d}\nu }\delta,\nonumber \\[5pt] a_2(\mu )&=a_2^0(\mu )+a_2^{1^2}(\mu ),\nonumber \\[5pt] a_2^0(\mu )&=\mathrm{e}^{\frac{1}{\varepsilon }\int _{\mu }^{\mu _{\mathrm{in}}}\lambda _2(\nu )\mathrm{d}\nu }\delta,\nonumber \\[5pt] a_2^{1^2}(\mu )&=\int _{\mu }^{\mu _{\mathrm{in}}}c_{1^2}(\nu )\mathrm{e}^{\frac{1}{\varepsilon }\Lambda (\nu )}\mathrm{d}\nu \,\delta ^2,\nonumber \\[5pt] \Lambda (\nu ;\,\mu )&=\int _\mu ^\nu \lambda _2(\rho )\mathrm{d}\rho + \int _\nu ^{\mu _{\mathrm{in}}}2\lambda _1(\rho )\mathrm{d}\rho. \end{align}

Here, we suppress the dependence on $\mu _{\mathrm{in}}$ in $\Lambda (\nu ;\,\mu )$ and notice that $\Lambda (\nu ;\,\mu )$ is smooth and maximal when $\nu =\mu _{1^2:2}$ so that we can expand the integrand near $\nu =\mu _{1^2:2}$ to find

(3.26) \begin{align} a_2^{1^2}(\mu )&=c_{1^2}\mathrm{e}^{\frac{1}{\varepsilon }\Lambda }\delta ^2 \int _{\mu }^{\mu _{\mathrm{in}}}\mathrm{e}^{\frac{1}{2\varepsilon } \Lambda _{\nu \nu }(\nu -\mu _{1^2:2})^2}\mathrm{d}\nu (1+\mathrm{O}(\sqrt \varepsilon ))\nonumber \\[5pt] &=c_{1^2}\mathrm{e}^{\frac{1}{\varepsilon }\Lambda }\delta ^2\sqrt{\frac{-2\varepsilon }{\Lambda _{\nu \nu }}}\mathrm{err}\,\left (\frac{\mu _{1^2:2}-\mu }{\sqrt{\varepsilon }}\right )\left (1+\mathrm{O}(\sqrt \varepsilon )\right ), \end{align}

where $\Lambda _{\nu \nu }=(\lambda _2-2\lambda _1)'\lt 0$ is the derivative in $\mu$ of the resonance condition, and $c_{1^2},\Lambda$ , and $\Lambda _{\nu \nu }$ are evaluated at $\nu =\mu _{1^2:2}$ , and $\mathrm{err}\,(z)=\int _{-\infty }^z \exp ({-}z^2)\mathrm{d} z$ .

First, one sees that $a_2^0\ll a_2^{1^2}$ . In order to identify the largest amplitude upon exiting a neighbourhood, thus determining if the local drop number is 1 or 2, we therefore need to compare $a_1$ and $a_2\sim a_2^{1^2}$ . The transition from a local drop number of 1 to a drop number 2 occurs when

(3.27) \begin{equation} a_2^{1^2}=\delta,\quad \text{and }\quad a_1=\delta. \end{equation}

Solving these two equations for $\mu _{\mathrm{in}}$ and $\mu_\mathrm{out}$ then identifies the transition parameter values $\mu _{\mathrm{in},1\to 2}$ and $\mu _{\mathrm{out},1\to 2}$ . In order to solve (3.27), we use the expression for $a_1$ to find $\int _{\mu _{\mathrm{out},1\to 2}}^{\mu _{\mathrm{in}}}\lambda _1(\nu )\mathrm{d}\nu =0$ . Substituting this result into the expression for $\Lambda$ in (3.25), we obtain

(3.28) \begin{equation} \Lambda (\mu _{1^2:2},\mu _{\mathrm{out},1\to 2})=\int _{\mu _{\mathrm{out},1\to 2}}^{\mu _{\mathrm{1^2:2}}}(\lambda _2(\rho )-2\lambda _1(\rho ))\mathrm{d}\rho. \end{equation}

We next assume that $\Delta \mu =\mu _{1^2:2}-\mu _{\mathrm{out},1\to 2}\gt 0$ is small so that we can Taylor expand

\begin{equation*} \Lambda (\mu _{1^2:2},\mu _{\mathrm {out},1\to 2})=-\frac {1}{2}\Lambda _{\nu \nu }(\mu _{1^2:2},\mu _{1^2:2})(\Delta \mu )^2+\mathrm {O}((\Delta \mu )^3).\end{equation*}

Inserting this into (3.26) and partially solving for $\Delta \mu$ , we find

(3.29) \begin{equation} \Delta \mu =\sqrt{\frac{2\varepsilon }{\Lambda _{\nu \nu }}\log \left ( c_{1^2}\delta \sqrt{-2\varepsilon/(\Lambda _{\nu \nu })} \mathrm{err}\,(\Delta \mu/\sqrt{\varepsilon })\right )}. \end{equation}

Assuming that at leading order $\Delta \mu \sim \sqrt{-\varepsilon \log \varepsilon }$ , we find $\mathrm{err}\,(\Delta \mu/\sqrt{\varepsilon })\to \sqrt{\pi }$ , so that, consistently, at leading order

(3.30) \begin{equation} \Delta \mu =\sqrt{\frac{2\varepsilon }{\Lambda _{\nu \nu }}\log \left ( c_{1^2}\delta \sqrt{-2 \pi \varepsilon/(\Lambda _{\nu \nu })}\right )}. \end{equation}

Away from the transition, the amplitude of $a_2(\mu )$ is given at leading order through

(3.31) \begin{equation} a_2(\mu )=\left \{\begin{array}{ll} \mathrm{O}(\varepsilon a_1^2(\mu )), & \mu \gt \mu _{1^2:2},\\[5pt] \exp \left ({\frac{1}{\varepsilon } \int _{\mu }^{\mu _{1^2:2}}\lambda _2(\rho )\mathrm{d}\rho +\frac{1}{\varepsilon } \int _{\mu _{1^2:2}}^{\mu _{\mathrm{in}}}2\lambda _1(\rho )\mathrm{d}\rho }\right ),& \mu \lt \mu _{1^2:2}, \end{array} \right. \end{equation}

or, neglecting $\log \varepsilon$ -terms,

(3.32) \begin{equation} b_2(\mu )=\varepsilon \log (a_2(\mu ))=\int _\mu ^{\mu _{\mathrm{in}}} \Lambda _{\mathrm{max}}(\rho )\mathrm{d}\rho, \qquad \Lambda _{\mathrm{max}}(\rho )=\mathrm{max}\,\{\lambda _2(\rho ),2\lambda _1(\rho )\}=\left \{\begin{array}{ll}2\lambda _1(\rho ),& \rho \gt \mu _{1^2:2},\\[5pt] \lambda _2(\rho ),& \rho \lt \mu _{1^2:2}. \end{array}\right. \end{equation}

For $\mu _{\mathrm{out}}\lt \mu _{\mathrm{out},1\to 2}$ , at leading order the local drop parameter value is given by setting $a_2^{1^2}(\mu )=\delta$ in (3.25), which gives

(3.33) \begin{equation} \int _{\mu _{\mathrm{out}}}^{\mu _{1^2:2}}\lambda _2(\rho )\mathrm{d}\rho + \int _{\mu _{1^2:2}}^{\mu _{\mathrm{in}}}2\lambda _1(\rho )\mathrm{d}\rho \stackrel{!}{=}0, \end{equation}

in agreement with our predictions in Section 3.1 and with direct simulations illustrated in Figure 8.

Figure 8. Numerical comparisons for the drop-by-1 to drop-by-2 transition. Top: local and global drop numbers depending on $\mu _{\mathrm{in}}$ (left) and $\mu _{\mathrm{out}}$ (right) for several values of $\varepsilon$ ; predicted drop as solid black curve. Bottom: $\mu _{\mathrm{out}}$ versus $\mu _{\mathrm{in}}$ for several values of $\varepsilon$ together with prediction for $\varepsilon =0$ (solid black) (left); drop values $\mu _{\mathrm{in}}$ both local and global as observed in top row shown here as a function of $\varepsilon$ , with linear fit and predicted values of drops as green dot. Extrapolated values at $\varepsilon =0$ from linear fits fall within 10% of the predicted value for local drops and within 2% for global drops (right). Throughout, $\mu =\mu _{\mathrm{c}}+\hat{\mu }$ indicates values relative to criticality.

Figure 9. Top: dynamics of centre manifold equations (3.20) with instability thresholds for $a_1$ and $a_2$ and $\lambda _{1/2}=0$ , respectively, and subsequent 1:1 and 1:2 resonances; $\mu$ decreases to the left, the vertical plane $a_1=0$ is invariant. Sample dark and bright red trajectories illustrating trajectories exiting in the $a_1$ and $a_2$ directions, respectively. Bottom: dynamics in projective coordinates (3.35) with invariant planes $a_2=0$ and $\xi =a_1^2/a_2=0$ ; trajectories follow the stable branch in a transcritical bifurcation in the horizontal plane with unstable manifold in purple and blue for nontrivial and trivial branch, respectively; trajectories exit in the $a_2$ direction along the unstable manifold, albeit at locations where $a_1^2/a_2\neq 0$ , so $|a_1|\sim \sqrt{|a_2|}\gg |a_2|$ (bright red), or when $a_1^2/a_2\sim 0$ and $|a_2|\gg |a_1|$ (dark red).

Rigorous asymptotics and a geometric view on resonances. We next present a geometric view that also explains why higher order terms are irrelevant. Consider the new variable $\xi =a_1^2/a_2$ , together with $a_2$ , and find the new equations

(3.34) \begin{align} \xi '=&(2\lambda _{1}(\mu )-\lambda _{2}(\mu ))\xi -c_{1^2}(\mu )\xi ^2 +\mathrm{O}(|a_2|)\nonumber,\\[5pt] a_2'=&\lambda _{2}(\mu )a_2+c_{1^2}(\mu )\xi a_2+\mathrm{O}(|a_2|^2), \end{align}

(3.35) \begin{align} \mu '=&\varepsilon. \end{align}

In the invariant plane $a_2=0$ , we find a slow passage through a transcritical bifurcation at the resonance, as $2\lambda _{1}(\mu )-\lambda _{2}(\mu )$ passes through zero. Before the resonance, $2\lambda _{1}(\mu )-\lambda _{2}(\mu )\gt 0$ , and the nontrivial equilibrium $\xi _*=(2\lambda _{1}(\mu )-\lambda _{2}(\mu ))/c_{1^2}$ is exponentially attracting; past the resonance, trajectories follow the now stable trivial equilibrium $\xi =0$ . Normal to this stable family of equilibria, trajectories first decay, then grow in the $a_2$ -direction with the eigenvalue $\lambda _{2}+c_{1^2}\xi _*=2\lambda _{1}$ ; after the passage through the transcritical bifurcation growth is governed by the normal eigenvalue $\lambda _{2}$ , thus reproducing the integral criterion (3.33). A refined desingularisation, analysing for instance the slow passage through the transcritical bifurcation using geometric blowup as in [39] should then lead to a leading-order expansion of the form obtained through the direct calculations outlined above; see Figure 9 for an illustration.

3.4 The drop-by-3 transition

Turning to predictions for higher drops, we need to account for more unstable eigenvalues and therefore higher resonances. We first adapt the asymptotic calculation for the 1:2-resonance criteria above to the drop-by-3 transition.

Once $\lambda _{3}$ becomes positive, we track the six-dimensional unstable manifold associated to the first three complex unstable modes in the static problem. Keeping track of what interactions of modes $\mathrm{e}^{\pm 3ix}, \mathrm{e}^{\pm 2 ix}, \mathrm{e}^{\pm ix}$ can produce the original modes $\mathrm{e}^{ix}, \mathrm{e}^{2ix}, \mathrm{e}^{3ix},$ we see that to leading order, the dynamics on this unstable manifold are given by

(3.36) \begin{align} \dot{a}_1 &= \lambda _{1}(\mu ) a_1 + c_{-1, 2} \bar{a}_1 a_2 + c_{3,-2} a_3 \bar{a}_2, \end{align}

(3.37) \begin{align} \dot{a}_2 &= \lambda _{2}(\mu ) a_2 + c_{1^2} a_1^2 + c_{3, -1} a_3 \bar{a}_1, \end{align}

(3.38) \begin{align} \dot{a}_3 &= \lambda _{3} (\mu ) a_3 + c_{1^3} a_1^3 + c_{1,2} a_1 a_2, \end{align}

(3.39) \begin{align} \dot{\mu } &= - \varepsilon, \end{align}

also taking into account the slow evolution of $\mu$ . Error terms can be seen to be small either heuristically or in the analogue of the geometric desingularisation shown in Figure 9; see the discussion below. We again consider this system with initial conditions $a_1 ({-}{T_{\mathrm{in}}}) = a_2 ({-}{T_{\mathrm{in}}}) = a_3 ({-}{T_{\mathrm{in}}}) = \delta$ , $\mu ({-}T_{\mathrm{in}})=\mu _{\mathrm{in}}$ , and look for the first time $T_{\mathrm{out}}$ for which $|a_j({T_{\mathrm{out}}})| = \delta$ for some $j = 1,2$ or $3$ .

The coefficients of the nonlinear terms depend on $\mu$ and hence are also slowly varying, but since the values of these coefficients will not affect our predictions at leading order, we suppress this dependence. Also note that we have included the cubic term $c_{1^3}a_1^3$ in the equation for $a_3$ , since the analysis of Section 3.3 suggests that for some time we have $a_2 \sim a_1^2$ , and so $a_1^3$ terms are in fact on the same order as $a_1 a_2$ for some relevant time. Following a similar line of reasoning, we see as in Section 3.3 that the back-coupling terms involving $a_2$ and $a_3$ in the equation for $a_1$ , and $a_3$ in the equation for $a_2$ , are higher order. We, therefore, neglect these terms and arrive at the system

(3.40) \begin{align} \dot{a}_1 &= \lambda _{1}(\mu ) a_1, \end{align}

(3.41) \begin{align} \dot{a}_2 &= \lambda _{2}(\mu ) a_2 + c_{1^2} a_1^2, \end{align}

(3.42) \begin{align} \dot{a}_3 &= \lambda _{3} (\mu ) a_3 + c_{1^3} a_1^3 + c_{1,2} a_1 a_2, \end{align}

(3.43) \begin{align} \dot{\mu } &= - \varepsilon. \end{align}

We can now solve this system recursively, first solving (3.40) for $a_1$ , then solving for $a_2$ via its variation of constants formula, and then inserting both expressions for $a_1$ and $a_2$ into (3.42) and solving via the variation of constants formula. We shall use the expressions for $a_1$ and $a_2$ derived above and are left with the equation for $a_3$ , which we write as a sum of the solution to the homogeneous equation $a_3^0$ , the solution $a_3^{1^3}$ with inhomogeneity $c_{1^3}a_1^3$ , and the solution $a_3^{1,2}$ with inhomogeneity $c_{1,2}a_1a_2$ . We find

(3.44) \begin{align} a_3^0(\mu )&=\mathrm{e}^{\frac{1}{\varepsilon }\int _{\mu }^{\mu _{\mathrm{in}}}\lambda _3(\nu )\mathrm{d}\nu }\delta,\nonumber \\[5pt] a_3^{1^3}(\mu )&=\int _{\mu }^{\mu _{\mathrm{in}}}c_{1^3}(\mu )\mathrm{e}^{\frac{1}{\varepsilon }\Lambda _{{1^3}}(\nu )}\mathrm{d}\nu \delta ^3,\nonumber \\[5pt] \Lambda _{1^3}(\nu )&=\int _\mu ^\nu \lambda _3(\rho )\mathrm{d}\rho + \int _\nu ^{\mu _{\mathrm{in}}}3\lambda _1(\rho )\nonumber \\[5pt] a_3^{1,2}(\mu )&=\int _{\mu }^{\mu _{\mathrm{in}}}c_{1,2}(\mu )\mathrm{e}^{\frac{1}{\varepsilon }\Lambda _{{1,2}}(\nu )}\mathrm{d}\nu \delta ^3,\nonumber \\[5pt] \Lambda _{1,2}(\nu )&= \int _\mu ^\nu \lambda _3(\rho )\mathrm{d}\rho + \int _\nu ^{\mu _{1^2,2}}(\lambda _1(\rho )+\lambda _2(\rho ))\mathrm{d}\rho + \int _{\mu _{1^2,2}}^{\mu _{\mathrm{in}}}3\lambda _1(\rho )\mathrm{d}\rho. \end{align}

Following the same strategy as in the analysis for near the $1^2$ :2-resonance, we evaluate the integrals to leading order near the maximum of the exponential, finding

(3.45) \begin{align} b_3^{1^3}&=\varepsilon \log |a_3^{1^3}|=\Lambda _{1^3}(\mu _{1^3:3})+\mathrm{O}(\varepsilon ^{1^-}), \end{align}

(3.46) \begin{align} b_3^{1,2}&=\varepsilon \log |a_3^{1,2}|=\Lambda _{1,2}(\mu _{1,2:3})+\mathrm{O}(\varepsilon ^{1^-}), \end{align}

where we used the short hand $1^-$ to denote any exponent less than 1, accounting for potential terms of the form $\varepsilon \log \varepsilon$ . Clearly, $b_3^{1,2} \gt b_3^{1^3} \gt b_3^0$ , so that the amplitude of $a_3$ is at leading order given through $\exp (b_3^{1,2}/\varepsilon )$ . The crossover occurs when $a_3\sim a_2$ , that is, at leading order, when

(3.47) \begin{align} &\int _{\mu _{\mathrm{out},2\to 3}}^{\mu _{1,2:3}} \lambda _3(\rho )\mathrm{d}\rho + \int _{\mu _{1,2:3}}^{\mu _{1^2:2}}(\lambda _1(\rho )+\lambda _2(\rho ))\mathrm{d}\rho + \int _{\mu _{1^2:2}}^{\mu _{\mathrm{in},2\to 3}}3\lambda _1(\rho )\mathrm{d}\rho \end{align}

(3.48) \begin{align} \stackrel{!}{=} &\int _{\mu _{\mathrm{out},2\to 3}}^{\mu _{1^2:2}} \lambda _2(\rho )\mathrm{d}\rho + \int _{\mu _{1^2,2}}^{\mu _{\mathrm{in},2\to 3}}2\lambda _1(\rho )\mathrm{d}\rho + \mathrm{O}(\varepsilon ^{1^-}) \end{align}

(3.49) \begin{align} \stackrel{!}{=} &\, \mathrm{O}(\varepsilon ^{1^-}). \end{align}

The two equalities determine $\mu _{\mathrm{out},\,2\to 3}$ and $\mu _{\mathrm{in},\,2\to 3}$ . Linearising at a solution with respect to these two variables, we find a Jacobi matrix with determinant

\begin{equation*} \left |\begin {array}{cc} -\lambda _3(\mu _{\mathrm {out}})+\lambda _2(\mu _{\mathrm {out}})& 3\lambda _1(\mu _{\mathrm {in}})\\[5pt] -\lambda _2(\mu _{\mathrm {out}})& 2\lambda _1(\mu _{\mathrm {in}}) \end {array}\right |=\lambda _1(\mu _{\mathrm {in}})\left ({-}2 \lambda _3(\mu _{\mathrm {out}})+ 5 \lambda _2(\mu _{\mathrm {out}})\right ). \end{equation*}

Assuming the additional non-resonance condition $2\lambda _3(\mu _{\mathrm{out}})\neq 5 \lambda _2(\mu _{\mathrm{out}})$ , we therefore expect to be able to solve for the variables $\mu _{\mathrm{in}}$ and $\mu _{\mathrm{out}}$ with the implicit function theorem with corrections to the transition value as $\mathrm{O}(\varepsilon ^{1^-})$ , rather than $\mathrm{O}(\sqrt{\varepsilon })$ for the 1-2-transition. The results agree with our summary in Section 3.1 and with numerical simulations shown in Figure 10.

The geometric picture for the drop-by-3 transition. We can follow the strategy for the analysis of the drop-by-2 transition and introduce projective variables that encode the quotients of amplitudes and resonances. For instance, set $ \xi _2=a_1^2/a_2, \xi _3=a_1a_2/a_3,$ writing for short $\lambda _j=\lambda _{j}(\mu )$ , suppressing $\mu$ -dependence in the nonlinear terms, and suppressing higher-order terms in $a_1$ , to find

\begin{align*} a_1'&=\lambda _1 a_1,\\[5pt] \xi _2'&=(2\lambda _1-\lambda _2)\xi _2 -c_{1^2:2} \xi _2^2,\\[5pt] \xi _3'&=(\lambda _1+\lambda _2-\lambda _3)\xi _3 - c_{1,2:3}\xi _3^2 - c_{1^3:3}\xi _3^2\xi _2. \end{align*}

We previously analysed the dynamics in the $(a_2,\xi _2)$ -subsystem. In the regime $\mu \gt \mu _{1^2:2}$ , $\xi _3$ follows the nontrivial stable branch with $\xi _3=(\lambda _1+\lambda _2-\lambda _3)/(c_{1,2:3} - c_{1^3:3}\xi _2)$ . For $\mu _{1,2:3}\lt \mu \lt \mu _{1^2:2}$ , $\xi _2\sim 0$ and $\xi _3=(\lambda _1+\lambda _2-\lambda _3)/c_{1,2:3}$ become the nontrivial stable branch, which at $\mu _{1,2:3}$ exchanges stability with the trivial branch $\xi _3=0$ in a transcritical bifurcation. In particular, $\xi _3\sim 0$ for $\mu \lt \mu _{1,2:3}$ . To reproduce the explicit computations above, we recover growth by analysing exponential growth in the normal direction of $a_1$ , and of $a_{2/3}$ through the values of $\xi _{2,3}$ .

Figure 10. Numerical comparisons for the drop-by-2 to drop-by-3 transition. Top: local and global drop numbers depending on $\mu _{\mathrm{in}}$ (left) and $\mu _{\mathrm{out}}$ (right) for several values of $\varepsilon$ ; predicted drop as solid black curve. Bottom: $\mu _{\mathrm{out}}$ versus $\mu _{\mathrm{in}}$ for several values of $\varepsilon$ together with prediction for $\varepsilon =0$ (solid black) (left); drop values $\mu _{\mathrm{in}}$ both local and global as observed in top row shown here as a function of $\varepsilon$ , with linear fit and predicted values of drops as green dot. Extrapolated values at $\varepsilon =0$ from linear fits fall within 1% of the predicted value (right). Throughout, $\mu =\mu _{\mathrm{c}}+\hat{\mu }$ indicates values relative to criticality.

3.5 Drop-by-4 and beyond

Predicting higher drops appears to be cumbersome in general. We outline here the calculations for the transition from a drop-by-3 to a drop-by-4. The relevant amplitude equations are

(3.50) \begin{align} \dot{a}_1 &= \lambda _{1}(\mu ) a_1, \end{align}

(3.51) \begin{align} \dot{a}_2 &= \lambda _{2}(\mu ) a_2 + c_{1^2} a_1^2, \end{align}

(3.52) \begin{align} \dot{a}_3 &= \lambda _{3} (\mu ) a_3 + c_{1,2} a_1 a_2+ c_{1^3} a_1^3, \end{align}

(3.53) \begin{align} \dot{a}_4 &= \lambda _{4} (\mu ) a_4 + c_{1,3} a_1 a_3+ c_{2,2} a_2^2 + +c_{1^2,2}a_1^2a_2+ c_{1^4} a_1^4, \end{align}

(3.54) \begin{align} \dot{\mu } &= - \varepsilon. \end{align}

Higher modes will be irrelevant until growth dominates the source terms from resonant interactions. In order to determine the onset of growth in the mode $a_4$ , one would inspect the lowest order source terms, $a_2^2$ and $a_1a_3$ , and find the associated resonances, $\mu _{2^2:4}$ and $\mu _{1,3:4}$ , which give a lower bound for the crossover to a drop-by-4.

It is convenient to track approximations to the logarithms of the amplitudes, $b_k=\varepsilon \log a_k$ , so that for instance

(3.55) \begin{align} b_1(\mu )&=\int _\mu ^{\mu _{\mathrm{in}}} \lambda _1(\mu )\mathrm{d}\mu,\nonumber \\[5pt] b_2(\mu )&=\int _\mu ^{\mu _{1^2:2}} \lambda _2(\mu )\mathrm{d}\mu +2b_1(\mu _{1^2:2}),\nonumber \\[5pt] b_3(\mu )&=\int _\mu ^{\mu _{1,2:3}} \lambda _3(\mu )\mathrm{d}\mu +b_1(\mu _{1,2:3}) +b_2(\mu _{1,2:3}). \end{align}

For $b_4$ , we need to consider all resonant terms in the equation for $a_4$ . At leading order, we find that $b_4$ is given as a maximum of the expressions obtained by treating all terms separately,

\begin{align*} b_4^{1,3}(\mu )&=\int _\mu ^{\mu _{1,3:4}} \lambda _4(\mu )\mathrm{d}\mu +b_1(\mu _{1,3:4})+b_3(\mu _{1,3:4}),\\[5pt] b_4^{2,2}(\mu )&=\int _\mu ^{\mu _{1,3:4}} \lambda _4(\mu )\mathrm{d}\mu +2b_2(\mu _{1,3:4}),\\[5pt] \ldots & \end{align*}

and similar expressions for cubic and quartic resonant terms. At any given $\mu \lt \mu _{1,3:4}$ , the amplitude is given by the maximum $b_4=\max _m b_j^m$ where $m$ runs through all resonant source terms.

It turns out that $b_4^{2,2}$ maximises the amplitudes, and we can find the value of the transition from a drop-by-3 to a drop-by-4 by solving $b^4_{2,2}(\mu _{\mathrm{out}})=0$ together with $b^3(\mu _{\mathrm{out}})=0$ for $\mu _{\mathrm{out}}$ and $\mu _{\mathrm{in}}$ . For finite $j$ , this equation can easily be solved numerically, and the predicted drop transitions compare well with numerical simulations as shown in Figure 11.

Generalising to higher transitions appears straightforward in principle but cumbersome. One obtains successively $\log$ -amplitudes $b_1,b_2,b_3,b_4,\ldots,b_\ell,\ldots$ for a given $\mu _{\mathrm{in}}$ as functions of $\mu$ , maximising at each $\ell$ and each $\mu$ over all potential resonant source terms. Solving for $b_\ell (\mu )=0$ as a function of $\mu$ given $\mu _{\mathrm{in}}$ then yields the local drop parameter value $\mu _{\mathrm{out}}$ . The argmax of all the exit values $\mu _{\mathrm{out}}$ then determines the local drop number for a given $\mu _{\mathrm{in}}$ . We caution, however, that at some point the skew product structure may not be preserved since, for instance, terms of the form $a_3\bar{a_2}$ in the equation for $a_1'$ may dominate the linear term $\lambda _1a_1$ .

Figure 11. Drop-to numbers for initial patterns with $j=3,\ldots,8$ varying $\mu _{\mathrm{in}}$ , plotted against $\mu _{\mathrm{in}}-\mu _{\mathrm{c}}$ (top) and $\mu _{\mathrm{out}}-\mu _{\mathrm{c}}$ (bottom). Shown are both local (blue) and global drops (yellow), which always yield an equal drop or a drop by one less than the local drop. The actual drop-to number is the integer part of the plotted value (which is slightly shifted to improve readability). Red markers indicate the theoretical prediction for drops.

3.6 Drop number transitions for large $\text{j}$

In the limit of large $j$ , the criteria for transitions can be evaluated explicitly. Using the expansion for eigenvalues (2.5), we can find the resonances as defined in (3.3) explicitly. Setting $\mu =3j^2-\frac{1}{2}+\hat{\mu }$ , we have in decreasing order,

\begin{equation*} \hat {\mu }_{1^2:2}=-3, \quad \hat {\mu }_{1^3:3}=-6, \quad \hat {\mu }_{1,2:3}=-\frac {15}{2}, \quad \hat {\mu }_{1,3:4}=-14, \quad \hat {\mu }_{1,4:5}=-\frac {45}{2}, \end{equation*}

and

\begin{equation*} \hat {\mu }_{1^4:4}=-10, \quad \hat {\mu }_{1^2,2:4}=-\frac {114}{10}, \quad \hat {\mu }_{2^2:4}=-\frac {27}{2}, \quad \hat {\mu }_{1,3:4}=-14, \ldots \end{equation*}

Evaluating the log-amplitudes in (3.5)–(3.8) explicitly for eigenvalues as in (2.5), we also find explicit crossover points $\hat{\mu }_{\mathrm{in},1\to 2}$ , $\hat{\mu }_{\mathrm{out},1\to 2}$ , for drop-by-1 to drop-by-2 and $\hat{\mu }_{\mathrm{in},2\to 3}$ , $\hat{\mu }_{\mathrm{out},2\to 3}$ for drop-by-2 to drop-by-3 as

(3.56) \begin{equation} \hat{\mu }_{\mathrm{in},1\to 2}=3,\quad \hat{\mu }_{\mathrm{out},1\to 2}=-3, \qquad \qquad \hat{\mu }_{\mathrm{in},2\to 3}=15,\quad \hat{\mu }_{\mathrm{out},2\to 3}=-12. \end{equation}

Note that the change to a drop-by-3 occurs well before the amplitude $a_4$ becomes relevant.

The drop values $ \hat{\mu }_{\mathrm{out},\ell }$ for a drop by $\ell$ as a function of $\hat{\mu }_{\mathrm{in}}$ are given by

(3.57) \begin{align} \hat{\mu }_{\mathrm{out},1}&=-\hat{\mu }_{\mathrm{in}}, &\qquad 0\lt \hat{\mu }_{\mathrm{in}}\lt \hat{\mu }_{\mathrm{in},1\to 2}, \end{align}

(3.58) \begin{align} \hat{\mu }_{\mathrm{out},2}&=\frac{1}{2} \left ({-}3 - \sqrt{-9 + 2 \hat{\mu }_{\mathrm{in}}^2}\right ), &\qquad \hat{\mu }_{\mathrm{in},1\to 2}\lt \hat{\mu }_{\mathrm{in}}\lt \hat{\mu }_{\mathrm{in},2\to 3}, \end{align}

(3.59) \begin{align} \hat{\mu }_{\mathrm{out},3}&=\frac{1}{6} \left ({-}24 - 2\sqrt{3} \sqrt{-33+ \hat{\mu }_{\mathrm{in}}^2}\right ), &\qquad \hat{\mu }_{\mathrm{in},2\to 3 }\lt \hat{\mu }_{\mathrm{in}}\lt \hat{\mu }_{\mathrm{in},3\to 4 }. \end{align}

For the drop-by-3 to drop-by-4 transition, we find, based on the source term $a_2^2$ and the associated 2,2:4-resonance,

(3.60) \begin{align} \hat{\mu }_{\mathrm{in},3\to 4}&=3 \sqrt{\frac{1}{2} \left (159+28 \sqrt{14}\right )}\nonumber \\[5pt] &=34.4521\ldots,\nonumber \\[5pt] \hat{\mu }_{\mathrm{out},3\to 4}&=\frac{1}{4} \left ({-}30-\sqrt{2 \left (9 \left (159+28 \sqrt{14}\right )-297\right )}\right )\nonumber \\[5pt] &=-23.6125\ldots. \end{align}

The transition associated with the resonant term $a_1a_3$ is located at

\begin{align*} \hat{\mu }^{13}_{\mathrm{in},3\to 4}&=\frac{1}{2} \sqrt{21 \left (147+8 \sqrt{201}\right )}\\[5pt] &=36.9757\ldots,\\[5pt] \hat{\mu }^{13}_{\mathrm{out},3\to 4}&=\frac{1}{4} \left ({-}30-\sqrt{21 \left (147+8 \sqrt{201}\right )-519}\right )\\[5pt] &=-25.0887\ldots, \end{align*}

a larger value of $\mu _{\mathrm{in}}$ . The drop values for finite $j$ have significant corrections, in particular for larger drop numbers or when the drop parameter value is close to the existence boundary.

4. Discussion

We analysed the slow passage through an Eckhaus bifurcation in a bounded domain. Under the conceptual assumption, well corroborated in simulations, that for perturbations of unstable patterns, the dominant Fourier mode of the perturbation determines the global drop, we reduced the analysis to the dynamics in a finite-dimensional centre-manifold. We found that the drop number then depends on how far the initial parameter value is from criticality. For small distances, one observes a drop-by-1, but for larger distances, the drop number increases. We derived criteria for the drop time and the drop number in general for drops up to 3 and found explicit formulas in the case of a large domain (large $j$ in our scaling). The formulas do not show an obvious pattern that would generalise to larger drop numbers, although our approach can in principle be pursued beyond drops by 3 and the transition to drops by 4. The main complicating factor in the analysis is the presence of spatio-temporal resonances, which lead to nonlinear coupling between amplitudes of modes for different drops. Spatial resonances allow for the presence of those nonlinear terms; temporal resonances determine cross-over points when growth in higher modes exceeds the resonant feeding from lower modes.

Some avenues that would be interesting to study further are a generalisation to higher drop numbers, a rigorous justification of the centre manifold analysis when spectral gaps fail, and the identification of relevant higher order terms for moderate values of $\varepsilon$ . On the other hand, it seems natural to pursue this analysis in an unbounded domain, where the Eckhaus instability is caused by essential spectrum. In this setting, the expansion of eigenvalues in the large- $j$ limit that we used throughout is universal beyond the Ginzburg–Landau equation and reflects the sideband nature with modulation equations of the form

(4.1) \begin{equation} a_t=-(a_{xx}+\hat{\mu }a - a^2)_{xx}, \end{equation}

for wavenumber corrections $a$ ; see for instance [58]. From this perspective, the precise form of parameter variation is irrelevant, and paths of the form $(\mu (\tau ),j(\tau ))$ would lead to equivalent results. The quadratic resonances, key to the calculation of drop times, are induced by the nonlinearity $(a^2)_{xx}$ . It would be interesting to relate the universal quadratic term here to the quadratic coefficients modelling resonances. In this direction, it appears tempting to analyse more generally potential resonance orderings and the ensuing sequence of transitions to higher drop numbers, as outlined in Section 3.5. Conceivably, this could be explicit in the limit as $j\to \infty$ and potentially conclude with a complete description of possible staircases as shown in Figure 3. On the other hand, the perspective of a wavenumber modulation equation (4.1) rather than an amplitude modulation such as Ginzburg–Landau demonstrates that local drop numbers should be largely independent of the model considered and predictions, suitably adjusted by computing the relevant eigenvalues $\lambda _j(\tau )$ along a parameter path would still hold. Exploring this systematically for instance in a simple Swift–Hohenberg equation

(4.2) \begin{equation} u_t=-(k(\tau )\partial _{xx}+1)^2 u + \mu (\tau ) u - u^3, \end{equation}

in a large periodic domain, with for instance $0\lt \mu (\tau )\equiv \bar{\mu }\ll 1$ and $k=k_0-\varepsilon \tau$ , $k_0\sim 1$ , or even more realistic reaction-diffusion models such as Klausmeier’s or the Gray–Scott model [35, 52] would be a natural next step.

In this setting of an unbounded domain, rather than describing the global evolution by a heteroclinic orbit, one could start to analyse the spatio-temporal evolution of the Eckhaus instability with frozen parameter in terms of spreading fronts, and, in a second stage, incorporate the effect of spatio-temporally slowly varying parameters; see [24] for the effect of slowly varying environments on spatial spreading, [26] for a review of the effect of spatiotemporal changes on pattern evolution, and [19] for the role of spatiotemporal resonances in spatial spreading.

In a different direction, one would also be interested in instabilities different from the Eckhaus side-band instability, in particular spatial period doubling [53], and in the effect of various types of spatio-temporal noise and variation. Even adding boundary conditions different from the simple periodic setting considered here may well destroy the simple pitchfork nature of the Eckhaus instability; see [42] for a conceptual analysis of the effect of boundaries.

Finally, imperfections in time and space would certainly affect predictions made here. A delay of the bifurcation is still expected but shortened when small temporal noise is present [5], and spatial noise may well change drop predictions in other ways.