Eckhaus instability of stationary patterns in hyperbolic reaction–diffusion models on large finite domains

We have theoretically investigated the phenomenon of Eckhaus instability of stationary patterns arising in hyperbolic reaction–diffusion models on large finite domains, in both supercritical and subcritical regime. Adopting multiple-scale weakly-nonlinear analysis, we have deduced the cubic and cubic–quintic real Ginzburg–Landau equations ruling the evolution of pattern amplitude close to criticality. Starting from these envelope equations, we have provided the explicit expressions of the most relevant dynamical features characterizing primary and secondary quantized branches of any order: stationary amplitude, existence and stability thresholds and linear growth rate. Particular emphasis is given on the subcritical regime, where cubic and cubic–quintic Ginzburg–Landau equations predict qualitatively different dynamical pictures. As an illustrative example, we have compared the above-mentioned analytical predictions to numerical simulations carried out on the hyperbolic modified Klausmeier model, a conceptual tool used to describe the generation of stationary vegetation stripes over flat arid environments. Our analysis has also allowed to elucidate the role played by inertia during the transient regime, where an unstable patterned state evolves towards a more favorable stable configuration through sequences of phase-slips. In particular, we have inspected the functional dependence of time and location at which wavelength adjustment takes place as well as the possibility to control these quantities, independently of each other.


Introduction
Self-organizing patterns are quite ubiquitous in nature and the study of the associated nonlinear spatial processes has become nowadays a sub-area of complexity science. Despite pattern formation is observed in different fields, such as biology, ecology, chemistry and medicine, the underlying phenomena share some common features [1,2]. For instance, it is widely accepted that regular periodic stationary patterns, that in the one-dimensional case are variously referred to as bands, rolls or stripes, arise from Turing-like instability. It results from the coupling between local interactions and dispersal, even in the absence of heterogeneity of the environment in which the process takes place. In some contexts, the large spatial scale and/or the slow evolution of patterns prevent replicability under laboratory conditions, so that theoretical investigations based upon mathematical models represent an efficient and inexpensive tool to acquire a rigorous understanding on the complex spatio-temporal dynamics. Mathematical models have indeed allowed to get valuable insights on many aspects of pattern dynamics: formation of regular and irregular mosaics, scaling phenomena, spatial modulation of pattern amplitude, transitions from a uniform state to a periodic one as well as between different periodic states [3][4][5][6][7][8]. In particular, it is known that, in a large domain, the equation ruling the spatio-temporal evolution of pattern amplitude close to onset is in the form of a Ginzburg-Landau (GL) equation [9][10][11][12]. This equation also allows to inspect the dependence of the stability of the resulting periodic states on wavelength. Indeed, patterns may undergo different destabilization mechanisms associated with changes in wavelength and amplitude, one of which is the Eckhaus instability (EI) [3][4][5][6][7]. EI acts on the roll phase to compress or dilate pattern wavelength and takes place when a given roll wavelength cannot be accommodated by the environment so that some of them are eventually created or eliminated. The process leading the system to be rearranged in a more favorable configuration is referred to as phase-slip [3,5,6].
A large body of literature has successfully addressed the bifurcation analysis of EI in both infinite and finite spatial domains and has also addressed the study of complex intriguing dynamics, such as those associated to quasi-periodic solutions and homoclinic snaking bifurcation structures [3-6, 8, 13-18]. However, none of those works has inspected the EI instability of stationary and quantized Turing patterns in the context of hyperbolic models. As widely reported in previous works [19][20][21][22][23][24][25][26][27], hyperbolic systems bring many advantages. Firstly, they take inertial effects explicitly into account and, thus, allow to overcome to paradox of propagation of disturbances at infinite speed, typical of parabolic models. They also appear more suitable to describe transient phenomena characterized by waves evolving in space over a finite time. Moreover, despite the mathematical convergence of parabolic and hyperbolic models is expected in the long-time limit due to the stationary nature of the excited patterns, the coupling between hyperbolicity and nonlinearity may generate richer transient dynamics.
In particular, the present manuscript addresses the analysis of one-dimensional stationary Turing patterns emerging in a class of two-compartments hyperbolic reaction-diffusion models, where both species undergo self-diffusion, defined on a large finite domain. In some previous works, linear and weakly-nonlinear stability analyses on the homogeneous steady states were carried out over small (in the presence of cross-diffusion [25]) or extended domains (with constant [19] and non-constant [20] inertial times). Here, we aim at addressing the stability and quantization of spatially-periodic patterns, in both supercritical and subcritical regimes.
With this in mind, we aim at deducing, first, the equation governing the evolution of pattern amplitude close to the onset of supercritical and subcritical bifurcation. To this aim, we apply multiple-scale weakly nonlinear expansion which is pushed up to the fifth perturbative order. Depending on the perturbative order at which the expansion is truncated, this procedure allows to build a cubic or cubic-quintic GL equation with real coefficients. Moreover, taking into account the quantization of modes due to the spatial confinement, we address a bifurcation analysis of the EI to accurately describe the existence and stability thresholds of all periodic branches [4]. This investigation is finalized at inspecting the qualitatively different dynamical features captured by the above GL equations.
Theoretical predictions here developed are then validated through a comparison with numerical simulations. To this aim, an illustrative example in which EI takes place in the context of dry-land ecology is here considered: the hyperbolic generalization of the modified Klausmeier model [19,20,[28][29][30][31]. It represents one of the simplest two-component systems able to predict the formation of vegetation bands in flat arid environments as the result of competition between soil water and vegetation biomass.
The comparison between analytical and numerical results made it possible to address several issues: (i) to validate theoretical predictions on bifurcation analysis developed in the frameworks of cubic and cubic-quintic real GL equations, in both supercritical and subcritical regime; (ii) to investigate how hyperbolicity affects not only transient regime dynamics between different patterned states but also the occurrence of phase slips observed during an Eckhaus-driven restabilizing process; (iii) the feasibility to control, independently of each other, the time and the spatial location at which such phase slips occur.
The manuscript is organized as follows.
In Sect. 2, we present the class of hyperbolic reaction-diffusion models under investigation and develop linear and weakly nonlinear stability analysis. The former allows to deduce the Turing instability threshold, the latter is used to deduce the envelope equations at different perturbative orders.
In Sect.3, we address the bifurcation analysis of the EI. First, we briefly review the most relevant results associated to the supercritical regime in the framework of the cubic GL equation. Then, we inspect the subcritical regime via both cubic and cubic-quintic GL equations. In particular, in all cases, we deduce the expressions of the most relevant dynamical features associated to primary and secondary quantized periodic states: stationary amplitude, existence and stability thresholds and linear growth rate. Finally, this section is concluded with a focus on the difficulties encountered in the theoretical prediction of the phase slip phenomenon.
In Sect. 4, we first validate the theoretical predictions developed in the previous sections by means of a comparison with numerical simulations carried out on the hyperbolic modified Klausmeier model. Then, we address further numerical investigations to inspect the spatiotemporal dependence of phase slip and to emphasize the role played by inertial times.
Concluding remarks are given in the last section.

Amplitude equations for hyperbolic reaction-diffusion models
We consider a class of 1D hyperbolic reaction-diffusion models for two species u(x, t) and w(x, t), both undergoing self-diffusion. We denote the w-by-u diffusion ratio by d, the kinetic terms by f (u, w) and g(u, w), the phenomenologically-introduced arbitrary constants by γ 0 and μ 0 which are related to the inertial times τ u and τ w of the two species by γ 0 = 1/τ u and μ 0 = d/τ w . According to Extended Thermodynamics (ET) theory [32], the diffusive fluxes J u (x, t) and J w (x, t) are considered as additional field variables obeying thermodynamically-consistent balance equations that, in the parabolic limit approximation, τ u → 0 and τ w → 0, reduce to classical Fick's law. In vector form, this model reads: being: where subscripts denote partial derivative with respect to the indicated variable. In the literature, the hyperbolic structure (1), (2) has been successfully employed in different contexts, such as plant ecology [19], epidemiology [33], air pollution [34] and chemistry [35]. In order to deduce the conditions for Turing instability let us briefly address linear stability analysis. Let B be a control parameter and U * = (u * , v * , 0, 0) a positive spatiallyhomogeneous steady-state satisfying N(U) = 0. By requiring that U * is stable against spatially-uniform perturbations but unstable with respect to non-homogeneous ones, we get the following restrictions: where the asterisk denotes that the function is evaluated at U * . From (3), it can be easily deduced that the critical values of control parameter B c and wavenumber q c at the onset of Turing instability are obtained by solving: Note that, since we deal with the formation of stationary patterns, the hyperbolic structure of system (1), (2) does not affect the expression of the critical parameters at onset, so that the occurrence of Turing instability is ruled by the same conditions found in classical parabolic models [2]. To describe the spatio-temporal evolution of pattern amplitude close to Turing threshold (4), (5), multiple-scale weakly nonlinear analysis is now addressed [3,6,8,36,37]. To this aim, we consider a small dimensionless parameter , expand the control parameter B around B c and the field U = U − U * as and introduce different time and spatial scales as follows: After substituting the above expansions into the governing system (1), (2), applying zeroflux boundary conditions over the physical domain x ∈ [0, D] and collecting terms of the same orders of , a set of linear equations for the U i is obtained: where the vectors F j ( j = 2, . . . , 5) are given by being ∇ = ∂/∂U, together with (10) Moreover, for a generic vector H, the expression (H · ∇) ( j) stands for the operator applied j times. The set of equations (8)-(10) has to be solved sequentially, as sketched in the Appendix. The removal of secular terms at O( 3 ) leads to an envelope equation ruling the spatiotemporal evolution of the pattern complex amplitude (X , T 2 , T 4 ) that takes the form of a real cubic GL equation: which preserves the structure found in classical parabolic models [3,4,6,11]. It is worth noticing that the real coefficients σ , L and ν inherit the dependence on the inertial times as follows: where the functions , and depend upon the w-by-u diffusion ratio d and the kinetic terms, together with their partial derivatives, but are independent of the inertial times. The explicit expressions of these functions are given in the Appendix for the illustrative example considered in Sect. 4. Therefore, the hyperbolic structure of the system may affect the spatiotemporal evolution of pattern amplitude. In particular, the growth rate σ and coefficient ν, that in the pattern forming region are always positive, play a significant role during transient regime. On the other hand, the sign of the Landau coefficient L determines the nature of the dynamical regime: L > 0 corresponds to the supercritical regime, whereas L < 0 to the subcritical one [19].
Remark It should be remarked that, to the best of our knowledge, weakly nonlinear analysis in hyperbolic reaction-diffusion systems on large finite domains has been pushed up to the third perturbative order [19,25], allowing to provide a suitable description of supercritical dynamics only.
From the removal of secular terms at O( 4 ), the following compatibility condition for the spatially evolution of the pattern amplitude is obtained [13]: Therefore, to investigate subcritical dynamics, we push weakly nonlinear analysis up to the fifth perturbative order where the removal of secular terms leads to the following envelope equation for the pattern amplitude: which is in the form of a real cubic-quintic GL equation. In (15), ∂/∂T = ∂/∂T 2 + 2 ∂/∂T 4 , the coefficients here appearing are real and represent second-order corrections of the coefficient involved in the cubic GL equation (12), namely σ = σ + 2 σ , L = L + 2 L, ν = ν + 2 ν and R = 2 R. Moreover, each of these coefficients encloses a non trivial dependence on the inertial times which, acting as additional degrees of freedom, may offer a richer scenario of spatio-temporal dynamics with respect to the parabolic counterpart. This statement holds true despite hyperbolic and parabolic models share the same structure of weakly inverted bifurcations to a stationary spatially-periodic state [14-17, 38, 39]. However, due to the cumbersome expressions here involved, conclusions can be only given through numerical simulations. For this reason, as well as to keep the length of the manuscript within reasonable limits, the full expressions of all the real coefficients involved in Eqs. (12)-(15) will be explicitly given in the Appendix for the illustrative example carried out in Sect. 4.

Bifurcation analysis of the Eckhaus instability
The analysis carried out in this section focuses on the formation and stability of stationary 1D Turing patterns originating in the hyperbolic reaction-diffusion system (1), (2) over a large finite domain D ∝ 1/ . In detail, the occurrence of the phenomenon of Eckhaus Instability (EI) will be analytically investigated in both supercritical and subcritical regime, with the main aim of quantifying some key dynamical features: stationary amplitude of patterns, existence and stability thresholds of each periodic bifurcating branch and linear growth rate.
To address how a non-favorable wavelength may lead to pattern instability, we adopt the standard procedure adopted in the literature [4,6].

Supercritical regime (a brief review)
To describe the spatio-temporal evolution of pattern amplitude close to onset of Turing instability in the supercritical regime, it suffices to push weakly nonlinear analysis up to the third perturbative order, where the cubic GL equation (12) is retrieved. We remind that the coefficients σ and ν are always positive whereas the Landau coefficient L is positive in the supercritical regime only. To reduce the number of coefficients here involved, let us apply the following rescaling: that allows to recast Eq. (12) as: For simplicity, let us drop the tilde notation and assume that rolls take the structure: where and Q describe, respectively, the amplitude and the phase of the envelope whereas θ is an arbitrary constant, to be determined according to boundary conditions. Substituting (18) into (17) gives, apart from the null solution = 0, the stationary amplitude: Patterned solution (18), (19) exists for ζ > Q 2 and is referred in the literature to as pure mode [4]. It has to be distinguished from the trivial solution = 0 which is termed conductive mode. Note that, according to (16), (18) and the structure of U 1 (see (A.13) in the Appendix), the field has total wavenumber (Q + Q c ). The finite domain, however, implies quantization of modes, i.e. not all wavenumbers are admitted but only those integer ones satisfying (Q n + Q c ) ∈ Z, n ∈ N 0 . Let us call the quantized amplitude of the n-th mode, that exists for ζ > Q 2 n := ζ e,n . Zero-flux boundary conditions, together with quantization of modes, restrict the possible values of θ to 0 and π only. For simplicity, we set θ = 0. To address linear stability of patterns, let us apply small perturbations in amplitude and phase of rolls as follows: where |ξ |, |φ| 1. Substituting (21) into (17), keeping the linear terms in ξ and φ and taking real and imaginary parts, the system ruling the evolution of perturbations results: Then, by assuming: beingφ andξ arbitrary constants and the dimensionless wavenumber k = 0, the following quadratic equation for the linear growth rate λ is obtained: whose roots are: Since the eigenvalue λ (2) n is always negative, the stability of the associated mode depends on λ (1) n which, taking into account (20), can be expressed as: Consequently, stable roll solutions in a finite domain satisfy: the last equality of which identifies the well-known Eckhaus parabola in the (Q, ζ ) plane [4,6]. By means of this approach, we have thus retrieved a well-known result: the primary branch 0 is always stable at onset (ζ ≥ ζ e,0 ) while all the other branches n (with n ≥ 1) have always n unstable directions at onset and undergo n secondary bifurcations (of pitchfork type) in order to become stable. The final, restabilizing, bifurcation corresponds to the Eckhaus threshold [4]. The analysis developed in the supercritical regime, where pattern amplitude obeys cubic real GL equation (17) can be thus summarized in tabular form, see Table 1.
Moreover, taking into account (13) and (16), we can conclude that the stationary amplitude of each pure mode n , the existence thresholds ζ e,n and the Eckhaus thresholds ζ E,n are unaffected by hyperbolicity, being the ratio σ/ν independent of inertial times. This result is compatible with the stationary nature of the emerging spatially-periodic patterns.

Subcritical regime
In this section we aim at describing the pattern amplitude close to onset of a subcritical bifurcation via a real cubic (12) or cubic-quintic (15) GL equation. Although it is known that these amplitude equations may capture different dynamical features, the goal of the Table 1 Dynamical features predicted by cubic GL equation (17) in the supercritical regime Stationary amplitude of pure modes Range of stability subsequent analyses is to detect and emphasize such differences. We keep in mind that the Landau coefficient L is negative in this regime.

Cubic Ginzburg-Landau equation
Starting from cubic GL equation (12), let us rescale the variables as in (16) except for the amplitude:˜ This allows to recast the cubic GL equation as follows: Removing tilde notation and adopting similar arguments as those addressed in the previous subsection, the quantized pure modes can be expressed as: and exist for ζ < Q 2 n := ζ e,n .
For each Q n , the previous expression is representative of two periodic branches that originate from the conductive state through a subcritical pitchfork bifurcation. It is easy to verify that the expression of the eigenvalue determining stability of each bifurcating branch λ (1) n (given in (26)) as well as the condition on pattern stability (given in (27)) are formally unchanged with respect to those found in the supercritical case. However, a key point has to be stressed according to the analysis here carried out: the existence (31) and stability (27) conditions are simultaneously fulfilled in the subcritical regime for n = 0 only. Results are summarized in Table 2. Therefore, if the pattern amplitude obeyed cubic GL equation (29) close to the onset of a subcritical bifurcation, the sole primary branch 0 would undergo Eckhaus instability [40]. This marks a notably difference with respect to the behavior outlined in the supercritical regime.

Cubic-Quintic Ginzburg-Landau equation
Let us now inspect whether the real cubic-quintic GL equation (15), obtained by pushing weakly nonlinear expansion up to fifth order, may provide a qualitatively different description Table 2 Dynamical features predicted by cubic GL equation (29) in the subcritical regime Stationary amplitude of pure modes Range of stability of subcritical dynamics. The signs of the coefficients are assumed to be σ > 0, L < 0, ν > 0 and R < 0 in order to guarantee that the the primary branch exhibits a weakly inverted bifurcation at onset (i.e., it shows hysteresis) and that such a bifurcation saturates to quintic order [14]. By using the following rescaling of variables: the cubic-quintic GL equation (15) may be recast as: being ∂/∂t = ∂/∂t 2 + 2 ∂/∂t 4 . Notice that the coefficients ζ and ρ are real and positive. Dropping tilde notation, considering "perfect" rolls structure as in (18) and using similar arguments as those addressed in the previous subsections, the quantized pure modes originating from the turning point can be expressed as n = n e i Q n x with: .
The large amplitude branch exists for whereas the small amplitude one exists for Considering (21), (23) and (33), the eigenvalue λ (1) n determining stability of pure modes is given by that, using (34), allows to deduce the following stability condition for the large amplitude branch: the last equality of which defines the Eckhaus threshold predicted by cubic-quintic GL equation. From the comparison between existence (35) and Eckhaus thresholds (38), it is possible to identify which bifurcating branch may undergo EI as: Inequality (39) points out that the restabilizing mechanism depends on the order of the bifurcating branch. Indeed, for n = 0, since quantization of modes implies Q 2 0 − 1/4 ≤ 0, the necessary and sufficient condition for the restabilization of the primary branch becomes: that represents a restriction on the allowed values of the coefficient ρ. On the other hand, for n ≥ 1, since Q 2 n − 1/4 > 0, the inequality (39) is fulfilled for any ρ > 0, so implying that all the secondary branches exhibit restabilizing Eckhaus bifurcations.
According to this analysis, the existence and consequently the amount of secondary bifurcations originating from periodic branches requires: with k ∈ N. Since ρ > 0, the radicand is always greater than 1, so that k * n ≥ 2|Q n |. Notice that, fulfillment of conditions (39) and (40), as well as the numerical estimation of the quantity k * n in (41), depends on the value of the coefficient ρ and thus, it has to be checked numerically for the specific model under consideration. This issue will be addressed in the next section.
Therefore, the cubic-quintic GL equation predicts at least four key differences with respect to the cubic counterpart: (i) the range of existence of patterns is bounded from below ζ ≥ ζ e,n := Q 2 n − 1/(4ρ); (ii) such an existence threshold is smaller than the bifurcation point (ζ n = Q 2 n ) at which periodic branches originate from the conductive state; (iii) the primary branch may undergo EI provided that (40) is satisfied whereas (iv) all the secondary branches undergo EI and the number of restabilizing bifurcations depends on value of ρ through condition (41).
Results contained in this subsection, that represent one of the novelties introduced in this manuscript, are summarized in Table 3.
In the literature [14], Brand and Deissler tackled a similar study on Eckhaus instability and deduced the analog of the Benjamin-Feir-Newell instability criterion for a weakly inverted bifurcation. However, some differences with respect to the present work should be pointed out. Indeed, those authors addressed the analysis on infinite domains only, so that quantization of modes was not investigated. Furthermore, the stability domain for finiteamplitude plane-wave solutions of cubic-quintic Ginzburg-Landau equation that are stable against the Eckhaus sideband instability was provided in terms of wavenumber instead of control parameter.

Spatio-temporal dependence of phase slip
The analysis carried out so far has revealed that the transition from an Eckhaus unstable state towards a more favorable, stable, patterned configuration may occur under different dynamical regimes. This process involves a sequence of transient states during which the wavelength of patterns is adjusted via the formation of amplitude defects and the appearance Table 3 Dynamical features predicted by cubic-quintic GL equation (33) in the subcritical regime Stationary amplitude of pure modes n ≥ 0 Range of stability Unfortunately, the global evolution of a given perturbation cannot be correctly described by a local analysis around a steady state, as the one performed in the previous sections, so that the occurrence in time and the location in space of phase slips is hard to be predicted. To the best of our knowledge, some good approximations of these quantities have been obtained by means of local theories under the assumption that phase slip occurs 'just' after an ad-hoc choice of initial data [5,41], but no explicit expressions have been provided in the general, global, case. Moreover, some works [18,40,43] pointed out the possibility to control and vary the time to phase slip only if some 'free' model parameters, not involved in the spectrum of the linearized problem, are available. At the same time, those works inspected the dependence of the time to phase slip on the linear growth rate but no mention was given on the functional dependence of the spatial location at which phase slip occurs. Therefore, the issue of controlling, independently of each other, time and space at which phase slip takes place appears to be still unaddressed. This issue will be investigated numerically in Sect. 4.3.

An illustrative example: the hyperbolic modified Klausmeier model
In order to validate the theoretical predictions carried out in the previous sections, let us now investigate, as an illustrative example, the occurrence of EI in the context of dry-land ecology. To this aim, let us take into account one of the simplest, conceptual, two-compartment models describing the formation of vegetation stripes in arid environments, i.e. the Klausmeier model [28,37,44,45]. In particular, the analysis will be focused on its modified version [29][30][31] that has been introduced in the literature to model pattern dynamics over flat terrains. Moreover, motivated by the observation of inertial effects [46][47][48][49][50] and long transient regimes [51][52][53] in vegetation patterned dynamics, an hyperbolic generalization was proposed in [19,20].
In this model, the two species introduced in Sect. 2 are representative of plant biomass u(x, t) and soil-water w(x, t). It is also assumed that d 1 measures the water-to-plant diffusion ratio [30], τ u = 1/γ 0 and τ w = d/μ 0 denote the inertial times associated to vegetation and water, respectively. Furthermore, the model encloses a per-capita rate of water uptake proportional to plant biomass, the plant growth rate proportional to water uptake, a linear dependence of plant loss with strength B and a mean annual rainfall represented by the parameter A. According to the above assumptions, the hyperbolic reaction-diffusion model is in the form (1), (2) with kinetic terms given by: In line with the analysis carried out in Sect. 2 , plant mortality B is taken as the main control parameter as it encloses variability due to natural, human and herbivory effects. According to literature data [28,54], realistic values of plant loss and rainfall belong to the ranges B ∈ (0, 2) and A ∈ (0, 3), respectively. To evaluate the spatio-temporal evolution of pattern dynamics, the governing system is integrated numerically over a finite domain of length D by means of COMSOL Multiphysics [55]. Moreover, the MATLAB package PDE2PATH [56] is used to build up the bifurcation diagrams in both supercritical and subcritical regimes.
In the Appendix we report all the explicit expressions of the quantities arising from linear and weakly nonlinear analysis for the illustrative example here considered.

Supercritical regime
In the supercritical regime, the following model parameters are chosen: water-to-plant diffusion ratio d = 10 3 and rainfall A = 2.8. According to (4) The bifurcation diagram obtained in the supercritical regime is depicted in Fig. 1, where solid red (dashed black) lines are representative of stable (unstable) branches. In agreement with theoretical predictions (see Table 1), it consists of a primary branch (n = 0) and several secondary branches (the cases n = 1, 2, 3 are here shown). The primary branch bifurcates supercritically at B e,0 giving rise to two stable branches for B > B e,0 . The secondary branches originate from the conductive state at B e,n and, for n > 1, undergo (n−1) unstable (pitchfork) bifurcations until they re-stabilize at the Eckhaus threshold B E,n . The bifurcation diagram is computed for different values of inertial times achieving the same result, as expected from the analysis carried out in Sect. 3.
The numerically-computed values of the existence and Eckhaus thresholds are summarized in Table 4 and compared to the theoretical ones given in Table 1. As it can be noticed, the resulting agreement between them is satisfying for all the investigated branches. The agreement gets slightly worse for higher-order branches, as the dimensionless distance from the Turing threshold 2 increases, as expected from weakly nonlinear approximation. It should be also remarked that the following inequalities on existence and Eckhaus thresholds, B e,n > B e,n−1 and B E,n+1 > B E,n , hold for any n ≥ 1.
In order to check the validity of the bifurcation diagram and, in turn, to confirm the theoretical predictions carried out in Sect. 3.1, let us now integrate numerically the governing systems (1), (2) by varying the control parameter B. In all simulations, the initial condition is assumed to take the form given by (6) 2 and (A.13), where the pattern amplitude is taken as a small perturbation of the stationary value for each considered branch n = ζ − Q 2 n e i Q n x , in line with Eqs. (18), (20). Results are shown in Fig. 2.  Zero-flux boundary conditions impose that only integer values of the wavenumber are allowed. Therefore, the system will only support periodic patterns whose overall wavenumber Q c + Q n are integer. Since q c = 0.3814, according to (16), the rescaled value of wavenumber becomes Q c = 24.28 and the closest integer is 24, so that the finite geometry shall select a primary mode with wavenumber Q c + Q 0 = 24. Considering that π(Q c + Q 0 )/D = 2π/ , with the pattern wavelength, the whole computational domain D shall thus accommodate 12 . Indeed, by choosing a value B > B e,0 , starting from a small perturbation of the branch 0 , the confirmation of the stability of the primary branch is achieved, as can be seen from Fig. 2a. Analogously, starting from B e,1 < B < B E,1 and setting as initial condition a small perturbation of the branch 1 , the system undergoes a phase slip, during which the number of wavelengths is reduced, and finally it converges towards the stable 0 branch (Fig. 2b). On the contrary, by leaving the previous initial condition unchanged, if the control parameter is set just above the Eckhaus threshold, i.e. B > B E,1 , the system now stabilizes along the stable branch 1 . In this latter case, in fact, the next closest integer is 25, so the finite geometry shall select a mode with wavenumber Q c + Q 1 = 25 that corresponds to D = 12.5 (see Fig. 2c). Similar conclusions can be drawn for higher-order branches. In detail, Fig. 2d represents the evolution from a small perturbation of the branch 2 (characterized by Q c + Q 2 = 23 and D = 11.5 ) by using B e,2 < B < B E,2 . Indeed, in this range of the control parameter, the initial branch is unstable, so the system converges to the stable branch 0 and a phase slip at the boundary of the domain occurs. Instead, if the value of the control parameter is slightly increased to overcome Eckhaus threshold, B > B E,2 , the branch 2 now stabilizes (see Fig. 2e). Finally, the last two figures correspond to simulations where the initial conditions are set as small perturbations of branch 3 with control parameter slightly smaller (B e,3 < B < B E,3 , Fig. 2f) or larger (B > B E,3 , Fig. 2g) than the Eckhaus threshold of the corresponding branch. As can be seen, in the former case, the initial branch is unstable and system experiences a transition towards the stable 0 branch whereas, in the latter case, the system stabilizes along the branch 3 , consistently with theoretical predictions.

Subcritical regime
Numerical investigations carried out in the subcritical regime make use of the following parameters: water-to-plant diffusion ratio d = 10 3 and rainfall A = 0.02. These parameters provide a critical value of plant loss B c = 9.378 × 10 −3 , critical wavenumber q c = 5.62 × 10 −2 and Landau coefficient L = −5.69 × 10 −5 < 0. As can be noticed, the much smaller value of critical wavenumber suggests that the excited wavelengths are much larger than those observed in the supercritical regime, consistently with previous results [19]. For this reason, the computational domain has been enlarged to D = 1000. Moreover, since subcritical pattern dynamics are expected to occur over much longer timescales, the spatio-temporal evolution is observed over a wider time window t ∈ [2 × 10 5 ].
Two main considerations can be drawn from the bifurcation diagram obtained in the subcritical regime, shown in Fig. 3. First, the primary branch undergoes a subcritical bifurcation at B 0 (corresponding to ζ = Q 2 0 ) where two unstable periodic branches originate from the conductive states, become stable at the Eckhaus threshold B E,0 < B 0 and still survive for B > B 0 (ζ > Q 2 0 ). Second, the secondary branches (n ≥ 1) undergo a restabilizing bifurcation for B E,n > B n (ζ E,n > ζ n ).
Therefore, two striking features appear markedly in contrast with the analysis developed in the framework of cubic GL equation: the existence of periodic solutions beyond the predicted threshold and the occurrence of restabilizing bifurcations for higher order branches too. Evidently, weakly nonlinear analysis developed in Sect. 3.2 fails to describe this phenomenon as it is based on the assumption that pattern amplitude behaves as O( ) close to onset. The approach based on cubic GL equation appears to be, thus, not suitable to predict the occurrence of restabilizing bifurcations at large amplitudes, whose description requires weakly nonlinear expansion to be pushed forward to higher orders, as done in Sect.3.2.
The features exhibited by the bifurcation diagram in Fig. 3 are, indeed, fully compatible with results obtained in the framework of cubic-quintic GL equation (summarized in Table  3). Several considerations may be addressed on this point. First of all, the occurrence of Eckhaus instability on the primary branch at ζ E,0 < Q 2 0 confirms that the numericallycomputed value of the coefficient R is negative and fulfills the restriction (40). Moreover, the secondary branches originating from the unstable conductive state do exhibit Eckhaus re-stabilizing bifurcations for values of the control parameter larger than existence threshold, i.e. B E,n > B e,n . Furthermore, each n-th order branch undergoes exactly (n − 1) secondary unstable bifurcations (as in the supercritical regime), in agreement with the numerical values of the coefficient k * n that always fall in the range n < k * n < n + 1. Let us finally comment that the bifurcation diagram drawn in Fig. 3 has been computed for different inertial times (by varying them over different order of magnitudes) obtaining  identical results, so verifying the independence of existence and Eckhaus thresholds on inertial effects. The theoretical values of Eckhaus threshold of the first three branches originating from the conductive state arising from cubic and quintic analysis are reported in Table 5. The previous theoretical values are then compared to the numerical ones, obtaining a satisfying agreement in all cases. Of course, theoretical values deviate away from numerical ones as far as the dimensionless distance from the threshold is increased, as expected from weakly nonlinear approximation. Let us now inspect the spatio-temporal evolution of patterns resulting from variation of the control parameter B along the different branches characterizing the bifurcation diagram of Fig. 3. To achieve this goal, the governing system is integrated numerically by using the same boundary and initial conditions as the ones adopted in the supercritical regime. Results of this analysis are shown in Fig. 4. According to the previous values of critical parameters together with the scaling of variables (32), the dimensionless value of critical wavenumber becomes Q c = 17.89. Therefore, the finite geometry shall select a primary mode with wavenumber Q c + Q 0 = 18, i.e. the computational domain shall accommodate D = 9 . Indeed, by choosing a value B E,0 < B < B 0 , starting from a small perturbation of the branch 0 , the primary branch results to be stable even for values of the control parameter smaller than the primary bifurcation threshold, so confirming the subcritical character (see Fig. 4a). Let us now set the initial condition as a small perturbation of the branch 1 , which is characterized by Q c + Q 1 = 17 and D = 8.5 , and choose a control parameter which falls into the range B 1 < B < B E,1 or B E,1 < B < B 2 , i.e. respectively slightly below or above the Eckhaus threshold of the first secondary branch (predicted by quintic GL equation). Results reveal that, in the former case, the system undergoes a phase slip and finally converges towards the stable 0 branch (see Fig. 4b). In the latter case, the system remains on the stable 1 branch, proving the occurrence of a restabilizing bifurcation (see Fig. 4c). The instability of the last secondary branch for values smaller than Eckhaus threshold can be proved by setting the initial condition as a small perturbation of 2 (Q c + Q 2 = 19 and D = 9.5 ) and choosing B 2 < B < B E,2 . As shown in Fig. 4d, the system selects a different wavelength and stabilizes along the stable 0 branch. On the contrary, if the control parameter is chosen in such a way that B > B E,2 , no phase slips occurs and the system stabilizes on the 2 branch, as depicted in Fig. 4e. In all cases, results agree with theoretical predictions developed in Sect. 3.2.1.  2 and e B > B E,2 . The different initial conditions here used are specified in the main text Let us finally comment that all the simulations reported in this section have been obtained by integrating numerically the hyperbolic model by using μ 0 = 10 5 and γ 0 = 1.1, that represent a good approximation of the behavior close to the parabolic limit in the subcritical regime. In the next section, the behavior far away from the parabolic limit is explored in more detail.

Spatio-temporal dependence of phase slip
Let us characterize more deeply the phenomenon of phase slip taking place during Eckhausdriven dynamics, with more emphasis on the role played by hyperbolicity.
As previously mentioned [5,18,41,42], a phase slip may be defined as a solution of the GL equation whose zeros of the pattern amplitude | (x, t)| vary as a function of space and time. Let us call t sli p and x sli p , respectively, the first time instant and the spatial location at which such an event takes place. Due to the impossibility to predict these values theoretically, ad-hoc numerical investigations are performed with the twofold aim of: (i) elucidating how inertial times may affect this phenomenon in the supercritical and subcritical regimes and (ii) establishing strategies to control these quantities independently of each other.
Let us now investigate how the time and the space at which phase slip takes place depend on the model parameters and the initial conditions. In this context it will be emphasized how inertial times represent, effectively, additional degrees of freedom that may be used to enrich the scenario of transient dynamics. Indeed, in the parabolic version of the Klausmeier model [28][29][30][31], the linearized problem giving rise to Turing threshold is completely determined by all model parameters, i.e. A, B and d, as deduced in (4), (5). These parameters also specify the stationary amplitude of each branch of quantized periodic state. Consequently, no free parameters are available to characterize the evolution of the system during the transient  (42) it is possible to manage two extra parameters, i.e. γ 0 and μ 0 , that are strictly correlated to the inertial times of the involved species and that are always present in any physical system. As shown in the previous sections, changes applied to inertial times do not yield any variation of existence and stability thresholds, but it is expected that they play a relevant role during transient regime.
To this aim, let us integrate numerically the hyperbolic system by fixing the parameters A, B, d and γ 0 , and varying the phenomenological coefficient μ 0 . In the first set of investigations, the parameters are chosen in such a way dynamics fall into the supercritical regime. In particular, the plant loss is chosen in the range B e,2 < B < B E,2 and the initial condition is set as small perturbation of the unstable branch 2 , so that the system is expected to undergo an Eckhaus-driven transition towards the stable branch 1 (as depicted in Fig. 2d). Results are shown in Fig. 5 where μ 0 is varied from 10 2 to 10 5 . It can be noticed that, as the parameter μ 0 increases, the phase slip takes place earlier, leaving its spatial location at the edges of the domain unaltered. Moreover, dynamics reported from (e) to (a) represent progressive deviations from the parabolic limit and allow to appreciate how much longer transient regimes may be experienced by varying a single inertial-related parameter.
To describe in more detail the control of the only time to phase slip, let us characterize its functional dependence. From the numerical standpoint, keeping in mind (6) 2 and (A.13), an useful way to approximate the occurrence of a zero of pattern amplitude consists of searching where and when the current state U(x, t) crosses the steady-state U S . To check the validity of the previous statement, let us depict the trajectory followed by the field variables from the initial state U in (taken as small perturbation of the unstable 2 branch) towards the stationary patterned state U end (stable patterned branch 0 ). As illustrative examples, the trajectories corresponding to three different values of μ 0 = 20, 10 2 and 10 5 are shown in Fig. 6a for the supercritical regime, whereas those corresponding to μ 0 = 10 2 , 10 3 and 10 5 for the subcritical one are depicted in Fig. 7a. It is found that, independently of the value of μ 0 , all the trajectories cross the steady homogeneous state U S at the same location: x sli p = 84 in Fig. 6 a, b Numerically-computed trajectories u(x sli p , t) corresponding to supercritical dynamics from a small perturbation of the unstable 2 branch (U in ) towards the stable patterned 0 branch (U end ), crossing the homogeneous state U S . Fixed parameters: B e,2 < B < B E,2 , d = 10  The above observations suggest to correlate the time to phase slip t sli p with the linear growth rate λ (1) n (see expressions (26) and (37)). Indeed, once the growth factor is recasted in the original time variables, say λ (1) n , one argues that it depends linearly on the parameters ν (in the supercritical regime) or ν (in the subcritical one), which both bring the contribution from inertia. Indeed, by plotting the dependence of t sli p as a function of λ (1) n , with n = 2, two different linear relationships are retrieved in log-log scale, as depicted in Fig. 6c for the supercritical regime (with adjusted r-squared value r 2 a = 0.999), and in Fig. 7c for the subcritical one (r 2 a = 0.995). These results confirm that the time to phase slip t sli p , close to onset, is strictly related to the linear growth rate λ (1) n in both dynamical regimes. Notice that the above analysis is carried out by considering values of inertial times falling either in the ranges τ u < τ w and τ w < τ u , obtaining the same functional dependence. However, the different signs of the proportionality coefficients deduced in the linear fits indicate that the time to phase slip is negatively correlated with the linear growth rate in the supercritical regime [40] Results described so far pointed out that the time at which phase slip occurs may be controlled by the linear growth rate, leaving the location unchanged. Let us now inspect whether it is possible to find a strategy to get the opposite scenario, namely the possibility to modulate the location x sli p at which wavelength adjustment takes place leaving the time to phase slip unaltered. To this aim, since perfectly-periodic initial conditions may represent an idealization of the reality, let us consider an ecologically-realistic scenario where a defect, localized at x ≡ x de f with amplitude a, is present in the initial data, as shown schematically in Fig. 8a. Let us then build the scatterplots representative of location x sli p to phase slip as a function of defect position x de f , obtained for a = 50. Results shown in Fig. 8b reveal that the spatial location of phase slip can be modulated linearly (r 2 a = 0.999) by varying the defect position, leaving the time to phase slip almost unchanged. To confirm the above statements, in Fig. 9 some examples of spatiotemporal evolution of patterns obtained for three different locations of a large defect are depicted. As it can be noticed, the time to phase slip is kept almost fixed at t 350.
Therefore, the numerical investigations here proposed suggest different strategies to control, separately, the time and the spatial location at which phase slip takes place.

Conclusions
In this paper, a theoretical study on the EI of stationary patterns in a hyperbolic reactiondiffusion system defined over a large finite domain was addressed in the supercritical and subcritical regimes. To this aim, linear and weakly nonlinear analysis were first addressed to deduce the equations governing the pattern amplitude close to criticality. Then, the bifurcation analysis of the EI was carried out to describe existence and stability thresholds of all bifurcating branches. Finally, the above analytical predictions were validated by comparison with results of numerical simulations performed on a conceptual model of interest in drylands ecology, the modified Klausmeier model, where stationary vegetation patterns emerge over flat arid terrains. These numerical investigations also allowed to get some insights into the phenomenon of phase-slip.
The main results can be summarized as follows: (i) It was shown that, while the cubic real GL equation (12) describes satisfactorily well the features associated to primary and secondary branches of periodic solutions in the supercritical regime, it is not able to predict the proper range of existence of such periodic solutions as well as the restabilizing mechanism that n-th order secondary branches, with n ≥ 1, undergo in the subcritical regime. The above restrictions were overcome by pushing weakly nonlinear analysis up to the fifth order, where the resulting envelope equation takes the form of cubic-quintic GL equation (15). To achieve this goal, it was necessary to deduce the explicit expressions of the main properties of primary and secondary quantized periodic branches (stationary amplitude n , existence ζ e,n and stability ζ E,n thresholds, linear growth rate λ (1) n ) characterizing the bifurcation diagram in the subcritical regime in the framework of quintic GL equation. All these quantities were collected in Tables 1-3. This approach provided the complete description of all (n ≥ 0) periodic branches appearing in the bifurcation diagram of subcritical modes, as can be noticed from the comparative theoretical-numerical results reported in Tables 4 and 5.
(ii) From the inspection of numerically-computed bifurcation diagrams (see Figs. 1, 3) and spatio-temporal evolution of patterns (see Figs. 2, 4), it was obtained the twofold goal of: validating the theoretical predictions arising from multiple-scale weakly-nonlinear analysis and elucidating the role played by inertial effects that are always present in any physical system. From the theoretical viewpoint, in particular, it was shown that hyperbolicity affects the expression of the linear growth rate but leaves the other quantities (stationary amplitude, existence and stability thresholds) unchanged. This result suggested that the hyperbolic model provides additional degrees of freedom that may be used to better characterize transient regimes, in particular the transition from an Eckhaus-unstable state towards a more favorable stable configuration through sequences of phase-slips.
(iii) The functional dependence of time and location at which wavelength adjustment takes place is investigated numerically with the main goal of finding strategies to control these quantities independently of each other. Our results revealed that the time to phase slip may be modulated by varying the inertial times leaving its location unaltered. Moreover, it was shown that t sli p strongly depends upon the linear growth rate λ (1) n , in both dynamical regimes (as depicted in Figs. 6 and 7). In particular, results indicate that t sli p and λ (1) n are negatively correlated in the supercritical regime but positively in the subcritical one. On the other hand, it was numerically shown that, by considering a localized defect into the initial conditions, it is possible to modulate the location to phase slip leaving the time to phase slip almost unchanged (see Figs. 8 and 9).
In future works, it is planned to extend bifurcation analysis of EI to the case of oscillatory periodic patterns in the more general framework of hyperbolic reaction-advection-diffusion models, where the instability thresholds are remarkably affected by inertial effects and pattern amplitude close to onset is ruled by complex GL equations. GNFM through "Progetto Giovani GNFM 2020" entitled "Analisi di biforcazione e teoremi di buona posizione in modelli matematici multi-scala di interesse biologico".
Data Availability All data generated or analysed during this study are included in this published article.

Conflict of interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.