On the stability of relativistic perfect ﬂuids with linear equations of state p = K (cid:2) where 1 / 3 < K < 1

For 1 / 3 < K < 1, we consider the stability of two distinct families of spatially homogeneous solutions to the relativistic Euler equations with a linear equation of state p = K ρ on exponentially expanding FLRW spacetimes. The two families are distinguished by one being spatially isotropic while the other is not. We establish the future stability of nonlinear perturbations of the non-isotropic family for the full range of parameter values 1 / 3 < K < 1, which improves a previous stability result established by the second author that required K to lie in the restricted range ( 1 / 3 , 1 / 2 ) . As a ﬁrst step towards understanding the behaviour of nonlinear perturbations of the isotropic family, we construct numerical solutions to the relativistic Euler equations under a T 2 -symmetry assumption. These solutions are generated from initial data at a ﬁxed time that is chosen to be suitably close to the initial data of an isotropic solution. Our numerical results reveal that, for the full parameter range 1 / 3 < K < 1, the density gradient ∂ x ρρ associated to a nonlinear perturbation of an isotropic solution develops steep gradients near a ﬁnite number of spatial points where it becomes unbounded at future timelike inﬁnity. This behavior of the density gradient was anticipated by Ren-dall (Ann Henri Poincaré 5(6):1041–1064, 2004), and our numerical results conﬁrm his expectations.


Introduction
Relativistic perfect fluids with a linear equation of state on a prescribed spacetime (M, g) are governed by the relativistic Euler equations 1   ∇i T ij = 0 ( where T ij = (ρ + p)ṽ i ṽj + pg ij is the stress energy tensor, ṽi is the fluid four-velocity normalized by gij ṽi vj = −1, and the fluid's proper energy density ρ and pressure p are related by p = Kρ.
Since K = dp dρ is the square of the sound speed, we will always assume 2 that 0 ≤ K ≤ 1 in order to ensure that the speed of sound is less than or equal to the speed of light.We further restrict our attention to exponentially expanding Friedmann-Lemaître-Robertson-Walker (FLRW) spacetimes (M, g) where M = (0, 1] × T 3 and 3 3) It is important to note that, due to our conventions, the future is located in the direction of decreasing t and future timelike infinity is located at t = 0. Consequently, we require that ṽ0 < 0 holds in order to guarantee that the four-velocity is future directed.For use below, we find it convenient to introduce the conformal four-velocity via (1.4) 1 Our indexing conventions are as follows: lower case Latin letters, e.g.i, j, k, will label spacetime coordinate indices that run from 0 to 3 while upper case Latin letters, e.g.I, J, K, will label spatial coordinate indices that run from 1 to 3.
2 While this restriction on the sound speed is often taken for granted, it is, strictly speaking, not necessary; see [7] for an extended discussion.
3 By introducing a change of time coordinate via t = − ln(t), the metric (1.2) can be brought into the more recognizable form g = −d t ⊗ d t + e 2 tδ ij dx I ⊗ dx J , where now t ∈ [0, ∞).
1.1.Stability for 0 ≤ K ≤ 1/3.It can be verified by a straightforward calculation that (ρ * , v i * ) = (t 3(1+K) ρ c , −δ i 0 ), t ∈ (0, 1], (1.5) defines a spatially homogeneous solution of the relativistic Euler equations (1.1) on the exponentially expanding FLRW spacetimes (M, g) for any choice of the parameter 0 ≤ K ≤ 1 and constant ρ c ∈ (0, ∞).From a cosmological perspective, these solutions are, in a sense, the most natural since they are also spatially isotropic and hence do not determine a preferred direction.The future, nonlinear stability of the solutions (1.5) on the exponentially expanding FLRW spacetimes was first established in the articles 4 articles [20,22] for the parameter values 0 < K < 1/3.Stability results for the end points K = 1/3 and K = 0 were established later 5 in [13] and [8], respectively.See also [6,10,11,14] for different proofs and perspectives, the articles [9,12] for related stability results for fluids with nonlinear equations of state on the exponentially expanding FLRW spacetimes, the articles [4,23,26] for analogous stability results on other classes of expanding cosmological spacetimes, and [19] for related, early stability results for the Einstein-scalar field system.One of the important aspects of these works is they demonstrate that spacetime expansion can suppress shock formation in fluids, which was first discovered in the Newtonian cosmological setting [25].This is in stark contrast to fluids on Minkowski space where arbitrary small perturbations of a class of homogeneous solutions to the relativistic Euler equations form shocks in finite time [3].
A consequence of these stability proofs is that the spatial components of the conformal four-velocity v i of small, nonlinear perturbations of the homogeneous solution (1.5) decay to zero at future timelike infinity, that is, lim for the parameter values 0 ≤ K < 1/3.This behaviour is, of course, consistent with the isotropic homogeneous solutions (1.5).On the other hand, when K = 1/3, the spatial components of the conformal four-velocity v i for perturbed solutions do not, in general, decay to zero at timelike infinity, and instead limit to a spatial vector field ξ I on T 3 , that is, lim This behaviour is consistent with a family of non-isotropic homogeneous solutions defined by 6 where (ρ c , ν c ) ∈ (0, ∞) × (0, ∞), which satisfy the relativistic Euler equations for K = 1/3.The known stability results for K = 1/3 imply the future stability of nonlinear perturbations of these solutions.
Noting that solutions of the type (1.6) can be made arbitrarily close to solutions of the type (1.5) for K = 1/3 by choosing ν c sufficiently small, from a stability point of view there seems to be verify little difference between the two classes of solutions for small ν c .Indeed, the future nonlinear stability of both classes of solutions, where ν c is sufficiently small, can be achieved via a common proof.However, as will become clear, the essential difference between these solutions is that, from an initial data point of view, stable perturbations of solutions of the type (1.5) are generated from initial (ρ, v I )| t=1 that is sufficiently close to (ρ * , v I * )| t=1 and satisfies min x∈T 3 (g IJ v I v J ) t=1 = 0, (1.7) while stable perturbations of solutions of the type (1.6) (1.8) 1.2.Stability for 1/3 < K < 1.Until recently, it was not known if any solutions of the relativistic Euler equations were stable to the future for the parameter values 1/3 < K < 1.In fact, it was widely believed that solutions to the relativistic Euler were not stable for these parameter values.This belief was due, in part, to the influential work of Rendall [18] who used formal expansion to investigate the asymptotic behaviour of relativistic fluids on exponentially expanding FLRW spacetimes with a linear equation of state.Rendall observed that the formal expansions become inconsistent for K in the range 1/3 < K < 1 if the leading order term in the expansion of g IJ v I v J at t = 0 vanished somewhere.He speculated that the inconsistent behaviour is the expansions could be due inhomogeneous features developing in the fluid density that would ultimately result in the blow-up of the 4 In these articles, stability was established in the more difficult case where the fluid is coupled to Einstein's equations.However, the techniques used there also work in the simpler setting considered in this article where gravitational effects are neglected. 5Again, stability was established in these articles in the more difficult case where the fluid is coupled to Einstein's equations. 6More generally, we could set the spatial components of the conformal four-velocity v I • to be any non-zero vector in R 3 and determine v 0 • via the conditions g ij v i • v j • = −1 and v 0 • < 0. However, for simplicity, we will assume here that v I • is chosen so that it is pointing in the direction of the coordinate vector field density contrast ∂ I ρ ρ at future timelike infinity.Speck [23, §1.2.3] added further support to Rendall's arguments by presenting a heuristic analysis that suggested uninhibited growth would set in for solutions of the relativistic Euler equations for the parameter values 1/3 < K < 1 .Combined, these considerations left the stability of solutions to the relativistic Euler equations in doubt for K in the range 1/3 < K < 1.
However, it was established in [16] that there exists a class of non-isotropic homogeneous solutions of the relativistic Euler equations that are stable to the future under small, nonlinear perturbations.This class of homogeneous solutions should be viewed as the natural continuation of the solutions (1.6) over the parameter range 1/3 < K < 1, and they are defined by where and u = u(t) solves the initial value problem (IVP) Existence of solutions to this IVP is guaranteed by Proposition 3.1.from [16], which we restate here: for all t ∈ (0, 1].Moreover, for each ρ c ∈ (0, ∞), the solution u determines a homogeneous solution of the relativistic Euler (1.1) equations via (1.4) and (1.9).
The main result of [16] was a proof of the nonlinear stability to the future of the homogeneous solutions (1.9) for the parameter values 1/3 < K < 1/2.It was also established in [16] that under a T 2 -symmetry assumption, future stability held for the full parameter range 1/3 < K < 1.An important point that is worth emphasising is that the initial data used to generate the perturbed solutions from [16] satisfies the condition (1.8) at t = 1, and furthermore, this positivity property propagates to the future in the sense that the perturbed solutions satisfy It is this property of the perturbed solutions from [16] that avoids the problematic scenario identified by Rendall.This article has two main aims: the first is to establish the nonlinear stability to the future of the homogeneous solutions (1.9) for the full parameter range 1/3 < K < 1 without the T 2 -symmetry that was required in [16].The second aim is to provide convincing numerical evidence that shows the density contrast blow-up scenario of Rendall is realized if the condition (1.8) on the initial data is violated.
Before stating a precise version of our stability result for the homogeneous solutions (1.9), we first recall two formulations of the relativistic Euler equations from [16].The first formulation, which was introduced in [15] and subsequently employed in [14] to establish stability for the parameter range 0 < K ≤ 1/3, involves representing the fluid in terms of the modified fluid density ζ defined via and the spatial components v I of the conformal fluid four-covelocity7 v i = g ij v j .In terms of these variables, the relativistic Euler equations (1.1) can be formulated as the following symmetric hyperbolic system: where ) and (1.24) The second formulation of the relativistic Euler equations is obtained by introducing a new density variable and decomposing the spatial components of the conformal fluid four-velocity as where u(t) solves the IVP (1.11)-(1.12).Then setting ψ =t 2µ + e 2 w1 , (1.30) where and Π = diag(0, 0, 1, 1). (1.41) For later use, we also define and observe that Π and Π ⊥ satisfy the relations An important point regarding the formulation (1.36) is that it is symmetrizable.Indeed, as shown in [16], multiplying (1.36) by the positive definite, symmetric matrix where it is straightforward to verify from (1.37)-(1.39) that the matrices are symmetric, that is, (A I ) tr = A I . (1.47) We are now in a position to state the main stability theorem of this article.The proof is presented in Section 2. Before stating the theorem, it is important to note that, due to change of variables defined via (1.25)-(1.28)and (1.35), the homogeneous solutions (1.9) correspond to the trivial solution W = 0 of (1.36).
) is the unique solution to the IVP (1.11)-(1.12)from Proposition 1.1 and ζ0 , w 0 J ∈ H k+1 (T 3 ).Then for δ > 0 small enough, there exists a unique solution to the initial value problem Moreover, (i) W = ( ζ, w J ) tr satisfies the energy estimate for all t ∈ (0, 1] where8 ) such that the estimate Ē(t) t µ−σ holds for all t ∈ (0, 1] where (iii) u and W = ( ζ, w J ) tr determine a unique solution of the relativistic Euler equations (1.1) on the spacetime region M = (0, 1] × T 3 via the formulas (iv) and the density contrast (1.55) 1.3.Instability for 1/3 < K < 1.It is essential for the stability result stated in Theorem 1.2 to hold that the initial data used to generated the nonlinear perturbations of homogeneous solutions of the type (1.9) satisfies the condition (1.8).This leaves the question of what happens when this condition is violated, which would be guaranteed to happen for some choice of initial data from any given open set of initial data that contains initial data corresponding to an isotropic homogeneous solution (1.5).To investigate this situation, we consider a T 2 -symmetric reduction of the system (1.15) obtained by the ansatz where ζ is as defined above by (1.25).It is not difficult to verify via a straightforward calculation that the relativistic Euler equations (1.15) will be satisfied provided that z and w solve9 (1.59) In Section 3, we numerically solve this system for specific choices of initial data Importantly, these choices include initial data for which w 0 crosses zero at two points in T 1 , and as a consequence, violates the condition (1.8).From our numerical solutions, we observe the following behaviour: (1) For all K ∈ (1/3, 1) and all choices of initial data (z 0 , w 0 ) that are sufficiently close to homogeneous initial data of either family of solutions (1.5) and (1.9), z and w remain bounded and converge pointwise as t 0.
(2) For each K ∈ (1/3, 1) and each choice of initial data (z 0 , w 0 ) that violates (1.8) and is sufficiently close to homogeneous initial data of the family of solutions (1.5), there exists a = (K) ∈ Z ≥0 such that sup This indicates an instability in the H -spaces for solutions of (1.58)-(1.59) that is not present, c.f. Theorem 1.2, in solutions generated from initial data satisfying (1.8).We also observe that the integer is a monotonically decreasing function of K with a minimum value of 1.For the initial data we tested, the blow-up at t = 0 in the derivatives occurs at a finite set of spatial points.
(3) For all K ∈ (1/3, 1) and all choices of initial data (z 0 , w 0 ) that are sufficiently close to homogeneous initial data of either family of solutions (1.5) and (1.9), solutions to (1.58)-(1.59)are approximated remarkably well, for times sufficiently close to zero, by solutions to the asymptotic system10 ∂ t z = 0, (1.60) everywhere except, possibly, at a finite set of points where steep gradients form in z, which only happens for K large enough and initial data violating (1.8).
(4) For each K ∈ (1/3, 1) and each choice of initial data (z 0 , w 0 ) that violates (1.8) and is sufficiently close to homogeneous initial data of the family of solutions (1.5), the density contrast ∂xρ ρ develops steep gradients near a finite number of spatial points where it becomes unbounded as t 0. This behaviour was anticipated by Rendall in [18], and it is not consistent with either the standard picture for inflation in cosmology where the density contrast remains bounded as t 0, or with the behaviour of the density contrast of solutions generated from initial data satisfying (1.8), c.f. Theorem 1.2.
1.4.Stability/instability for K = 1.When the sound speed is equal to the speed of light, i.e.K = 1, it is well known that the irrotational relativistic Euler equations coincide, under a change of variables, with the linear wave equation.Even though the future global existence of solutions to linear wave equations on exponentially expanding FLRW spacetimes can be inferred from standard existence results for linear wave equations, a corresponding future global existence result for the irrotational relativistic Euler equations does not automatically follow.This is because the change of variables needed to interpret a wave solution as a solution of the relativistic Euler equations requires the gradient of the wave solution to be timelike.Thus an instability in the irrotational relativistic Euler equations can still occur for K = 1 if the gradient of the wave solution starts out timelike but becomes spacelike somewhere in finite time.This phenomena was shown in [5] to occur in the more difficult case where coupling to Einstein's equations with a positive cosmological constant was taken into account.In fact, it was shown in [5] that all wave solutions generated from initial data sets that correspond to a sufficiently small perturbation of the FLRW fluid solution (i.e.(1.5) in our setting) become spacelike in finite time.This proves that the self-gravitating version of the isotropic homogeneous (1.5) are unstable, and in the irrotational setting at least, characterizes the cause of the instability.What is not known is if the other family of homogeneous solutions (1.9) or their self-gravitating versions remain stable for K = 1.
1.5.Future directions.The most natural and physically relevant generalization of the the stability result stated in Theorem 1.2 would be an analogous stability result for the coupled Einstein-Euler equations with a positive cosmological constant for K satisfying 1/3 < K < 1.We expect that establishing this type of stability result is feasible by adapting the arguments from [14].This expectation is due to the behaviour of the term t −2 ρv i v j , which is the only potentially problematic term that could, if it grew too quickly as t 0, prevent the use of the arguments from [14].However, by Theorem 1.2, we know that ).This shows that t −2 ρv i v j decays quickly enough as t 0 to expect that it should not be problematic.We are currently working on generalizing Theorem 1.2 to include coupling to Einstein's equations with a positive cosmological constant, and we will report on any progress in this direction in a follow-up article.We are also planning to investigate numerically, under a Gowdy symmetry assumption, if a similar behaviour, as described in Section 3, occurs for initial data that violates (1.8) when coupling to Einstein's equations with a positive cosmological constant is taken into account.Multiplying this equation through by t µ gives Applying Π ⊥ to (1.36), we further observe, with the help of (1.43), that Next, we decompose the term Π ⊥ A I ∂ I W in (2.2) as follows Inserting this into (2.2) and multiplying the resulting equation on the left by 3) It is worth noting at this point that it is the use of this equation to control Π ⊥ W instead of (2.2) that is responsible for the improvement of the range of the parameter values for which stability holds from 1/3 < K < 1/2 in [16] to 1/3 < K < 1 in this article.Now, multiplying (2.3) by and adding the resulting equation to (2.1) yields Setting W := Π ⊥ W + t µ ΠW = ( ζ, w 1 , t µ w 2 , t µ w 3 ) tr , (2.5) it is then not difficult to verify that the above equation can be expressed as where ) a short calculation using (1.41)-(1.42),(1.44), and (2.4) shows that the matrices (2.7) are given by and

.10)
From these formulas, it is clear that the matrices B i are symmetric, that is, We proceed by differentiating (1.36) spatially to get ) we can write this as Multiplying the above equation on the left by A 0 and recalling the definitions (1.46), we find that WJ satisfies Finally, combining (2.6) and (2.13) yields the Fuchsian system where ) and As will be established in Step 2 below, the Fuchsian system (2.14) satisfies assumption needed to apply the Fuchsian global existence theory from [2]; see, in particular, [2, Thm.3.8.]and [2, §3.4.].This global existence theory will be used in Step 3 of the proof to establish uniform bounds on solutions to the initial value problem (1.48)-(1.49)under a suitable small initial data assumption.These bounds in conjunction with a continuation principle will then yield the existence solutions to (1.48)-(1.49) on (0, 1] × T 3 as well as decay estimates as t 0.
(2.36) Additionally, by (1.41)-(1.43),we observe that P satisfies P 2 = P, P tr = P, ∂ t P = 0 and ∂ I P = 0, (2.37) while the symmetry of the matrices A i , that is,  2.31) of the matrices A 0 , S, A 0 , B 0 , A I , B I , A I and the source term G, the lower bound (2.27) on B 0 , and the identity (2.32), it is not difficult to verify that for each µ ∈ (0, ∞) that there exists a constant θ > 0 such that  16 ).It is also clear that we can view (2.14) as an equation for the variables W = ( W , W J ), with W = ( ζ, w 1 , w2 , w3 ) and W J = ( ζJ , w 1J , w2J , w3J ), where the maps A i and F depend on the variables (t, W ) and (t, W ), respectively.
Taken together, (i) the variable definitions (1.29), (2.5), (2.12), (2.20) and (2.30), (ii) the smoothness properties (2.21), (2.26), (2.28), (2.29) and (2.31) of the matrices A 0 , S, A 0 , B 0 , A I , B I , A I and the source term G, (iii) the identities (2.34)-(2.35)and the lower bound (2.23) satisfied by matrix A 0 , (iv) the definitions (2.17) and (2.19) of the matrices A I and the source term F , (v) the properties (2.37) of the projection map P, and (vi) the bounds (2.39) and (2.40) on ∂ t A 0 and F , respectively, imply that for any11 µ ∈ (0, ∞) and R > 0 chosen sufficiently small, there exist constants θ, γ 1 = γ1 , γ 2 = γ2 > 0 such that the Fuchsian system (2.14) satisfies satisfies all the assumptions from Section 3.4 of [2] for following choice of constants: under the time transformation12 t → t p , the transformed version of (2.14) will satisfy all of the assumptions from Section 3.1 of [2].Moreover, since the matrices A I have a regular limit as t 0, the constants b and b from Theorem 3.8 of [2] vanish.This fact together with β 1 = 0 and κ = κ = µ implies that the constant 13 z from Theorem 3.8 of [2] that determines the decay rate is given by z = µ.

2.3.
Step 3: Existence and uniqueness.By (1.44) and (1.47), we know that the matrices A i are symmetric.Furthermore, from the analysis carried out in Step 2 above, we know that the matrices A i and the source term A 0 G depend smoothly on the variables (t, w J ) for t ∈ (0, 1] and w J in an open neighbourhood of zero, and that the matrix A 0 is positive definite on this neighbourhood.As a consequence, the system (1.48) is symmetrizable and can be put in the symmetric hyperbolic form (1.45) by multiplying it on the left by the matrix A 0 .Since k ∈ Z >3/2+1 and W 0 := ( ζ0 , w 0 J ) tr ∈ H k+1 (T 3 , R 4 ), we obtain from an application of standard local-intime existence and uniqueness theorems and the continuation principle for symmetric hyperbolic systems, see Propositions 1.4, 1.5 and 2.1 from [24,Ch. 16], the existence of a unique solution to IVP (1.48)-(1.49)where T * ∈ [0, 1) is the maximal time of existence.From the computations carried out in Step 1 of the proof, this solution determines via (2.5) and (2.12) a solution where we observe that (2.44) On the other hand, by Step 2 we can apply 14 Theorem 3.8.from [2] to the time transformed version of (2.14) as described in [2, Section 3.4] to deduce, for δ > 0 chosen small enough and the initial data satisfying (2.44), difficulties as we can change between these two conventions by using the simple time transformation t → −t. 13 This constant is denoted by ζ in the article [2].We denote it here by z because ζ is already being used denote the modified density. 14It it is important to note that the regularity k ∈ Z >3/2+1 of the initial data (2.44) is less than what is required to apply Theorem 3.8.[2] to the Fuchsian system (2.14).The reason that we can still apply this theorem is that the matrices A I in (2.14) do not have any 1/t singular terms; see Remark A.3.(ii) from [1].
the existence of a unique solution to the IVP (2.42)-(2.43)with T * = 0 that satisfies the following properties: (1) The limit of P ⊥ W * as t 0, denoted P ⊥ W * (0), exists in H k−1 (T 3 , R 16 ).(2) The solution W * is bounded by the energy estimate for all t ∈ (0, 1], where the implied constant depends on δ. (3) For any given σ > 0, the solution W * satisfies the decay estimate for all t ∈ (0, 1], where the implied constants depend on δ and σ.By uniqueness, the two solutions W and W * to the IVP (2.42)-(2.43)must coincide on their common domain of definition, and consequently, W (t) = W * (t) for all t ∈ (T * , 1].But this implies by (2.41), the energy estimate (2.45), and Sobolev's inequality [17,Thm. 6 By shrinking δ if necessary, we can, by (2.44), make W 0 H k as small as we like, which in turn, implies via the above estimate that we can bound W by W (t) W 1,∞ ≤ R 2 for all t ∈ (T * , 1], where R > 0 is as determined in Step 2 of the proof.This bound is sufficient to guarantee that the matrices A i and the source term A 0 G from the symmetric hyperbolic system (1.45) remain well defined and that the matrix A 0 continues to be positive definite.By the continuation principle and the maximality of T * , we deduce that T * = 0, and hence that W (t) = W * (t) for all t ∈ (0, 1].From this and the energy estimate (2.45), it then follows with the help of the definitions (1.41)-(1.42),(2.5), (2.12), (2.18) and (2.41) that where We further obtain from the decay estimate (2.46) and (2.36) the existence of functions ζ * , w * 1 ∈ H k−1 (T 3 ) and w * 2 , w * 3 ∈ H k (T 3 ) such that Ē(t) t µ−σ (2.47) for all t ∈ (0, 1] where We also note by (1.4), (1.14), (1.17To complete the proof, we find from differentiating (1.50) that the density contrast can be expressed as . (2.48) Since µ > 0, we can choose σ > 0 small enough so that µ − σ > 0. Doing so then implies by (2.47) that ζ and w 1 converge in H k−1 (T H k−2 = 0, which completes the proof.

Numerical solutions
3.1.Numerical setup.In the numerical setup that we use to solve the system (1.58)-(1.59), the computational domain is [0, 2π] with periodic boundary condition, the variables z and w are discretised in space using 2 nd order central finite differences, and time integration is performed using a standard 2 nd order Runge-Kutta method (Heun's Method ).As a consequence, our code is second order accurate 15 . 15Strictly speaking one also needs to enforce the CFL condition to ensure convergence.In this case we have used the tightened 4/3 CFL condition for Heun's Method which is discussed in [21].
3.1.1.Convergence tests.We have verified the second order accuracy of our code with convergence tests involving perturbations of both types of homogeneous solutions (1.5) and (1.9).In our convergence tests, we have evolved the system (1.58)-(1.59)staring from the the two initial data sets (z 0 , w 0 ) = (0, 0.1 sin(x)) ) We have estimated z| t=0 , w| t=0 , w H | t=0 by taking the values of the functions at a time-step close to t = 0 and calculated w H (t) using second order central finite differences.As shown in Figure 5, the numerical solutions clearly replicate the above decay rates suggesting the code is correctly implemented.3.2.Numerical behaviour.Beyond the convergence tests, we have generated numerical solutions to the system (1.58)-(1.59)from a variety of initial data sets (z 0 , w 0 ) for which w 0 satisfies the conditions (1.7) and (1.8).We employed resolutions ranging from 1000 to 160,000 grid points in our simulations.For initial data satisfying (1.7), we chose functions w 0 that cross the x-axis at least twice, 16 while for initial data satisfying (1.8), w 0 does not cross the x-axis at all.
All of the solutions in this article displayed in the figures are generated from initial data of the form (z 0 , w 0 ) = (0, a sin(x + θ) + c) for some particular choice of the constants a, c, θ ∈ R. From our numerical solutions, we observe, for the full parameter range 1/3 < K < 1 and all choices of the initial data with a sufficiently small, that z and w remain bounded and converge pointwise as t 0; see Figures 6 and 7.
3.2.1.Derivative blow-up at t = 0.While z and w remain bounded, our numerical simulations reveal that derivatives of the solutions of sufficiently high order blow-up at t = 0 for the parameter values 1/3 < K < 1 and initial data satisfying (1.7).In Table  3.2.2.Asymptotic behaviour and approximations.For the full range of parameter 1/3 < K < 1 and all choices of initial data, we observe that our numerical solutions display ODE-like behaviour near t = 0.In particular, these solutions can be approximated by solutions of the asymptotic system (1.60)-(1.61) at late times using the following procedure: (i) Generate a numerical solution (z, w) of (1.58)-(1.59)from initial data (z 0 , w 0 ) specified at time t 0 > 0.
(ii) Fix a time t0 ∈ (0, t 0 ) when the numerical solution (z, w) first appears to be dominated by ODE behaviour.(v) Compare the numerical solution (z, w) to the asymptotic solution (z, w) on the region (0, t0 ) × T 1 .
Using this procedure, we find that numerical solutions (z, w) of the system (1.60)-(1.61)can be remarkably well-approximated by solutions (z, w) of the asymptotic system.In particular, by setting t = 0 in (3.9) and noting that we can solve for w| t=0 to get wf := w| t=0 = sgn(w 0 )|w 0 | where sgn(x) is the sign function, we have, with the help of (3.8), that (z, w)| t=0 ≈ (z 0 , wf ). (3.11) It is worth noting that this ODE-like asymptotic behaviour of solutions generated from initial data satisfying (1.8) is expected by Theorem 1.2.What is interesting is that this behaviour of solutions persists for initial data that violates (1.8).
To illustrate how well solutions (z, w) of (1.58)-(1.59)can be approximated by solutions (z, w) of the asymptotic system (1.60)-(1.61)near t = 0, we compare in Figure 8 the plot of wf = w| t=0 , for a fixed choice of t0 and w0 (see (3.10)), with that of w(t) at times close to zero.From the figure, it is clear that the agreement is almost perfect for times close enough to zero. 1−K e (1+K)z where ρ c ∈ (0, ∞).Differentiating this expression, we find after a short calculation that the density contrast is given by .12)Using this formula to compute the density contrast for numerical solutions of (1.58)-(1.59),we observe from our numerical solutions that density contrast displays markedly different behaviour depending on whether or not it is generated from initial data satisfying (3.2).For solutions generated from initial data satisfying (3.2), we find that the density contrast remains bounded and converges as t 0 to a fixed function, which is expected by Theorem 1.2.An example of this behaviour is provided in Figure 9. On the other hand, the density contrast of solutions generated from initial data violating (3.2) develop steep gradients and blows-up at t = 0 at isolated spatial points; see Figure 10 for an example of this behaviour.
As in Section 3.2.2,we can compare the density contrast of the full numerical solutions with the density contrast computed from a solutions of the asymptotic equation.We do this by evaluating (3.12) at t = 0 and using (3.11) to approximate the density contrast at t = 0 by This formula identifies, at least heuristically, that the blow-up at t = 0 in the density density contrast is due the vanishing of w.Once again the agreement between the numerical and asymptotic plots is close enough that the two are practically indistinguishable as can be seen from Figure 11.Plots of density contrast, ∂xρ ρ , calculated from numerical results (Blue) and the asymptotic map (Green).K=0.45, N = 160000, (z 0 , w 0 ) = (0, 0.1 sin(x)).Points near w 0 = 0 in the asymptotic map have been removed to emphasise agreement of the plots away from the singularities.
, w 0 ) = (0, 0.1 sin(x) + 0.15) (3.2) using resolutions of N = 200, 400, 800, 1600, 3200, and 6400 grid points.The initial data (3.1) and (3.2) satisfy the conditions (1.7) and (1.8), respectively, and the solutions generated from this initial data represent perturbations of the homogeneous solutions (1.5) and (1.9), respectively.To estimate the error, we took the base 2 log of the absolute value of the difference between each simulation and the highest resolution run.The results for are shown in Figures1, 2, 3 and, 4 from which the second order convergence is clear.

Table 1 .
1, we list, for a selection of K values, the corresponding minimum value of for which supx∈T 1 |∂ x z(t, x)| + |∂ x w(t, x)| ∞ as t 0.From these values, it appears that is a monotonically decreasing function of K. Observed value of for various K