Stabilization of uni-directional water-wave trains over an uneven bottom

We study the evolution of nonlinear surface gravity water-wave packets developing from modulational instability over an uneven bottom. A nonlinear Schr\"odinger equation (NLSE) with coefficients varying in space along propagation is used as a reference model. Based on a low-dimensional approximation obtained by considering only three complex harmonic modes, we discuss how to stabilize a one-dimensional pattern in the form of train of large peaks sitting on a background and propagating over a significant distance. Our approach is based on a gradual depth variation, while its conceptual framework is the theory of autoresonance in nonlinear systems and leads to a quasi-frozen state. Three main stages are identified: amplification from small sideband amplitudes, separatrix crossing, and adiabatic conversion to orbits oscillating around an elliptic fixed point. Analytical estimates on the three stages are obtained from the low-dimensional approximation and validated by NLSE simulations. Our result will contribute to understand dynamical stabilization of nonlinear wave packets and the persistence of large undulatory events in hydrodynamics and other nonlinear dispersive media.


I. INTRODUCTION
Modulational instability (MI) is an ubiquitous phenome-non for wave-packets propagating in a weakly nonlinear medium [1]. It consists in the appearance of sidebands growing around a uniformly-modulated carrier and was observed in deep water waves, nonlinear optics, Bose-Ein-stein condensates, and plasma physics [2,3].
If the envelope of the wave-packet is narrowbanded, the nonlinear stage of the evolution (i.e., when the sidebands start to grow at amplitudes comparable to the unstable stationary background) can be modeled by means of the universal nonlinear Schrödinger equation (NLSE). This integrable equation exhibits exact solutions, e.g., stationary envelope solitons and pulsating breathers of Kuznetsov-Ma-, Peregrine-and Akhmediev-type [4][5][6][7] The Akhmediev breather (AB) is the prototype of the nonlinear evolution of MI: in the time-like NLSE, an initially slightly modulated time-periodic train of pulses reaches its peak value at a given point in space, as a result of the exponential sideband growth, as is followed by the recovery of the initial state known as Fermi-Pasta-Ulam recurrence [8]. Because of this characteristic feature, i.e., waves appearing from nothing and suddenly disappearing [9], it is also a candidate solution for the explanation of rogue waves in water and other nonlinear systems.
The NLSE is a framework not only valid for deep water, but also for intermediate depth cases, as is well known from the literature [10][11][12].
The depth is thus an important degree of freedom to tune the dispersion and nonlinear coefficients during wave propagation, thus allowing the possibility to dynamically control the MI gain. In optics, an adiabatic variation of fiber dispersion is well known to provide an effective path to soliton compression [13,14]. Moreover, the transition from two fibers of different dispersion was recently proposed to control an AB at its peak focusing point [15]. A standard fiber has a large cross-section, thus the nonlinear coefficient does not change much (as it depends mostly on the core area and Kerr nonlinear refractive index); the dispersion is instead much more sensitive to geometry [16]. The opposite is true for surface gravity waves in water: the group velocity always decreases with frequency, while the nonlinearity can be tuned to positive or negative values [17].
Here, we propose a theoretical framework for the control of breathing water wave-packets over a smoothly varying uneven bottom. A three-wave truncation [18,19] allows us to formulate the conditions required for stabilization, as well as the limits of our approach.
We rely on a mechanism similar to autoresonance, in which the change of an external parameter in the system allows one to lock it in a stable and stationary oscillating state of large amplitudes, starting from an initial condition close but not exactly matching the resonant condition. This theory finds its origin in accelerator and plasma physics [20][21][22] and was recently applied also to optical frequency conversion [23][24][25][26].
Sec. II recalls the generalized NLSE model of Ref. [27] and the description of the nonlinear stage of MI by means of a three-wave truncation approach. In Sec. III we discuss the conditions for stabilization and report improvements on the implementation of the abrupt transition as proposed in [15]. Numerical results are presented in Sec. IV. Sec. V is devoted to result summary and outlook.

A. Generalized NLSE
In [27], a NLSE-like equation is derived for the onedimensional evolution of the envelope of surface water waves on an uneven bottom of depth h at frequency ω = √ gkσ, with σ ≡ tanh κ and κ ≡ kh, k being the local wavenumber, which varies with h, while ω is fixed.
The 2D Laplace equation ∂ 2 ∂x 2 + ∂ 2 ∂z 2 Φ = 0 for the velocity potential Φ in the longitudinal and depth coordinates (x, z) is solved with the usual kinematic and dynamic boundary conditions at the free surface [12], whereas the bottom boundary condition reads as It is required that the bottom slope is small enough to prevent wave-reflections due to wavenumber mismatches, i.e., h ′ (x) = O(ε 2 ), with ε ≡ ka ≪ 1 the wave steepness (a is the local carrier wave amplitude). By employing the standard method of multiple scales up to third-order in ε [12], the following evolution equation was derived [27] i ∂U ∂ξ where ξ ≡ ε 2 x, τ ≡ ε x 0 dζ cg(ζ) − t are the coordinates in a frame moving at the group velocity of the envelope, c g ≡ ∂ω ∂k = g 2ω σ + κ(1 − σ 2 ) ; moreover β, γ, and µ ≡ µ 0 dκ dξ represent the dispersion, cubic nonlinearity and shoaling coefficient, respectively. The first two are simply the coefficients of the NLSE on arbitrary depth, see [10], and are functions of κ only; detailed expressions can be found in Appendix A; µ results from wave-energy conservation arguments as µ 0 ≡ 1 2ωcg d [ωcg] dκ , i.e., µ is the logarithmic derivative of c g . At variance with [27], we include also a homogeneous loss term, ν due to, e.g., viscosity or friction with bottom and sidewalls, which is appropriate at the NLSE order [28,29].
Let g = ω = 1 for definiteness. It is well-known that β < 0 for all values of κ (provided that only surface gravity waves are considered) [blue solid line in Fig. 1(a)], while γ ≥ 0 for κ ≥ 1.363 [red dashed line in Fig. 1(b)]. Recall also that c g [red dashed line in Fig. 1(a)] is maximum for κ ≈ 1.20.
The form of µ allows us to simplify Eq. (2); let U = V exp − ξ 0 µ(y) dy − νξ . Eq. (2) can be rewritten as i.e., a NLSE with varying parameters, with The effect of shoaling is clear from Eq. (4): in the focusing regime, βγ < 0, it slightly increases the effective nonlinearity, because c g monotically decreases, see the red dashed line in Fig. 1(a). The effect of ν is to decrease the impact of nonlinearity as the wave propagates. It is easy to verify that the perfect compensation of ν by shoaling is impossible for increasing depth. For the sake of simplicity, we will take ν = 0 in what follows, apart for the last part of Sec. IV.
In the framework of field theory, Eq. (3) conserves the total mass N ≡ ∞ −∞ |V | 2 dτ and the momentum P ≡ Im ∞ −∞ V * ∂V ∂τ dτ . We use them in our numerical calculations to assure the precision of solutions. No other conserved quantity is present, if coefficients have no specific functional dependence.

B. Modulation instability
Eq. (3) possesses a steady-state solution V s (ξ) = V 0 exp −iV 2 0 ξ 0γ (y)dy . For βγ < 0, this solution is modulationally unstable for a detuning Ω ∈ [0, Ω C ] from the central frequency ω, with Ω C ≡ 2 γ β V 0 . The linear MI gain is g = |βΩ| Ω 2 C − Ω 2 , with peak at Ω M ≡ ΩC √ 2 . This is the result of the conventional linear stability analysis, but it can also be thought of as the nonlinear phase-matching condition between the steady-state solution and the two sidebands, a sort of nonlinear resonance condition.
The main parameter of our problem,γ/β, is shown as blue solid line in Fig. 1(b). For comparison, we also include the ratio γ/β, as a dotted black line, to show that the effect of shoaling on MI is quite small (less than 5%) in the focusing regime. As this parameter is changed the same sideband frequency can turn from modulationally stable to unstable or experience a different instability gain along the MI curve. In Fig. 1(b), it is clear that the range of variation is quite limited, compared to optical fibers, because both β and γ tend to their deep-water limits as κ → ∞. The choice of the reference value κ = 2 (marked in Fig. 1) is a good trade-off for having strong enough nonlinear effects, while avoiding high-order corrections appearing whenγ ≈ 0, see for instance [30,31].
The MI gain is a linear approximation, beyond which the nonlinear behavior demands a more detailed analysis.

C. Nonlinear regime
A thorough understanding of the problem can come from a low-dimensional analysis. We follow the threewave truncation proposed in Ref. [18], that was proven effective also in higher-order generalizations of the NLSE [19,32].
where Ω is the angular detuning in normalized units, and A n , with n ∈ −1, 0, 1 are complex variables, the phases of which are denoted by φ n . It is easy to reduce Eq. (3) to a one degree-of-freedom (d.o.f.) Hamiltonian system [18]. The canonical variables are the conversion rate to sidebands η ≡ |A1| 2 +|A−1| 2 E and the relative phase ψ ≡ is a conserved quantity of the truncated system, as well as the sideband imbalance χ ≡ |A 1 | 2 − |A −1 | 2 . Compared to [18], we consider a slightly different set of variables, more suitable to our goals.
The Hamiltonian function is where the prime denotes the derivative with respect to ξ. More details are found in Appendix B.
A final transformation to X ≡ E ξ 0γ (y)dy allows us to simplify the Hamiltonian function to with α ≡ − βΩ 2 γE + 1 = Ω ΩM 2 − 1 = −4a AB + 1, with a AB the well known parameter of the AB. Noẇ where the dot denotes the derivative with respect to X.
The system is modulationally unstable for |α| ≤ 1, the peak gain is for α = 0, while the MI cut-off is for α = 1.
The two different topologies of the phase-plane (for α ≷ 1) are exemplified in Fig. 2  For what follows, it is also useful to recall that, for α < 1, the trajectory emanating from (ψ 0 ,η 0 ) is homoclinic and is referred to as a separatrix. By direct inspection of Eq. (7), it is easy to see that trajectories always turn anticlockwise for X > 0, irrespective of α. This implies that the separatrix exits the origin in the second and fourth quadrants and rejoins it in the third or first, respectively. We recall also that conventionally, trajectories outside (resp. inside) the separatrix are named period-two (resp. one) solutions. This is apparent in Fig. 2(b) and corresponds to the classification of NLSE solutions being either periodic or double-periodic [33]. The separatrix turns out to correspond to an AB, while the centers (ψ 2 ,η 2 ) to the steady state dn-oidal solution [34] It is also important that H (X) (ψ 0 ,η 0 ) = 0 for all α. For α ≥ 1, H (X) > 0 everywhere in the whole unit disk, while for α < 1, H (X) ≷ 0, outside or inside the separatrix, respectively. This is obvious, by noticing that In general, as the bathymetry and thusγ/β vary, the change of α lets H (X) (or H (ξ) ) vary across 0, see Eq. (6). This additional degree of freedom provides the flexibility to explore the stabilization regime we will present in the next section.
We will refer to the results of the present section as truncated or three-wave model, while the solutions of Eq. (3) are referred to as simulations.

III. STABILIZATION OVER AN UNEVEN BOTTOM
It is well known from classical mechanicsthat a trajectory oscillating around an elliptic fixed point keeps on following the same type of oscillatory trajectories if an internal parameter is changed adiabatically, i.e., the speed of variation is much smaller than the oscillation frequency [35]. For Hamiltonian system, a quantity, called the adiabatic invariant, is conserved all along the transition; this is the classical counterpart of Ehrenfest theorem in quantum mechanics. In order to solve our problem, we have to go beyond this result and recall the theory of autoresonance [21,22,36]. Two possible regimes can occur. Either the trajectory starts close to an equilibrium and a parameter is changed adiabatically, so that the adiabatic invariant is conserved; or it is forced to cross the separatrix and phase-locks in the close proximity to an equilibrium and the adiabatic invariant is not conserved. We explain below that the second solution is much more practical if the total transition length is constrained and for the flexibility in initial conditions. Thus, we focus here on how to physically apply the second approach to our model. As we showed above in Sec. II C, our system has only elliptic fixed points for α > 1 and both unstable hyperbolic and elliptic points for α < 1. Our aim is to stabilize the trajectory around (ψ 2 ,η 2 ) starting from small oscillations aroundη 0 , by varying α. The trajectory must thus cross the separatrix: a sign change of the Hamiltonian is associated to this transition.
Thus, three different aspects have to be considered: (i) the initial stage where the system behaves almost linearly, (ii) the separatrix crossing stage, and (iii) small oscillations around an equilibrium adiabatically shifted towards a larger η. We describe the three successively below.

A. Linear stage
We start from η 0 ≡ η(0) ≪ 1 and α 0 ≡ α(0) > 1, and linearize the system of Eq. (7) in order to understand its behavior when we tune the parameters to cross the bifurcation point α = 1 from above. By letting R ≡ √ ηe iψ , we reduce Eq. (7) tȯ which can also be obtained by linearizing the complex system reported in App. B, Eq. (B1), directly. The validity of Eq. (8) is limited to η ≪ 1; nevertheless, we can obtain some useful information about the full dynamics. We let R = u + iv and split Eq. (8) in real and imaginary part to getu If we divide these two equations term by term, we see that the solution is of the form v 2 = C − α−1 α+1 u 2 , which is either an ellipse or a hyperbola, for resp. α ≷ 1. For α > 1, it entails periodic oscillations, albeit, Λ lin (α) > Λ nl (α, H (X) ), defined as the periods predicted by Eqs. (8) and (7) respectively, see App. C and D. For α < 1, Eq. (8) gives exponentially divergent solutions. For α = 1, we have a pair of straight lines v = ± √ C, i.e., the horizontal semi-axis of the ellipse diverges. If |u| ≫ C at the same X, ψ → mπ =ψ 2 . Thus we can define this stage as the phase-locking stage.
The trajectories of the full nonlinear system turn anticlockwise, so do necessarily the solutions of its linearized version [the first of Eqs. (9) clearly shows that]. In order to follow the separatrix and then cross it and approach the centers located at ∓η 2 = ∓ 2 7 (1 − α), (u, v) are required to lie in the second or fourth quadrant, respectively: at α ≈ 1, we thus imposeuu > 0 (or, equivalently, uv < 0). Otherwise, the solution moves away from the elliptic fixed points and oscillates outside the separatrix.
In order to find suitable initial conditions, we resort to a local approximation in power series, shown in App. E. We conclude that, for α 0 close to 1 and ψ 0 = ± π 2 trajectories evolve to the correct quadrant and phase-lock to, respectively, π or 0, while ψ 0 = 0 does not.
A lower limit to ∆α i must be imposed. α(X * ) = 1 gives X * = α0−1 ∆α i . We require that X * ≪ Λ nl /2, i.e., the MI band is crossed before the system reaches the peak η. Otherwise, the trajectory would point back and could not enter the separatrix as this last appears. We conclude that ∆α i ≫ 2(α0−1) , as shown in App. C and D.
In order to lie close but near the separatrix as it appears, we require η 0 ≪ 1. The Hamiltonian takes thus the value H (X * ) ≈ 1 + (X * ) 2 v 2 0 at the bifurcation point.

B. Intermediate regime
Suppose that the the solution of Eq. (7) behaves at X * as a trajectory close to the separatrix, Eq. (C6) in App. C, and grows away fromη 0 . After an initial exponential growth, η slows down and its growth rate starts soon decreasing. The homoclinic orbit appears at α = 1 and expands linearly in width with decreasing α. From Eq. (6), as αη decreases, H (X) will change sign, thus separatrix crossing occurs. The analytic treatment to characterize the solution near this point is very involved for the system given by H (X) [37] and does not provide hints about the dynamics of Eq. (3). Nevertheless, we estimate the optimal variation of α and the distance at which it can be achieved by following a simpler argument, similar to what reported in Ref. [38]. Starting at X * , the optimal transition is such that H (X) (X * * ) = H min , where X * * marks the adiabatic stage start. In this way, the orbit reaches closely to (ψ 2 ,η 2 ). From Eq. (6), we have where we assumed, as before, that α decreases linearly with slope ∆α t . We can thus approximately integrate Eq. (10) and write η * ≡ η(X * ) is known from the linear stage above, and we take η * * ≡ η(X * * ) =η 2 , to enforce the proximity to the center at some given distance. We thus require that i.e., the start and end of the intermediate stage are separated by roughly a fourth of a period of an external orbit close to the separatrix, computed at X * , see App. C. This is justified by the fact that we have period-two solutions outside the separatrix. By plugging these values into Eq. (11), we obtain the optimal slope for changing α in the intermediate stage, We notice that the farther we start from the separatrix, the larger the variation of α is required.

C. Adiabatic conversion stage
Suppose that the separatrix is crossed and, at distance X * * , the system is close to the center (η 2 ,ψ 2 ) computed at the current value of α(X * * ). Suppose we can approximate α(X) = α(X * * ) − ∆α f (X − X * * ). The trajectory will keep on oscillating around the equilibrium, which in turn varies with α, provided that an adiabaticity condition on ∆α f is satisfied. We estimate it by resorting to the same approach of Ref. [36].
It is easy to check that κ 2 grows monotonically for − 3 2 < α < 1, thus the most stringent upper bound on ∆α f occurs at X * * .

A. Initial conditions
We suppose for simplicity that κ is changed linearly all over the domain: α 0 = 1.56, i.e., Ω = 1.6Ω M , and κ varies from 2 to 5. In practice, this means a linear variation of h, see the magenta dashed dotted line in Fig. 3(a) (the scale on the right axis). The effect on α is instead a faster variation in the beginning and slower after ξ ≈ 200. This is a particularly favorable situation for the locking into the elliptic fixed point, according to the previous discussion.

B. Simulation results
In Fig. 3, we clearly identify the three stages described above: (i) the linear, around α = 1, where η grows and ψ approaches π (red-shaded area); (ii) the intermediate, starting at ξ ≈ 200, where the growth slows down, the separatrix is crossed and ψ locks to π; (iii) and the adiabatic, starting at ξ ≈ 450, where η adiabatically follows the equilibrium up to η ≈ 0.5 (green-shaded area). The residual oscillations in both amplitude and phase are below 5% and 1% in relative terms, see blue solid lines in panels (a) and (b), respectively. The second-order sideband fraction, defined as η (2) represents less than 10% of the total mass N [red dashed line in panel Fig. 3(a)]. They are generated via nonlinear processes of the sort 0 ± Ω ± Ω → ±2Ω, which are thresholdless and oscillating. They partially account for the discrepancy between the numerical solution andη 2 (black dotted line). Systematically simulations end up oscillating around a larger η than predicted by the threewave truncation [18,19].
We compare simulations (solid line with changing hue) to the truncated model (dashed blue line) in the phaseplane, Fig. 4. In both cases η grows, the phase is locked and the residual oscillations are very small. Notably in the simulation, the oscillations around the average are limited to less than 0.025.
The asterisk marks the α = 1 transition, after which the linear approximation soon breaks down. The circle denotes instead the separatrix crossing, H (ξ) = 0. Notice that the trajectory turns away from the horizontal axis just after a close approach to an elliptic equilibrium (equivalent to (ψ 2 ,η 2 )). This occurs at η = 0.28, and the phase is then locked, see Fig. 3(b).
The three-wave solution (dashed line in Fig. 4) exhibits larger oscillations than the simulated ones (the horizontal and vertical scales differ much): in fact, the final value of κ is chosen to minimize these latter. The former meets its optimal conversion effectiveness at κ(ξ = 1000) ≈ 5.5, which combines the fast locking condition with the adiabatic following of the center: ∆α i ≫ 0.5 at ξ = 0, while ∆α f ≪ 5 at ξ = 450. For such a κ, the simulation turns out to oscillate more, which we explain by the faster displacement of the elliptic fixed point of the NLSE compared toη 2 , i.e., the conditions (12) and (14) are stricter for the NLSE than for the truncated model. This is again inherent to the three-wave approximation. In principle there is no limit on how large the fraction of N can be funneled into η. The physical range ofγ/β is nevertheless limited, see Fig. 1(a). Finally, once the total length of the system is constrained, though, ∆α f is bounded. The condition of Eq. (12) looks quite more stringent, but we verified numerically that, provided the separatrix is crossed, the behjaviour is very similar to the optimal one: the blue dashed line in Fig. 3, pertaining to the three-wave model, shows indeed the typical behavior.
A third alternative representation is available. Recall that Eqs. (5)-(7) are equivalent to a particle moving in a potential well, as explained in App. C. Notice that the potential well W (η) depends on the initial value of H (X) , thus its minima do not correspond to equilibria, in general. We let H (X) vary and update it at each integration step, according to the X-dependent expression Eq. (6), see App. C for more details. In Fig. 5 we show the map of the accessible values of the potential well −W (η) ≤ 0: this is very shallow and narrow at the beginning (where α > 1), then it becomes broader and deeper. Again, we see that the linear approximation, dark green dotted line, diverges at ξ ≈ 200. After this linear stage, the well smoothly widens and deepens. In the last stage, from ξ ≈ 450, the potential well gets deeper and deeper and the results of the simulation (blue solid line) is clearly trapped into it, as expected by the adiabatic following of the elliptic fixed point, proven above, and in spite of the systematic difference with the three-wave results.
Finally, a further limitation inherent to water waves is that nonlinear effects cannot be increased arbitrarily, because they scale as ε 3 , an AB envelope peaks are roughtly two-to-three times the background amplitude, and wavebreaking occurs if ε 0.4 [39]. The physical soundness of our approach is confirmed by representing the evolution of ε attained by U -the solution of Eq. (2), which represents the envelope of physical surface elevation, see Fig. 6. We notice that the proposed stabilization technique almost completely suppresses oscillations of U ; this reflects in negligible ε overshoots, never larger than 0.3, which guarantees that the train of pulses wont break.

C. A glimpse into a physical realization in hydrodynamics
In the previous section we use a quite conservative set of parameters, in order to assure the validity of the NLSE and the non-breaking of the wavetrain. The question arises if the stabilization can be achieved in a laboratory setting.
We consider a 100 m long wavetank. Our ξ ∈ [0, 1000] domain can be rescaled to it, once the carrier frequency is f = 1.  This is an idealization, because damping occurs. From Eq. (4), we notice that the shoaling partially compensates dissipation. Nevertheless, shoaling becomes negligible for larger κ, while the wave field keeps on damping exponentially: mathematically it is impossible to have exact compensation because c g decreases with κ. Moreover, we showed above that phase-locking is kept only if α is changed slowly. Thus, it is not possible to simply choose an arbitrary larger κ so thatγ reaches the same values of the undamped case; in fact, this will lie outside of the accessible parametric range because, for κ > 5, γ is almost constant, see Fig. 1.
The pulse train is thus meta-stable: for large enough damping, the separatrix will be crossed again, a periodtwo solution will be observed, and eventually the wave will vanish completely [40,41].
The analytic treatment is as involved as the one required to describe the second stage of the stabilization.
We found from simulations that keeping every parameter as before a total loss of 20% can be tolerated. For the wavetank length specified above, this corresponds to ν ≈ 2 × 10 −3 m −1 , which is reasonable when the effect of sidewall dissipation is taken into account [42].

V. CONCLUSIONS
In this work we study the nonlinear stage of evolution of modulational instability in surface water-waves over a water body of gradually increasing depth. We show that this stage can be stabilized and results in a uniform train of pulses on a background. The initial condition needs not be restricted to an exact NLSE solution (e.g., an Akhmediev breather), but just a harmonic perturbation with a given small amplitude.
Based on a three-wave truncation, we studied how a linear depth change naturally leads to a virtually frozen state (which can be considered close to dn-oidal solution of the NLSE), provided that suitable initial conditions (frequency lying just outside the instability margin and with a relative phase facilitating separatrix crossing) are chosen.
Within these restrictions, still a wide range of carrier frequencies and depth variation allows one to achieve the stabilization, even in spite of the unavoidable viscous damping.
Although the flexibility available to vary parameters in the hydrodynamics of surface water-waves is much less than in other physical systems, such as optical fibers, our results will help clarify the possibility to dynamically control the breathing evolution of water wave-packets and to understand the impact of bathymetry on the persistence (or lifetime) of rogue waves.

ACKNOWLEDGMENTS
We acknowledge the Swiss national Science foundation (SNF grant 200020 175697) and the University of Sydney-University of Geneva Partnership collaboration award.
Then, let A n = √ ζ n exp iφ n , with ζ n and φ n real functions. By replacing these variables in Eq. (B1), we notice that φ n appear only in the relative phase ψ defined in the main text. Moreover, it is easy to observe that the total intensity E ≡ ζ 0 + ζ 1 + ζ −1 as well as the sideband imbalance χ ≡ ζ 1 − ζ −1 are conserved. It is thus practical to define η as in the main text, so that η ∈ [0, 1].
Appendix E: Local solution of Eq. (9) We consider the behavior of Eq. (9) around α = 1. At this point the evolution of v has an essential singularity; nevertheless, we numerically find that the solution is regular and this fact is the key to phase-locking.
Let α = 1 − ∆α i (X − X * ), with ∆α i > 0. At distance 0 < X = X * ≪ Λ (0) lin we reach the MI band edge Ω C . We look for a solution of Eq. (9) of the kind v(X) = ∞ n=0 a n (X − X * ) n and u(X) = with a 0 and b 0 arbitrary constants. From these expressions, it is easy to verify that the phase-locking conditions stated in the main text-uv < 0 for X X * -are equivalent to a 0 b 0 < 0. Indeed, v(0) = a 0 is a positive minimum (negative maximum), for b 0 ≶ 0. The other extremum of v if for X − X * = b0 2a0 < 0. As far as u is concerned, it has a single maximum (minimum) at X − X * = 2 ∆α i > 0. Now, we can find the best initial conditions for achieving phase-locking. If u(0) = 0 and v(0) = v 0 , we obtain a 0 ≈ v 0 and b 0 ≈ −v 0 X * : a 0 b 0 < 0. If, instead, u(0) = u 0 and v(0) = 0, we obtain a 0 ≈ ∆α i u 0 (X * ) 2 and b 0 ≈ u 0 : a 0 b 0 > 0 and the conditions for crossing the separatrix are violated. The approximation signs are valid if ∆α i (X * ) n ≪ 1, for n ≥ 1. Finally, notice that tan ψ(X * ) = a0 b0 , thus tan ψ → 0, i.e. phase-locked trajectories, only for the former condition. Only the unshaded regions inside the unit circle correspond to trajectories leading to phase-locking.
We graphically illustrate these results in Fig. 7. We show two different trajectories v 2 + α−1 α+1 u 2 = C, with C = η 0 = 0.025 and α ∈ {1.1, 1.01}. The orbits continuously move from one ellipse to another of bigger horizontal semi-axis. In order for the initial conditions to permit phase-locking, we require that they cross into the second or fourth quadrants before α = 1. This intuitively justifies also the lower bound on ∆α i discussed in the main text.
An alternative local solution is to consider a second order equation for u, which reads u +α α + 1u + (α 2 − 1)u = 0, and gives the same expressions of Eq. (E1).