On the speed of propagation in Turing patterns for reaction-diffusion systems

This study investigates transient wave dynamics in Turing pattern formation, focusing on waves emerging from localised disturbances. While the traditional focus of diffusion-driven instability has primarily centred on stationary solutions, considerable attention has also been directed towards understanding spatio-temporal behaviours, particularly the propagation of patterning from localised disturbances. We analyse these waves of patterning using both the well-established marginal stability criterion and weakly nonlinear analysis with envelope equations. Both methods provide estimates for the wave speed but the latter method, in addition, approximates the wave profile and amplitude. We then compare these two approaches analytically near a bifurcation point and reveal that the marginal stability criterion yields exactly the same estimate for the wave speed as the weakly nonlinear analysis. Furthermore, we evaluate these estimates against numerical results for Schnakenberg and CDIMA (chlorine dioxide-iodine-malonic acid) kinetics. In particular, our study emphasises the importance of the characteristic speed of pattern propagation, determined by diffusion dynamics and a complex relation with the reaction kinetics in Turing systems. This speed serves as a vital parameter for comparison with experimental observations, akin to observed pattern length scales. Furthermore, more generally, our findings provide systematic methodologies for analysing transient wave properties in Turing systems, generating insight into the dynamic evolution of pattern formation.


Introduction
While patterns induced by self-organisation, for example by Turing's mechanism of diffusion-driven instability, are traditionally considered in terms of stationary solutions (Murray, 2003), there has also been extensive interest in spatio-temporal behaviours arising from such systems, for instance the propagation of patterning from a localised disturbance (Tarumi and Mueller, 1989;Myerscough and Murray, 1992;Liu, Maini and Baker, 2022).As a specific example, we consider the common example of Schnakenberg kinetics (Schnakenberg, 1979) on a one-dimensional spatial domain: Here a, b, D 1 and D 2 are positive parameters and the boundary conditions are zero flux.This system has one spatially homogeneous steady state (2) 1 arXiv:2403.09247v1[nlin.PS] 14 Mar 2024 For a localised perturbation about this steady state, a stationary Turing pattern eventually forms via a transient pattern propagating across the domain at an approximately constant speed after initiation, as highlighted in Figure 1.Here we focus on this transient wave and its properties, especially its asymptotic speed, amplitude and profile.We also consider a modelling representation, albeit simplified, of the CDIMA (chlorine dioxide-iodine-malonic acid) reaction kinetics taken from Konow, Somberg, Chavez, Epstein and Dolnik (2019), which is based on the two-variable version of the kinetics (Lengyel and Epstein, 1991), where the short range species u 1 inhibits the production of u 2 , with the latter promoting the production of u 1 .This constitutes a well known example of pure kinetics Murray (2003), in distinct contrast to the cross kinetics of Schnakenberg, with the opposite interactions, whereby u 1 catalyses the production of u 2 , which inhibits the production of u 1 .The speed of such a wave is amenable to the application of the marginal stability criterion of Dee and Langer (1983); Tarumi and Mueller (1989) and Myerscough and Murray (1992).Originally formulated as a hypothesis (Dee and Langer, 1983), the application of the marginal stability criterion is subject to various heuristics to motivate its use.Marginal stability has been motivated as corresponding to the transition between convective and absolute instabilities (Rovinsky and Menzinger, 1992;Tobias, Proctor and Knobloch, 1998;Sandstede and Scheel, 2000;Sherratt, Dagbovie and Hilker, 2014;Ponedel, Kao and Knobloch, 2017) while it is argued that the application of marginal stability relies on the observation that there is a velocity of the moving reference frame at which the leading edge of a perturbation neither grows nor decays, as can be assessed by a Fourier representation of the solution in a fixed reference frame (Ben-Jacob, Brand, Dee, Kramer and Langer, 1985).The exact conditions under which this criterion 2 can be applied are still an open question; furthermore, marginal stability provides no information on the amplitude or profile of the transient wave.Thus we will also consider such a propagating front in terms of weakly non-linear analysis and envelope equations (Hoyle, 2006), which is justified by weakly non-linear analysis and multiple scale asymptotics, and additionally can provide estimates for the wave shape and amplitude.However, we acknowledge at the outset that both techniques are unlikely to be valid if a subcritical bifurcation is present; for instance marginal stability is considered to be limited to pulled fronts (Ponedel et al., 2017;Avery, Holzer and Scheel, 2023), while weakly non-linear analysis is not sufficient to resolve complexities such as the impact of multiple steady states influencing the dynamics.With a focus on supercritical instabilities, the primary aim of this paper is to provide systematic methodologies for determining the properties of propagating patterning transients for Turing systems, such as those observed in Figure 1b.Noting that it is simpler to implement, the marginal stability approach is presented first.Then we develop a weakly non-linear analysis to determine the envelope equation, which takes the form of a real Ginzburg-Landau equation (GLE) to determine estimates, not only for the wave speed, but also for the wave amplitude and profile, which we compare with numerical simulations.In addition, we also compare wave speed predictions and analytical expressions from both methods, with the objective of ascertaining whether the simpler, more tractable, marginal stability method can provide similar accuracy to that of envelope methods for predicting the speed of patterning transients, and whether envelope methods can provide a quantitative justification for the use of marginal stability wave-speed estimates.
2 Analytics -the foundations for determining the characteristic speed of patterning We consider a one-dimensional reaction-diffusion (RD) system of two species, subject to homogeneous Neumann (zero flux) boundary conditions, where α denotes a bifurcation parameter such that α = 0 at a bifurcation point corresponding to the loss of stability of the homogeneous steady state (HSS), u * (α), i.e.R(u * (α), α) = 0. Without loss of generality, the choice of sign of α entails the HSS is stable for α < 0 and unstable for α > 0. We shall denote vectors as u and matrices as A throughout the text.

Marginal stability
The marginal stability approach (Dee and Langer, 1983;Ben-Jacob et al., 1985;Van Saarloos, 1988;Tarumi and Mueller, 1989;Myerscough and Murray, 1992;Chomaz and Couairon, 2000) identifies the location in time and space where the instability arises in the form of a travelling wave due to a localised perturbation.We define g(k where k is the wavenumber and σ(k) is the classical dispersion relation (Murray, 2003;Krause, Gaffney, Maini and Klika, 2021;Klika, 2017).The marginal stability criterion then states that 0 = dg dk , ℜg = 0, evaluated at the leading edge.However, the diverse arguments supporting the formulation of marginal stability conditions are all heuristic (Dee and Langer, 1983;Ben-Jacob et al., 1985;Van Saarloos, 1988;Tarumi and Mueller, 1989;Myerscough and Murray, 1992;Chomaz and Couairon, 2000;Van Saarloos, 2003) and it is hypothesised that there is a change in the separation of variables of the solution from travelling wave (TW) coordinates to pattern-generating envelope fronts with fixed periodicity in space (Dee and van Saarloos, 1988).
For the complex derivative dg dk to exist, we assume g to be analytic in k, so that the Cauchy-Riemann conditions hold and the marginal stability conditions can be rewritten as where Note that c decouples from this set of equations (c = 1 k I ℜσ) and we are left with two equations for the two unknowns k R , k I .
In practical terms, equations (4) are solved as follows.First, the dispersion relation σ(k) (recall that both σ, k ∈ C) is obtained as a solution to the following quadratic equation with complex coefficients where J denotes the linearised kinetics, that is J ij = R i,uj and we denote partial derivatives via the comma notation, so that R i,uj = ∂Ri ∂uj u * .Consequently, the following set of two equations (quadratic in k), implicitly defining the real and imaginary parts of σ(k), is obtained: This pair of equations, together with the first two equations in ( 4), where we rewrite all the powers of k in terms of its real, k R , and imaginary, k I , parts constitute a set of four relations for four unknowns.Finally, the travelling wave speed is given by c = 1 k I ℜσ while the wavenumber of the pattern can be estimated as k MS = ℑg/c = k R + ℑσ/c.The latter follows from the observation that the imaginary part of the shifted dispersion relation to the moving frame of the travelling wave is the frequency of the pattern scaled by the wave velocity c.

Marginal stability close to bifurcation point
In order to have a better understanding of the marginal stability approach, we asymptotically solve the marginal stability conditions close to the bifurcation point.To this end, we rescale the bifurcation parameter as α = ϵ 2 , expand the linearised kinetics as J = J 0 + ϵ 2 J 1 + ϵ 4 J 2 , where, for example, . Standard linear analysis reveals that the bifurcation point is characterised by the condition and the critical wavenumber k c is given by the condition for a repeated root for k 2 (Murray, 2003) As is usually the case with asymptotic solutions, the appropriate scaling can only be determined by calculation and the existence of a dominant balance in all considered orders.When considering regular asymptotic expansions in the small parameter ϵ, this leads to (by trial and error) the choice of where κ R , κ I play the role of the unknown variables while we denote ℜσ = ℜσ 0 + ϵ 2 ℜσ 1 + O(ϵ 4 ) and The nonlinear leading order problem is satisfied exactly at the bifurcation point, i.e. when and, as a result, the leading order TW speed is zero.
The first subleading problem is solved sequentially as follows.First, the implicit equations for the real and imaginary parts of σ, eqns.( 5) and ( 6), yield and ℑσ 1 = 0.
From the imaginary part of the first condition for marginal stability, ∂ℑσ ∂k R = −c, we obtain the following relation for speed c (with the expected scaling with the distance from the bifurcation point): Further, the second condition, ℜg = 0, gives the relation for κ I : Therefore, the marginal stability conditions close to the bifurcation point give the following estimate of the TW speed (in t, x dimensional coordinates) To obtain the correction to the wavenumber and to check that we indeed have a plausible asymptotic solution with the assumed scaling, we calculate κ R from the real part of the first condition in (4) for marginal stability noting that the Cauchy-Riemann conditions are satisfied at both considered orders.

Evaluation of marginal stability criterion performance
We consider Schnakenberg and CDIMA kinetics (as two representatives of kinetics which are in phase, CDIMA, and out of phase, Schnakenberg) to compare the marginal stability results with numerical solutions both near and far from the bifurcation point.To this end, we consider a primed Turing system, that is, we choose parameters from the Turing space (diffusion-driven instability region).Both models involve two morphogens.The Schnakenberg kinetics are specified above in Eq. ( 1) while the CDIMA kinetics are where a, b, µ are (positive) model parameters.
For the case of zero flux boundary conditions, the spatially homogeneous stationary solution of the Schnakenberg model has been given in Eq. (2) while, for the CDIMA kinetics, we have We now consider the following set of parameters in the Schnakenberg kinetics above: D 1 = 1, D 2 = 20, a = 0.05, b = 1.4 and a 1D domain with size L = 200 (Schnakenberg I).This set of parameters is within the Turing space as indicated in Figure 2, where we plot the Turing space in the a, b parameter space for D 1 , D 2 fixed at 1, 20 respectively.We also highlight the localisation of the nearest bifurcation point b c = 1.712 for the chosen bifurcation parameter α = b c − b, which leads to ϵ = √ α = 0.559.
We use Wolfram Mathematica 12.0 to solve the full PDE system with zero flux boundary conditions and a local perturbation of the homogeneous steady state u * in the second component u 2 using the method of lines and 10,000 points for spatial discretisation.Pattern formation is initiated from a local perturbation well inside the domain, at x = 60, and travels in both directions at an apparently fixed speed, see Figure 1b, in agreement with the previous observations in the literature (Liu et al., 2022).
The choice of magnitude for the local perturbation, as characterised by ρ (see Fig 1a), impacts the duration of the transient time before the travelling wave velocity asymptotes to its final speed, as may be observed in Fig 3 for the Schnakenberg system in that the deviation between the travelling wave speed and its large ρ asymptote in the plot is due to the duration of the transient.Since we solve the PDE numerically on a bounded domain and for different kinetics and various parameter values, the transient time will inevitably vary.Thus, we opt for an initial condition of significant magnitude, as demonstrated in Fig. 3, such that the front velocity has adequately asymptoted within the considered The horizontal axis (in log scale) expresses ρ > 0, the magnitude of the perturbation in the second component as described in Fig. 1a, while on the vertical axis we plot the numerically calculated c RD using the same methodology as described in the text for the Schnakenberg I system.To this end, we use the two black vertical lines in Fig. 1b highlighting the locations of two cross-sections for the assessment of the pattern amplitude A RD , wavenumber k RD and speed of the travelling wave c RD .Note that the intermediate transition around ρ = 0.2 is due to finite size effects, when slight variations in ρ cause an abrupt insertion of an additional half-wave of a pattern into the domain, that is, a slight variation of the pattern wave number k RD occurs.
bounded domain.In accordance with this objective, henceforth we select an initial perturbation for all numerical computations to be three times the homogeneous stationary solution in the second variable, that is ρ = 3, with the perturbation of the first component remaining at zero.Then the numerically determined velocity can provide a reasonable estimation of the asymptotic rate of pattern propagating in the environment.Such an estimation can then be used for comparison.
From the numerical solution determined in this manner, we may extract all the observed characteristics from the (numerical) localisation of the highlighted maxima of the pattern, as in Figure 1b, with: i) the amplitude A RD = 0.982; ii) the pattern wavenumber k RD = 0.571, denoting 2π divided by the distance between the neighbouring maxima; and iii) the speed of the travelling wave c RD = 1.288 from numerically calculating the times when the travelling pattern has reached half of the final amplitude at the two highlighted locations.
We now compare the results of the numerical solutions to those of the marginal stability analysis.To this end, we denote the predictions of the wavenumber and travelling wave speed from marginal stability by k MS and c MS , respectively.Numerical solution of Eqns.(4) for the above case gives k MS = 0.555 and c MS = 1.375 for Schnakenberg I, for example.
We repeat this process for other choices of parameter values and we also consider the CDIMA kinetics.The results are summarised in Table 1, where absolute and relative errors with respect to the numerical solution of both wave number and front speed are shown (more details for the CDIMA calculations are given in Appendix A).Note that the longest transient time, and hence the greatest numerical error in calculated properties k RD and c RD , can be expected near the bifurcation point, when the perturbation growth rate is very slow and hence the requirement for a sufficiently large region for determining the velocity of front propagation is greatest.Hence, the larger relative error near the bifurcation point is not necessarily a sign of a poorer performance of the marginal stability criterion.Table 1: A summary of marginal stability results (index MS) and numerical estimation (index RD) for the key properties of the TW and resultant pattern for CDIMA and Schnakenberg kinetics.We compare each case using relative and absolute error in both the front propagation speed and the wavenumber of the pattern.The domain is one-dimensional with length L = 200 and D 1 = 1, with b c defined as in Fig 2 and Fig 7a.
Next, we compare the TW velocity prediction over a broader parameter space around the bifurcation point for Schnakenberg kinetics.We use the idea of continuity of solutions of the algebraic marginal stability equations on the bifurcation parameter to estimate the localisation of their solutions for nearby parameter values and then compute the exact solution.This allows calculation autonomy, for example, in the initial estimates of the roots of the algebraic problem (and there might be more than one solution).However, the key for comparison is the estimation of the TW speed from the numerical solution of the full RD system, c RD .When the bifurcation parameter is varied, there is a change in the characteristic time of pattern formation as well as the localisation of the two maxima used to determine the travelling wave velocity.Again, to automate the process, we store the locations of the maxima as indicated by the two black vertical lines in the density plots, Fig. 1b, and use the continuous dependence of the solution on the parameter to find their localisation after a small variation of the bifurcation parameter in their neighbourhood.Finally, we use 10k points for spatial discretisation (further refinement does not lead to more than 1% difference in the values of the final pattern).
The results of this approach are shown in Fig 4 with plots for the TW velocity in x, t variables.It can be seen that we indeed observe a behaviour corresponding to the square root of the distance from the bifurcation point in its neighbourhood and that there is a match between the asymptotic expression   for the travelling wave speed close to the bifurcation point and the numerical solution of the marginal stability criterion.Finally, note that the equations for marginal stability may have multiple solutions, as can be seen in Fig 4, where multiple branches of c MS appear (even close to the bifurcation point).

Envelope equation, GLE
We now consider the amplitude equation approach.We perturb the HSS v = u − u * and expand to identify the approximate evolution equation of the perturbation: where the linear operator L I = D ∂ 2 ∂x 2 + J is the linearised RD problem with J ij = R i,uj evaluated at u * , which is the case for higher derivatives below as well.As we assume smooth kinetics R, L II and L III are symmetric and we have denoted the nonlinear terms in the following way hence L II represents the quadratic and L III the cubic kinetics approximation.Now, recalling the bifurcation parameter in general is α = ϵ 2 , the regular expansion of the perturbation near the bifurcation point α = 0 Taking advantage of the known relation in scaling of time and space near a bifurcation point (Hoyle, 2006) to balance two spatial derivatives with the one temporal derivative, we introduce the following set of slow variables: τ = ϵ 2 t, y = ϵx, while we consider the original independent variables t, x as fast.Using the method of multiple scales, writing v j as functions of both slow and fast variables, that is v j (t, τ, x, y), where and by equating powers of ϵ, we get: while we require strict periodicity in the spatial fast variable, corresponding to the microscale oscillations.
As we search for an envelope equation for a travelling wave with a stationary profile in t close to a bifurcation point, as motivated in the Introduction, we look for a solution of the form (Hoyle, 2006) where k c is the critical wavenumber corresponding to the bifurcation (onset of instability at the bifucation point), and V exp(ik c x) is the Turing (Fourier) eigenmode of the linearised problem, i.e.
where the explicit expression for k c is Eq. ( 8).This follows from the requirement that the larger of the two eigenvalues of −J D , where eigenvector of J D associated with the zero eigenvalue.We remark that this assumed form of solution, Eq. ( 14), implicitly defines the microscale and hence the periodic boundary conditions are imposed at (0, 2π/k c ) in the fast spatial variable x.
By looking at the O(ϵ 2 ) problem (13b) while taking into account the form of v 1 , see Eq. ( 14), one can see that v 2 has to be of the form where c.c. denotes complex conjugate.In particular, by comparing the coefficients of eigenvectors (invoking their orthogonality, see below), we obtain noting that the relations for the complex conjugate are exactly the same.From the first and last equations we directly have: while for the unknown Z 1 we obtain the singular equation as (k 2 c D − J 0 ) = J D and the critical wavenumber k c is such that J D has a zero eigenvalue.As we shall see below, after the discussion of the solvability condition for v 3 , the solvability condition for this singular set of equations is equivalent to the definition of the critical wavenumber k c and hence Z 1 is well defined (although not unique).
We now return to the second subleading order problem (13c).As its left-hand side is the same as the leading order problem which has a solution for zero right-hand side, we know from the Fredholm alternative (recalling that we assume v j is independent of the fast variable t) that the right-hand side of Eq. (13c) has to be perpendicular to the solution f of the adjoint homogeneous problem (L I 0 ) * f = 0. From the definition of the adjoint operator, ⟨L I 0 f, g⟩ = ⟨f, (L I 0 ) * g⟩, and recalling the strict periodicity on (0, 2π/k c ), we have and hence (L I 0 ) * g = D ∂ 2 g ∂x 2 + J T 0 g.Therefore, the solvability condition is the requirement of perpendicularity to the solution to (L I 0 ) * g = 0. We claim that such solutions are the left eigenvectors of J D , i.e. g = A(τ, y)W exp(ik c x) + c.c., where, WLOG, W T • V = 1 (considering simple eigenvalues of J).
Choosing V = (1, v) T and W = 1 1+vw (1, w) T , we have W T • V = 1, while for V to be the (right) eigenvector of J D corresponding to the zero eigenvalue we require Similarly, w follows from the fact that W has to be the (right) eigenvector of J D T corresponding to the zero eigenvalue, i.e.
Hence, we may now write the first order correction v 2 (15) explicitly as where we set where z is arbitrary and the solvability condition for Z 1 , eq. ( 16c), reads D 1 +vwD 2 = 0 (being equivalent to the definition of k c ).Note that the arbitrariness of z, and hence non-uniqueness of Z 1 , is manifested by the fact that any multiple of the vector (1, v) T = V can be added to Z 1 .However, this corresponds to the fact that we can add any product V exp(ik c x) to v 2 , i.e. the solution of the equation with zero right-hand side.This part of the solution is already included in the leading order solution, i.e. in the solution of the problem (13a), and hence we can set z = 0 without loss of generality.Finally, the solvability condition for the second subleading problem (13c) then requires (again, the c.c. terms generate the same expression) noting that W T • D • V = 0 from the solvability condition of Z 1 .Therefore the scalar envelope equation is the real Ginzburg-Landau equation where d = 4k 2 c D1D2w (1+vw)(J0)12 and the coefficient of the linear term is while the cubic term is ) .

Travelling wave analysis
First, let us rewrite the amplitude equation ( 21) in a more suitable form for the travelling wave analysis.Note that in our studied case, we know that c 1 > 0 because we are considering a primed Turing system where a travelling wave develops due to a localised perturbation.Specifically, we are concerned with the system with parameter values in Turing space and hence the homogeneous steady state solution corresponding to A = 0 is linearly unstable, i.e. c 1 > 0. Similarly, we naturally expect the existence of a fixed profile of the desired travelling wave, hence the existence of a positive stationary and homogeneous root of the amplitude equation A * > 0. Thus we have c 3 > 0 to allow for such a solution.We now rewrite the amplitude equation ( 21) as where 0 < A * = c1 c3 and we look for a non-negative solution (amplitude).We introduce the wave variable ξ = y/ √ d − cτ , consider a particular direction c > 0, and look for a travelling wave solution with a fixed profile travelling with a speed c = c√ d.Phase plane analysis to show the existence of the travelling wave is then standard.We rewrite the amplitude equation ( 22) as a first order system in ξ where ′ denotes the derivative with respect to the wave variable ξ.The corresponding boundary conditions (A, W ) = (A * , 0) at ξ = −∞, (A, W ) = (0, 0) at ξ = +∞ represent the sought heteroclinic connection between (unstable) A = 0 and (stable) A = A * .
The fixed point (0, 0) is a focus for c > 2 √ c 1 and a spiral for c ∈ (0, 2 √ c 1 ), while (A * , 0) is always a saddle with the unstable manifold in the direction 1, − c 2 1 − 1 + 8c 3 (A * /c) 2 .From these observations we can argue that there is a critical speed, , below which a travelling wave does not occur.However, there is no unique speed c of a travelling wave as for any c > c * one can construct a heteroclinic connection between the two fixed points satisfying the desired boundary conditions.
We now use the well known Kolmogorov (1937) results about the asymptotics of the speed of a compactly supported initial condition entailing a travelling wave.Namely, as the cubic kinetics satisfy the requirement of two zeros (which can be scaled to be 0 and 1, respectively), while being positive between them and with the highest derivative at zero, the observed travelling wave speed matches the critical speed (corresponding to the transition from a spiral to a focus) for sufficiently large times.
The analytic profile of the travelling wave corresponding to the critical wave speed is not available.However, we can readily find the profile of a heteroclinic connection for c = cA = 3 c1 2 = 3A * c3 2 as matching the desired boundary conditions.As c * = cA with δ ≈ 0.057, we can rewrite the equation for the travelling wave profile as and hence expect the profile to be closely represented by the analytic profile A A (ξ), eq. ( 23), for δ = 0.It can be shown that the appropriate form of the asymptotic expansion for A is and collecting powers of δ after substituting this into Eq.( 24) gives A 0 (ξ) = A A (ξ) and As the solution to the first equation is A 1 (ξ) = KA ′ 0 (ξ), which satisfies the required homogeneous Dirichlet boundary conditions at ±∞, the solvability condition for the second equation is Using the expected form of solution A 1 = KA ′ 0 , and invoking the exact form of A 0 (ξ), Eq. ( 23), a direct calculation yields Finally, the profile of the travelling wave has the following analytical approximation Note that this profile corresponds to a heteroclinic connection for a cubic reaction kinetics between an unstable node and a saddle that is unique in the above sense due to Kolmogorov's result.
In short, the concept of a characteristic speed of pattern propagation in a reaction-diffusion system is well founded (thanks to the dominance of the critical wave speed) and corresponds to  17), (18).The travelling wave profile can be approximated by Eq. ( 25).Crucially, it can be shown that the expression for the TW speed obtained from the envelope equation, Eq. ( 26), and from the marginal stability conditions close to the bifurcation point, Eq. ( 10), are exactly equivalent, when we take into account the conditions that hold at the bifurcation point (7) and D 1 + vwD 2 = 0 (being equivalent to the definition of k c ).To see this, let us compare the ratios of the coefficients (J1)12 (J1)11 , (J1)21 (J1)11 and (J1)22 (J1)11 in the squared expression for the TW velocity near a bifurcation point in the marginal stability approach, Eq. ( 10).If we show that they are in the ratio v, w and vw then, from Eq. ( 26), we will be left with proving that the coefficients of (J 1 ) 11 match in both expressions.First, the ratio of the coefficients of (J 1 ) 22 and (J 1 ) 11 is −D 1 /D 2 , which is equivalent to vw, following from the above expression D 1 = −vwD 2 .Next, the ratio of the coefficients of (J 1 ) 12 and (J 1 ) 11 is c D2 via the expression for the critical wavenumber k c .Finally, the ratio of the coefficients of (J 1 ) 21 and (J 1 ) 11 is D1(J0)22−D2(J0)11 , which is equivalent to w = − (J0)12 (J0)22−k 2 c D2 when using the expression for the critical wavenumber k c and the relation w = − D1 D2 v. Therefore, all that remains is to show that the coefficients of the (J 1 ) 11 terms are equal in both expressions (10) and ( 26).This follows from the fact that the bifurcation point enforces the condition (D 2 (J 0 However, in addition to the prediction of the travelling front wave speed and the spatial frequency of the resultant pattern, the amplitude equation approach estimates the amplitude of the final pattern.Hence, in addition to the above numerical verification of the asymptotic version of the marginal stability criterion, we shall assess the estimation of the amplitude of the pattern.

Assessment of envelope equation performance
It is worth mentioning that several key characteristics follow from the amplitude equation itself, as estimated by Eqn. ( 25).Firstly the asymptotic amplitude of the leading order pattern in v is equal to 2ϵ multiplied by the positive stationary solution of the envelope equation, i.e. 2ϵA * , with A * = c 1 /c 3 , noting the associated component of V = (1, v) and the factor of two arising from a complex conjugate, as follows from Eq. ( 14).We may also compare this leading asymptotic amplitude directly with the amplitude of the first component u 1 of the solution to the full problem.Second, the speed of the travelling wave is 2ϵ √ c 1 d and we have an analytical expression for this in terms of the model parameters, as given by the expressions following Eq.( 21).Third, an estimate of the wave number is immediately available, since the weakly nonlinear theory is based on the localisation of the bifurcation point, which is analytically determinable from linear stability theory, see Eq. ( 8).
We again consider Schnakenberg and the CDIMA kinetics with the same parameter values as above when testing the marginal stability criterion, hence we consider a primed Turing system.Note that the physical parameter we vary as the bifurcation parameter, denoted α, is subject to choice and hence can potentially affect the precision of the analytical results.Nonetheless, the methodology of deriving the envelope equation is not affected by this choice, though the expansion of the Jacobian given by J = J 0 + αJ 1 + O(α 2 ) generates expressions that ultimately are dependent on the choice of α.Finally, on comparing the full problem with the envelope equation, note that the initial perturbation is of the same magnitude (rescaled by ϵ) and location for both problems, allowing us to use the envelope equation for comparison.However, one distinction is that the perturbation is about a non-trivial homogeneous steady state for the full problem but about the trivial solution for the envelope equation (simply being the zero amplitude of the pattern corresponding to the homogeneous steady state), though this does not invalidate direct comparison between the two solutions.
As a particular example, we consider Schnakenberg kinetics with zero flux boundary conditions, the parameter values D 1 = 1, D 2 = 20, a = 0.05, b = 1.4 and a 1D domain with size L = 200 (Schnakenberg I), see Appendix A.2 for the other studied cases.Using the above analysis, we obtain, for this choice of reaction kinetics and parameter values, the following form of the envelope equation and that ϵ = 0.559.Hence, we may simply read out the predicted amplitude of the pattern A env = 2ϵA * = 1.015 and using Eq. ( 26) we have an estimate of the travelling wave speed c env = 1.297.The wavenumber k env is equal to the critical wavenumber k c , Eq. ( 8), and has the value k env = k c = 0.628.
In addition, we solve this scalar reaction diffusion equation with the same zero flux boundary conditions and the same initial condition as above in the full problem, that is, a localised perturbation of the same magnitude, ρ = 3, and location x = 60 about the unstable trivial homogeneous steady state.Plotting this solution reveals a good match with the solution to the full problem, see Figures 5a, 5b, where we plot the solutions in their natural coordinates.
In addition, from the numerical solution to the envelope equation, we determine the numerically observed speed of the travelling wave c num env as described above in the two highlighted locations giving c num env = 1.288, showing a good match with the critical wave speed determined analytically from the Kolmogorov asymptotic arguments.Finally, we plot the solutions, see Figures 6a, 6b, to both the full problem and the envelope equation for the two highlighted cross-sections in Fig. 5b corresponding to x = 91.46 and x = 170.92.In this way we can compare the predicted travelling wave profile together with its amplitude and velocity.We also plot the asymptotic estimate of the travelling wave profile, (25), shifted in time by t 0 to match the arrival of the wave at the first cross-section with the numerical solution to the envelope equation.This time shift t 0 is necessary as the determination of the travelling wave profile does not consider initiation and development of the front and thus is free up to a translational shift, which we fix by specifying t 0 .Its value is determined manually by choosing a value which results in the best match (visually overlapping) between the analytical profile, eq (25), and plotted as a dotted curve, together with a plot of the numerical solution to the envelope problem at the first cross-section (dashed).In particular, we present the function (u * ) 1 + 2ϵA(ξ) where A(ξ) is the identified asymptotic approximation for the amplitude in (25) evaluated at a point ξ(t, x) = ϵ √ d (x − c env ϵ(t + t 0 )).Therefore, there is a single fitting parameter, the time shift t 0 .Therefore, the precision of the analytical estimate of the front velocity c env is visually immediate as the difference between the time of arrival of the wave at the second cross-section profile, Fig 5b.
We also show the numerical solution to the full RD problem where there is no need to provide a time shift correction.Thus, from the figure we are able to determine the velocities observed in all    25).This shift in the plot of the analytical profile is necessary because it was obtained from phase space in TW coordinates and therefore is subject to the translational invariance of the travelling wave, see text for more details.As one can observe from the comparison of the two panels, all the three speeds c RD , c env , c num env of the travelling wave are similar and the approximate analytical profile matches the numerically calculated ones, see the text for more details.Finally, note that there is a disparity in the predicted amplitude of the pattern, the error being 3% of the numerically calculated amplitude from the full problem.three approaches: (i) the full reaction diffusion system with numerically estimate speed, c RD ; (ii) the analytical estimate for the asymptotic travelling wave speed from the envelope equation c env ; and (iii) the numerically estimated asymptotic travelling wave speed from the numerical solution of the envelope equation, c num env .In addition, we can also assess the closeness of fit for the front profiles and the time taken for the pattern to develop from the small localised disturbance.

Discussion and conclusions
The marginal stability criterion was put forward as a conjecture forty years ago (Dee and Langer, 1983) and while there are several hypotheses for its derivation (Dee and Langer, 1983;Ben-Jacob et al., 1985;Van Saarloos, 1988;Tarumi and Mueller, 1989;Myerscough and Murray, 1992;Chomaz and Couairon, 2000;Van Saarloos, 2003), for example, being a transition between convective and absolute instabilities (Dee and Langer, 1983;Rovinsky and Menzinger, 1992;Tobias et al., 1998;Sandstede and Scheel, 2000;Sherratt et al., 2014;Ponedel et al., 2017), its validity still remains an open problem.In addition, the marginal stability equations may not entail a unique speed of the travelling front, see Fig 4, and do not provide information on the profile or amplitude of the wave.
As an alternative route, we derived the envelope equation for the spatial pattern deposited after the front but it is restricted in validity to be close to the primary bifurcation point.However, we were able to show that there is a uniquely selected TW speed via Kolmogorov's argument of dynamical selection via a linear mechanism, also known as the pulled case.Note, however, that the selection of TW speed along the lines discussed above, when the TW speed is determined only by what is happening at the front, is limited only to the pulled front case, where its dynamics are driven by linear kinetics around the unstable homogeneous steady state, see Chomaz and Couairon (2000); Van Saarloos (1988); Berestycki and Hamel (2007).Pushed fronts require different approaches and will be the subject of future research.
The amplitude equation approach offers several benefits.It not only estimates the front velocity and pattern wavenumber, but also the magnitude of the deposited pattern amplitude and front shape.We derived explicit analytical expressions for all of these properties from the original model parameters.Moreover, the asymptotic solution to the marginal stability criterion yields exactly the same estimate near the bifurcation point.Our analysis establishes a firm basis for the notion of the characteristic speed of pattern propagation in a reaction-diffusion system, while strengthening the motivation for using the marginal stability criterion and provides a footing for the further exploration of a wave of competence (Liu et al., 2022).
In particular, we have shown that the marginal stability criterion and the amplitude equation give exactly the same estimates of the travelling wave speed near primary bifurcation points.Hence, Figure 4 can be also viewed as a comparison of the marginal stability criterion and the amplitude equation approach.This exact match is surprising given the very distinct foundations of the two frameworks but it serves as a confirmation of the marginal stability criterion at least in the vicinity of primary bifurcation points.
Note that other choices of the bifurcation parameter α may yield different predictions as the measure of the distance to the bifurcation point varies.We leave a more thorough discussion of this effect for future research.Nevertheless, we remark that, for example in the CDIMA kinetics case, the choice α = (a − a c )/a c , where a c = 11.058 is the critical parameter value, yields ϵ = 0.292 and hence different predictions of the characteristic features of the travelling wave solutions.Similarly, we also consider the choice α = (µ c − µ)/µ c , where µ c = 63.813 is the bifurcation point.The associated predictions are listed below in Table 2 and further details including the envelope equations are given in Appendix A.
We now summarise the results in the above considered scenarios in Table 2, restricting the results to relative errors for brevity.Note the role of the choice of the bifurcation parameter while keeping the parameter values the same, lines 1-3 in the table.The linear marginal stability criterion showed a good estimation of the TW speed and wavenumber in the studied examples.Similar results were then obtained for these two properties using a weakly nonlinear analysis and the envelope equation, which we showed to exactly satisfy the marginal stability criterion close to the bifurcation point.Furthermore, while one might expect a loss of accuracy of the amplitude equation approach with increasing distance from the bifurcation point, the marginal stability condition is not restricted to the neighbourhood of the bifurcation point and shows good estimates even outside this neighbourhood.However, as we have shown, there is the problem of multiple possible solutions, among which the correct solution corresponding to the characteristic TW speed cannot be selected without further analysis or insight.The envelope equation approach predicted that the evolution of the travelling wave in the envelope equation is slightly delayed compared to the actual full reaction-diffusion problem, the amplitude approach allows us to capture, qualitatively, the shape of the solution and, in many cases, is also reasonably accurate quantitatively, see Table 2.
This study demonstrates the existence of a characteristic speed of pattern propagation determined from detailed reaction kinetics together with diffusion and thus is a hallmark of each Turing pattern formation system.Thus, in turn, each Turing system has an associated wave propagation speed that may be compared to experimental results in the same way as the length scale characteristic of the observed pattern (Míguez, Dolnik, Munuzuri and Kramer, 2006;Konow et al., 2019;Konow, Dolnik and Epstein, 2021;Glover, Wells, Matthäus, Painter, Ho, Riddell, Johansson, Ford, Jahoda, Klika et al., 2017;Glover, Sudderick, Shih, Batho-Samblas, Charlton, Krause, Anderson, Riddell, Balic, Li et al., 2023).Table 2: A summary of marginal stability results (index MS), numerical estimation (index RD) and envelope equation results (index env) of the key properties of the TW and the deposited pattern for D 1 = 1.Note that these models are defined in detail and explored in the main text and in the Appendix.We discuss each case via the relative error in both the front propagation speed and the wavenumber of the deposited pattern.In addition, we show the comparison of pattern amplitude and the role of the bifurcation parameter choice.Finally, we plot the solutions, see Figures 9a, 9b, to both the full problem and the envelope equation for the two highlighted cross-sections in Fig. 7b corresponding to x = 98.83 and x = 165.9.In this way we can compare the predicted travelling wave profile together with its amplitude and velocity.We also plot the asymptotic estimate of the travelling wave profile, (25), shifted in time by t 0 to match the arrival at the first cross-section with the numerical solution to the envelope equation.This time shift t 0 is necessary as the determination of the travelling wave profile does not consider initiation and development of the front and thus is free up to a translational shift, which we fix by specifying t 0 .Its value is determined manually by choosing a value which results in the best match (visually overlapping) between the analytical profile, eq (25) and plotted as a dotted curve, together with a plot of the numerical solution for the envelope problem at the first cross-section (dashed).In particular, we plot the function (u * ) 1 + 2ϵA(ξ) where A(ξ) is the identified asymptotic approximation for the amplitude in (25) evaluated at the point ξ(t, x) = ϵ √ d (x − c env ϵ(t + t 0 )).Hence, there is a single fitting parameter, the time shift t 0 .Therefore, the precision of the analytical estimate of the front velocity c env is visually immediate as the difference between the times of arrival of the waves (full, dashed and dotted curves) at the second cross-section profile, Fig 9b.
We also show a solution to the full RD problem where there is no need to provide a time shift correction.Thus, from the figure we are able to determine the velocities observed in all three approaches: (i) the full reaction diffusion system with numerically estimate speed, c RD ; (ii) the analytical estimate for the asymptotic travelling wave speed from the envelope equation c env ; and (iii) the numerically estimated asymptotic travelling wave speed from the numerical solution of the envelope equation, c num env .In addition, we can also assess the closeness of fit for the front profiles and the time taken for the pattern to develop from the small localised disturbance.
In this situation, from the envelope equation analysis, we obtain and ϵ = 0.346, k env = k c = 0.843, A env = 1.082, c env = 1.575, c num env = 1.477.In Figure 10 we again compare the travelling wave nature, profile, amplitude, and speed as a solution to the full problem and envelope equation, while we also plot the estimated analytical profile of the travelling wave.We again determine the shift t 0 of the analytical amplitude profile as described above in CDIMA I.As a result, we can directly compare the travelling wave velocity and the profile from the three considered approaches, see Figure 10  The envelope equation analysis yields and ϵ = 0.125, k env = 0.915, A env = 0.796, c env = 0.686, c num env = 0.645.We do not plot the solution to this problem but we list the key characteristics in Table 2.

A.3 Other choices of bifurcation parameters
As mentioned in the main text, the choice of the bifurcation parameter may play a role in the weakly nonlinear analysis of behaviour near the bifurcation point.
We remark that the choice α = (a − a c )/a c in CDIMA I model, where a c = 11.058 is the critical parameter value, gives the following results   25).This shift in the plot of the analytical profile is necessary because it was obtained from the phase space in TW coordinates and therefore is subject to the translational invariance of the travelling wave which therefore must be fixed, see text for more details.As one can observe from the comparison of the two panels, all the three speeds c RD , c env , c num env of the travelling wave are similar and the approximate analytical profile matches the numerically calculated ones, see the text for more details.Finally, note that there is a disparity in the predicted amplitude of the pattern, the error being 53% of the numerically calculated amplitude from the full problem.We expect this error to reduce with ϵ.   25).This shift in the plot of the analytical profile is necessary because it was obtained from the phase space in TW coordinates and therefore does not reflect the dynamics of travelling wave formation, see text for more details.Note that all the three speeds c RD , c env , c num env of the travelling wave are similar and the approximate analytical profile matches the numerically calculated ones, see text for more details.Finally, note that there is a disparity in the predicted amplitude of the pattern but it is smaller than in the main text, the error being 13%.and ϵ = 0.465, A env = 1.869, c env = 1.609, c num env = 1.521.In this situation, the wavenumber is estimated as k c = 0.915.

∂A
We do not plot the solutions with these choices but we list the key characteristics for comparison in the first three rows of Table 2.One can see that even though the estimates are carried out for exactly the same system (the same kinetics, diffusion and parameter values), the choice of the bifurcation parameter α affects the accuracy of predictions from the envelope approach.

A.4 Comparison of all methods
As a final illustration, we show the comparison of marginal stability criterion results (and its multiple roots), envelope method, and the characteristics stemming from the numerical solution to the full problem, as an analogue to Fig 4 for the Schnakenberg case in the main text.In Fig. 11, we plot the TW velocity in x, t variables.We again observe a behaviour corresponding to the square root of the distance from the bifurcation point in its neighbourhood and that there is a match between the asymptotic expression for the travelling wave speed close to the bifurcation point and the numerical solution of the marginal stability criterion.Finally, note once more that the equations for marginal stability may have multiple solutions, where multiple branches of c MS appear (even close to the bifurcation point).
A localised perturbation of the homogeneous steady state in the second component u2, which is a piecewise linear function with a support of size 4.It is localised at x = 60 and given in multiples of the homogeneous steady state u * 2 which we express using ρ = u2(t = 0, x = 60)/u * 2 − 1.In this case, we consider ρ = 3.0.(b) We initiate the simulation with a local perturbation in the second component as shown in the left panel.The density plot of the first component of the solution, u1, shows the resulting transient wave of patterning spreading across the domain in both directions (left and right) with an approximately constant velocity, prior to the establishment of a stationary Turing pattern.After an initial transient period, we can indeed observe the asymptotic establishment of a constant travelling wave speed as indicated by the straight dashed line.

Figure 1 :
Figure 1: Illustration of pattern formation behind a wave travelling at a constant speed after a localised perturbation to the steady state.We consider Schnakenberg system, Eqn.(1), with zero-flux boundary conditions and parameter values D 1 = 1, D 2 = 20, a = 0.05, b = 1.4 and a 1D domain with size L = 200 (which we denote as Schnakenberg I).The profile of u 2 is out of phase with the profile of u 1 .

Figure 2 :Figure 3 :
Figure 2: The Turing space, the set of points which allow Turing instability, for Schnakenberg kinetics with parameters a, b, is highlighted in grey, where we fix the parameter values to D 1 = 1, D 2 = 20 and solve on a 1D domain with size L = 200.The full dot represents the chosen point for numerical solution of the full problem (Schnakenberg I).The open circle denotes the closest bifurcation point b c = 1.712 (corresponding to α = 0) for the chosen bifurcation parameter α = b c − b.
a = 12, b = 0.31, Convergence of all the methods to the same square root behaviour close to the bifurcation point bc = 1.712 is apparent.Note that the marginal stability approach shows multiple roots (multiple values of predicted TW speeds) as there are three branches of cMS, two of them existing near the bifurcation point.
A magnification of the figure in the left panel closer to the bifurcation point showing the convergence of the marginal stability asymptotics to the marginal stability criterion results.

Figure 4 :
Figure 4: Travelling wave speed (given in x, t units) as predicted from (i) the full solution of the RD problem with Schnakenberg kinetics, c RD ; (ii) from the marginal stability conditions, c MS ; and (iii) from the marginal stability asymptotics, c MSasympt .Parameter values were chosen as D 1 = 1, D 2 = 20, a = 0.05, and 1D domain with size L = 200 (Schnakenberg I).Note that the choice of the magnitude of the localised initial condition in the finite domain is important, as demonstrated in Fig 3 and we consider ρ = 3.Note that the marginal stability approach can yield multiple solutions, as illustrated here for the whole range of parameter b.Finally, the c RD predictions were only calculated from the bifurcation point b c to b ≈ 0.6 due to finite size effects such as the discretely decreasing number of maxima used to estimate c RD and the longer duration of transient effects before the pattern stabilises.

c
in the original x, t dimensional variables.We recall that d = 4k 2 D2) , see Eqns. ( (a) Density plot of the numerical solution to the envelope equation Eq. (27).The black region is the set of points (y, τ ) where A(τ, y) < Aenv/2 while the white region shows its complement.(b)Density plot of u1 of the numerical solution to the full problem.The black region is the set of points (x, t) where A(t, x) < ARD/2 while the white region shows its complement.

Figure 5 :
Figure 5: Schnakenberg kinetics with D 1 = 1, D 2 = 20, a = 0.05, b = 1.4 and a 1D domain with size L = 200 (Schnakenberg I).The numerical solution for the derived envelope equation is presented in panel (a), and for the full problem, in panel (b).Again, the vertical black lines highlight the positions of cross-sections for determining travelling wave speed and amplitudes.Note the similarity in both the time and location of the pattern initiation and also the development of the spatial pattern behind a front travelling with a fixed velocity (neglecting the initial transition behaviour and boundary effects).We plot the solution to the envelope equation in the corresponding coordinates τ = ϵ 2 t, y = ϵx (with ϵ = 0.559 for the listed parameter values).Hence the two plots are directly comparable and the two thresholded solutions are almost overlapping.

Figure 6 :
Figure 6: Schnakenberg kinetics with D 1 = 1, D 2 = 20, a = 0.05, b = 1.4 and a 1D domain with size L = 200(Schnakenberg I).Temporal profiles for the solutions at the highlighted positions -vertical black lines in the preceding density plots.The full curve corresponds to the numerical solution u 1 (t, x) to the full problem.The dashed curve is the numerical solution to the corresponding envelope problem (u * ) 1 + 2ϵA(ϵ 2 t, ϵx), and the dotted curve (significantly overlapping with the dashed) is the shifted in time estimated analytical profile (u * ) 1 + 2ϵA ϵ √ d (x − c env ϵ(t + 41)) from Eq. (25).This shift in the plot of the analytical profile is necessary because it was obtained from phase space in TW coordinates and therefore is subject to the translational invariance of the travelling wave, see text for more details.As one can observe from the comparison of the two panels, all the three speeds c RD , c env , c num env of the travelling wave are similar and the approximate analytical profile matches the numerically calculated ones, see the text for more details.Finally, note that there is a disparity in the predicted amplitude of the pattern, the error being 3% of the numerically calculated amplitude from the full problem.
(a) The Turing space, the set of points which allow Turing instability, for CDIMA kinetics in a, b parameters, is highlighted in grey, given the other parameter values are D1 = 1, D2 = µ, µ = 50 and a 1D domain with size L = 200.The full dot represents the chosen point in the a, b parameter space for numerical solution of the full problem (CDIMA I).The open circle denotes the closest bifurcation point bc (corresponding to α = 0) for the chosen bifurcation parameter α = bc − b.(b) We initiate the simulation by a local perturbation in the second component as shown in the left panel.The density plot of the first component of the solution, u1, shows the resulting transient wave of patterning spreading across the domain in both directions (left and right) with an approximately constant velocity, prior to the establishment of a stationary Turing pattern.After an initial transient period, we can indeed observe the asymptotic establishment of a constant travelling wave speed as indicated by the straight dashed line.

Figure 7 :
Figure 7: Illustration of pattern formation behind a wave travelling at a constant speed after a localised perturbation to the steady state.We consider the CDIMA system, Eqn.(1), with zero-flux boundary conditions and parameter values D 1 = 1, D 2 = µ, a = 12, b = 0.31, µ = 50 and a 1D domain with size L = 200 (CDIMA I), using an initial perturbation of the form of that presented in Fig 1a with ρ = 3.The profile of u 2 is in phase with the profile of u 1 .

Figure 8 :
Figure 8: CDIMA kinetics with D 1 = 1, D 2 = µ, a = 12, b = 0.31, µ = 50 and a 1D domain with size L = 200 (CDIMA I) using the initial conditions given by Fig 1a with ρ = 3.The numerical solution for the derived envelope equation is presented in panel (a), and for the full problem, panel (b).Again the vertical black lines highlight the positions of cross-sections for determining travelling wave speed and amplitudes.Note the similarity in both the time and location of the pattern initiation and also the development of the spatial pattern behind a front travelling with a fixed velocity (neglecting the initial transient behaviour and boundary effects).We plot the solution to the envelope equation in the corresponding coordinates τ = ϵ 2 t, y = ϵx (with ϵ = 0.293 for the listed parameter values).Hence the two plots are directly comparable and the two solutions are almost overlapping.
and the discussion in the main text.CDIMA III.As a final CDIMA example, we consider parameter values D 1 = 1, D 2 = µ, a = 12, b = 0.38, µ = 50 and L = 200 (CDIMA III).This is a variation of CDIMA I where we have moved the bifurcation parameter b closer to the bifurcation point b c while keeping the remaining parameters fixed.
∂τ = 4.494 ∂ 2 A ∂y 2 + 2.578A − 0.209A 3 .The temporal profile for the numerical solutions and the estimated analytical profile at the first location given by x = 98.83 in Figs.7b and 8b.All solutions are shown in the original t, x variables.The temporal profile for the numerical solutions and the estimated analytical profile at the second location given by x = 165.9 in Figs.7b and 8b.All solutions are shown in the original t, x variables.

Figure 9 :
Figure 9: CDIMA kinetics with D 1 = 1, D 2 = µ, a = 12, b = 0.31, µ = 50 and a 1D domain with size L = 200 (CDIMA I) using the initial conditions given by Fig 1a with ρ = 3. Temporal profiles of the solutions at the highlighted positions -vertical black lines in the preceding density plots, namely in Figs.7band 8b.The full curve corresponds to the solution u 1 (t, x) to the full problem.The dashed curve is the numerical solution to the corresponding envelope problem (u * ) 1 + 2ϵA(ϵ 2 t, ϵx), and the dotted curve (significantly overlapping with the dashed) is the shifted in time estimated analytical profile(u * ) 1 + 2ϵA ϵ √ d (x − c env ϵ(t + 29)) from Eq. (25).This shift in the plot of the analytical profile is necessary because it was obtained from the phase space in TW coordinates and therefore is subject to the translational invariance of the travelling wave which therefore must be fixed, see text for more details.As one can observe from the comparison of the two panels, all the three speeds c RD , c env , c num env of the travelling wave are similar and the approximate analytical profile matches the numerically calculated ones, see the text for more details.Finally, note that there is a disparity in the predicted amplitude of the pattern, the error being 53% of the numerically calculated amplitude from the full problem.We expect this error to reduce with ϵ.
The temporal profile for the numerical solutions and the estimated analytical profile at the first location given by x = 96.81.All solutions are shown in the original t, x variables.The temporal profile for the numerical solutions and the estimated analytical profile at the second location given by x = 170.22.All solutions are shown in the original t, x variables.

Figure 10 :
Figure 10: CDIMA kinetics with D 1 = 1, D 2 = 2µ, a = 10.5, b = 0.4, µ = 13 and a 1D domain with size L = 200 (CDIMA II) using the initial conditions given by Fig 1a with ρ = 3. Temporal profiles for the solutions at the highlighted positions.The full curve corresponds to the solution u 1 (t, x) to the full problem.The dashed curve is the numerical solution to the corresponding envelope problem (u * ) 1 + 2ϵA(ϵ 2 t, ϵx), and the dotted curve (significantly overlapping with the dashed curve) is the shifted in time estimated analytical profile (u * ) 1 + 2ϵA ϵ √ d (x − c env ϵ(t + 31)) from Eq. (25).This shift in the plot of the analytical profile is necessary because it was obtained from the phase space in TW coordinates and therefore does not reflect the dynamics of travelling wave formation, see text for more details.Note that all the three speeds c RD , c env , c num env of the travelling wave are similar and the approximate analytical profile matches the numerically calculated ones, see text for more details.Finally, note that there is a disparity in the predicted amplitude of the pattern but it is smaller than in the main text, the error being 13%.
Convergence of all the methods to the same square root behaviour close to the bifurcation point bc = 0.396 is apparent.Note that the marginal stability approach shows multiple roots (multiple values of predicted TW speeds) as there are two branches of cMS for b < 0.12.
A magnification of the figure in the left panel closer to the bifurcation point showing the convergence of the marginal stability asymptotics to the marginal stability criterion results.

Figure 11 :
Figure 11: Travelling wave speed (given in x, t units)) as predicted from (i) the full solution of RD problem with CDIMA kinetics, c RD ; (ii) from the marginal stability conditions, c MS ; and (iii) from the marginal stability asymptotics, c MSasympt .Parameter values were chosen as D 1 = 1, D 2 = µ, a = 12, µ = 50 and 1D domain with size L = 200 (CDIMA I).Note that the choice of the magnitude of the localised initial condition in the finite domains is important as demonstrated in Fig 3 and we consider ρ = 3.Note that the marginal stability approach can yield multiple solutions, as illustrated here for b ∈ (0.08, 0.12) where two branches appear.