A new instability in clustering dark energy?

In this paper, we study the effective field theory (EFT) of dark energy for the $k$-essence model beyond linear order. Using particle-mesh $N$-body simulations that consistently solve the dark energy evolution on a grid, we find that the next-to-leading order in the EFT expansion, which comprises the terms of the equations of motion that are quadratic in the field variables, gives rise to a new instability in the regime of low speed of sound (high Mach number). We rule out the possibility of a numerical artefact by considering simplified cases in spherically and plane symmetric situations analytically. If the speed of sound vanishes exactly, the non-linear instability makes the evolution singular in finite time, signalling a breakdown of the EFT framework. The case of finite (but small) speed of sound is subtle, and the local singularity could be replaced by some other type of behaviour with strong non-linearities. While an ultraviolet completion may cure the problem in principle, there is no reason why this should be the case in general. As a result, for a large range of the effective speed of sound $c_s$, a linear treatment is not adequate.


Introduction
In the near future, cosmology will benefit from numerous high precision observations [1][2][3], probing the Universe at different epochs, from the very early times when cosmic microwave background radiation (CMB) photons started to propagate and the Universe was around 400,000 years old to today, when the Universe is 13.8 billion years old and has entered a phase of accelerating expansion. One of the main goals of ongoing and future cosmological surveys is to elucidate the physical mechanism behind the late-time accelerated expansion of the Universe that has now been established by several independent observations [4][5][6].
Over the past years, a wide range of theories has been developed by cosmologists and particle physicists with the aim to address the question of the accelerated expansion of the Universe, either by modifying the theory of gravity or by considering an additional fluid component with a negative pressure usually called dark energy (DE) [7][8][9]. Among these theories, the effective field theory (EFT) of dark energy [10][11][12] has become quite popular since it allows to describe the dark energy phenomenology occurring at low energies with a reasonable number of free parameters in a generic way. In principle, the free parameters in the effective theory are connected to (i.e. determined by) a fundamental theory at high energy and respect the low-energy symmetries [13]. In practice however, they can be considered as parameters to be measured in low-energy experiments, similar to the moduli describing an elastic material in the context of material science.
In this paper, we focus on the subset of EFT of DE models with only two free parameters, α K (or equivalently c 2 s = δp/δρ) at the perturbation level and w = p/ρ at the background level, which is equivalent to the well-known theory of k-essence. The k-essence theory was first proposed in 2000 [14][15][16] to naturally, without any fine-tuning, explain the accelerated expansion of the Universe. At linear level k-essence has been explored well and is considered a viable theory for the late-time accelerated expansion. Like the cosmological constant, which is reached in the limit w → −1, this theory can explain all cosmological observations to date. However, by increasing the precision of the observations, we can hope to place much more stringent constraints on the space of theories. In the near future this will require a good understanding of the behaviour at non-linear scales which requires developing proper N -body simulations [17][18][19][20][21][22]. To capture the non-linear behaviour of k-essence we developed the code k-evolution based on the relativistic N -body code gevolution [23][24][25]. In previous studies [26][27][28][29] we maintained linearity of the k-essence field equations and studied the evolution when coupled to a non-linear N -body system. Here, we use the equations derived in [30] for the non-linear evolution of k-essence as an effective field theory, parametrised with the equation of state w and the speed of sound c s . The free parameters appearing in the field picture, e.g., α K , can be interpreted when writing the theory in the fluid picture. In Appendix A of [30], we showed that the fluid description and the field picture are equivalent and one can easily change the picture by well-defined transformations.
In this paper, we show that the EFT of DE for k-essence, in the limit of low speed of sound, suffers from a new instability triggered by one of the non-linear terms in the EFT expansion. In Sec. 2 we discuss the equations that describe the model. In Sec. 3 we present the numerical results for 3+1 D in the cosmological context where we solve the full 3+1 D partial differential equation for the k-essence scalar field numerically, using the EFT framework. We show that for low speed of sound, the numerical solution to this partial differential equation (PDE) blows up at some time before the current age of the Universe. In Sec. 4 we study a simplified PDE, using either planar or spherical symmetry to reduce the dimensionality to 1+1 D. We show analytically that the instability is present and therefore expected to appear in the full 3+1 D case. We also comment on how the solution becomes singular and when the solution ceases to exist, and we show how increasing the speed of sound could stabilise the system. In the final section, we conclude with a short discussion of the results.

Field equations
In this section we write down the equations for the k-essence scalar field parametrised with w and c 2 s , expanded around the background employing the weak-field expansion. The equations of motion as well as the stress energy tensor for clustering DE are obtained and discussed in detail in [30]. There we showed the results for clustering DE where we only keep linear terms in the DE scalar field equations. In order to study the evolution of perturbations we use the Friedmann-Lemaître-Robertson-Walker (FLRW) metric in the conformal Poisson gauge, where Ψ and Φ are the temporal and spatial scalar perturbations of the metric and correspond to the Bardeen potentials, B i is the transverse gravitomagnetic vector perturbation with two degrees of freedom, and h ij is the traceless transverse tensor perturbation with two degrees of freedom. In this gauge the PDE for the k-essence scalar field, including the non-linear corrections in the weak-field regime, reads In this equation π is the DE scalar field and ζ is an auxiliary field written in terms of the scalar field π, its time derivative ∂ τ π with respect to conformal time and the gravitational potential Ψ. Moreover, ( ∇π) 2 ≡ ∇π · ∇π, ∇ is the spatial gradient using partial derivatives and ∇ 2 is the corresponding Laplace operator. As is typical for a weak-field expansion inside the horizon we keep any higher-order terms only if they contain at least two spatial derivatives for each power of a perturbation variable beyond the first order. A simple heuristic argument comes from observing that ( ∇Ψ) 2 /H 2 ∼ v 2 ∼ Ψ and ∇ 2 Φ/H 2 ∼ δ ∼ 1, which implies that each spatial derivative effectively counts as −1/2 order. More details can be found in [30]. As discussed in Appendix A of [30], Eq. (2.3) is equivalent to the continuity and the Euler equations, is the anisotropic stress tensor, δ ij is the Kronecker delta, ∇ i ≡ ∂ i denotes the partial derivative and we have used the following definitions, and where ∇ −2 is the inverse Laplace operator. The effective fluid variables read as follows in terms of the field variables, However, in the implementation chosen for k-evolution, we solve the equations written in the field language and we solve a second-order PDE to update the scalar field π and ζ. Numerical results from k-evolution for the full non-linear PDE show that there exists a critical speed of sound c * s , such that for speed of sound c s smaller than c * s the solution of the PDE becomes singular in finite time. In the following section we show numerical results for examples of small and large speeds of sound where the solution is, respectively, singular and regular, and afterwards we justify the numerical results by studying the equations in simpler setups with spatial symmetries that allow a corresponding reduction of the dimensionality of the problem.

Results from cosmological simulations
In this section we show the numerical results from a first set of cosmological simulations with k-evolution that include non-linear terms in the k-essence field equations. For simplicity we start with standard linear perturbations in matter that are derived from the ΛCDM model and set the two additional fields π and ζ to zero initially. At very high redshift their contribution to the energy density is negligible and hence the matter solution is indeed the one of ΛCDM. The solution of π and ζ will contain a decaying mode due to the way we set the initial conditions, but this will have no relevant effect on the final evolution if our initial redshift was chosen high enough.
Starting the simulation at some initial time and solving the full equations of motion in the weak-field approximation for low speed of sound, one finds numerically that the scalar field π diverges and as a result the simulation breaks down in finite time. In Fig. 1 on the right we show the absolute value of the scalar field in the x − z plane taken at the y-position of the point with the maximal second derivative of π. On the left side we render the same 2D section as a 3D plot, where the height shows the value of the field for better illustration. In these images we see how, over a short period of time in the simulation, an instability is formed around the minimum with largest second derivative and blows up. Since other quantities are coupled to the scalar field, like for example the gravitational potential, they will also diverge and eventually the simulation breaks down. The instability is local and thus if we look at regions far away from the point with maximal second derivative of the scalar field, at the same redshift we see no hint of instability until the simulation itself fails. The reason why the instability is first formed around the minimum with highest curvature will become clear when we study the system in a simplified symmetric setup analytically in Sec. 4. The worked example shown in Fig. 1 is for illustrative purposes and is obtained from a simulation with c 2 s = 10 −7 and w = −0.9. Studying the behaviour for different values of the speed of sound systematically, using k-evolution with the full implementation of clustering DE according to Eqs. (2.2) and (2.3), we find that the solutions are unstable and diverge in finite time for low speed of sound only. Indeed, for otherwise fixed cosmological parameters and w 0 = −0.9, our numerical studies indicate that the system only blows up when c 2 s 10 −4.7 . In Fig. 2 we show the "blow-up redshift" z b for different speeds of sound for a fixed resolution of 0.58 Mpc/h. As the figure suggests, when increasing c 2 s there is a critical value for c 2 s where the system becomes stable. In Appendix A we discuss the effect of precision parameters (temporal and spatial resolution) on our results. We show in Appendix A.1 that increasing the time resolution of the solver does not change the blow-up time significantly. The dependence of our results on the spatial resolution is interesting and can be traced back to the fact that increasing the resolution also enhances the maximum amplitude of perturbations in the initial conditions, discussed further below. This dependence can be understood quantitatively using the extreme value theorem, as we discuss in Appendix A.2.
We study the robustness of the critical c 2 s on the chosen cosmology in Appendix A.3 where we show results from simulations that remain matter dominated forever. We find that the limit for the stability of the equation is not changed significantly even if we let the simulations run far into the future. This rules out the possibility that high speed of sound simulations are only stable due to the limited time available until the gravitational potential starts to decay when matter domination ends and DE takes over.
In  There are two limits, namely high speed of sound for which the system does not blow-up, and very low speed of sound when the system blows up at a redshift close to the blow-up redshift for c 2 s → 0. It is worth mentioning that the blow up redshift depends on the resolution of simulation (see Fig. 8), which also affects somewhat the minimal speed of sound for which the system is stable.
the simulation. As we are going to explain in detail in Appendix A.2, for a fixed box size increasing the resolution of the grid and number of particles would result in increasing the initial density of the perturbations as we are probing smaller scales. But since the blow-up happens first at the point with the highest curvature of the potential, corresponding to the point with the highest density, this in turn affects the blow-up redshift. Our results show that in matter domination there is a linear relation between the blow-up redshift and the initial density. In the next section we will validate this proportionality using a simplified spherically symmetric setup. In Appendix A.2 we discuss the relation between max(δ) and the resolution of a simulation.
Interestingly we also know that the limit of low speed of sound is the limit where important linear terms in the dynamics of the scalar field are suppressed and we end up with a highly non-linear evolution of the field. The critical value c 2 s ∼ 10 −4.7 in Fig. 2 can be understood by comparing the two most important terms in the dynamics of the scalar field. As we are going to show, stability of this system is ensured when the term c 2 s ∇ 2 π dominates over the non-linear term − 1 / 2 H 5c 2 s + 3w − 2 ( ∇π) 2 . Following the spirit of the weak-field expansion mentioned earlier, we can use a dimensional analysis and write ∇ 2 π ∼ L −2 π and ( ∇π) 2 ∼ L −2 π 2 , where L is a characteristic length scale.
Furthermore, in matter domination we have H = 2/τ , and the first-order perturbative solution for π ≈ 1 / 3 Ψτ calculated in Eq. (D.6) in [30]. The two terms will therefore be of similar size if the speed of sound fulfils the relation With w < 0 this relation is satisfied when c 2 s ∼ Ψ, and in cosmology we typically have Ψ ∼ 10 −5 which gives a critical value of c 2 s commensurate with our numerical measurement. Given the fact that the two relevant terms of the equations of motion are independent of the gravitational coupling, the instability is not a result of scalar field self-gravity. This can be checked numerically by turning off the scalar field's contribution to the gravitational field equations (the stars in Figs. 7 and 8 as well as the discussion in Appendix A.2), effectively turning it into a spectator field. Even in this case we find that the perturbations in the scalar field, sourced by the gravitational potentials of dark matter alone, become unstable at almost the same time as in the case with full gravitational coupling. Motivated by these numerical results in the cosmological context we study a simplified version of the full equations analytically in a symmetric setup.

Analytical results
In this section we discuss the most important terms in the PDE governing the scalar field dynamics. In particular, we discuss the non-linear dynamics in 1+1 dimensions if we consider either spherically symmetric or plane-symmetric solutions. We show that the non-linear PDE is suffering from an instability with similar behavior of what we found in the realistic 3+1 dimensional case. Studying the full dynamics even in 1+1 D remains a difficult task and is beyond the scope of this paper. However, the specific term in the equation of motion which we have recognized as the root cause of the instability of the system, namely ( ∇π) 2 [31], is studied thoroughly in [32][33][34][35] in a mathematical context. Here, we study this term with a more physical approach.
We rewrite Eq. (2.3) as a second-order PDE which is more appropriate for analytical studies, where N (π, ∂ τ π, ∇π, ∂ τ ∇π, ∇ i ∇ j π) includes all the non-linear terms, This equation is called a "non-linear damped wave equation" [36] in the mathematics literature. It has been studied by mathematicians in 1+1 D for some types of non-linearities, mainly of the form N (π, ∂ τ π) [37][38][39] but not in the general form appearing in the EFT of DE which is N (π, ∂ τ π, ∇π, ∇ i ∇ j π). In fact, the important remark is that for large speeds of sound and using the fact that π, ∂ τ π and ∇π are small, the dominant term in the dynamics would be the linear part of the PDE, whereas in the limit c 2 s → 0 the term c 2 s ∇ 2 π, which is the restoring force in the linear wave equation, vanishes and therefore the non-linear terms become relevant.
In this limit, i.e. c 2 s → 0, the equation also has a new symmetry whenever the scale factor is a power law in τ and therefore H ∝ τ −1 : it becomes invariant under the rescaling Ψ → λ 2 Ψ, π → λπ, τ → λ −1 τ . The scale invariance also approximately holds at finite c 2 s using c s → λc s as long as c 2 s remains sufficiently small under the rescaling, i.e. for scales much larger than the sound horizon. Scale invariance is indeed expected if the physical problem lacks a characteristic scale such as a sound horizon or a break in the power law in τ . This immediately leads to an interesting conclusion if we consider π as a spectator field in matter domination where it would be sourced by the gravitational potential Ψ produced by nonrelativistic matter. If we assume that Ψ is independent of time (a good approximation when matter is in the linear regime) and have a solution π(τ ; Ψ), we can generate new solutions π(τ ; λ 2 Ψ) = λ −1 π(λτ ; Ψ). Evidently, if for a given gravitational source field Ψ the spectator field π diverges at a certain value z b + 1 ∝ τ −2 b , that value is directly proportional to the overall amplitude of Ψ, i.e. z b + 1 ∝ Ψ. Here we use the fact that z + 1 ∝ τ −2 in matter domination. As we will see shortly, the divergence is actually sensitive to the maximum curvature of Ψ which is of course proportional to Ψ itself in these scaling solutions with constant λ.

Spherical symmetry
In this subsection we study the PDE (4.1) in the regime c 2 s 1 assuming spherical symmetry. We start from Eq. (4.1) and assume a spherically symmetric scenario, i.e. all fields are functions of τ and r only. Furthermore, we choose N = ln a as the new time coordinate, and also rescale π =π/H. The equation reads where we neglect c 2 s against coefficients of order unity (including w) and assume Ψ ∼ Φ in that context. In other words, the only term for which the value of c 2 s is important (once it is assumed that c 2 s 1) is the linear restoring force. Note that if we neglect radiation in the late Universe, the Hubble function is such that Moreover, if we neglect gravitational backreaction from the scalar field π and assume that Ψ is generated by a linear matter perturbation, we can also write the time evolution equation for Ψ, (4.7) We can easily infer that the nonlinear right-hand side is exponentially suppressed at early times as N → −∞ and therefore the initial solution should approach its linear expressioñ π → 2 / 3 Ψ for scales larger than the sound horizon. Note in particular that the linear solution (for w < 0) would be stable at all times if the non-linear self-coupling ofπ is neglected, as was the case for previous numerical studies mentioned in Sec. 1.
Further insight can be found if we consider the solution close to an extremum of the potential. Let us write where the factor H 2 0 is introduced to render the coefficient Ψ 2 dimensionless. We can write the solution forπ as the asymptotic one plus a correction, i.e. π(N, r) = 2 Inserting these ansätze into Eq. (4.7) we find that it neatly separates into two independent parts, one that has no r-dependence, 10) and one that scales as r 2 and reads The second equation is independent of the value of c 2 s . The asymptotic solution for 2 at N → −∞ can be inferred by recognising that 2 is of higher perturbative order than Ψ 2 and therefore becomes subdominant in the nonlinear contribution. One finds that 2 → 8 / 45 e N Ψ 2 2 as N → −∞. At this point it is also worth noting that this asymptotic solution is entirely governed by the first two non-linear terms on the right-hand side of Eq. (4.7), i.e. the terms quadratic in gradients. The other two terms are asymptotically subdominant, and hence do not efficiently trigger the non-linear evolution ofπ. This justifies our claim that the term ( ∇π) 2 is the most relevant non-linear term when discussing the instability. The other terms do, however, have a small effect on the precise time when the divergence of 2 occurs.
In Fig. 4 we show the fully non-linear solution for 2 in the two cases where Ψ 2 = ±1. A divergence only occurs in the case where Ψ 2 is positive (around minima of the potential wells). For Ψ 2 < 0 the solution as N → +∞ approaches the exact solution 2 = − 2 / 3 Ψ 2 − √ −Ψ 2 e −N/2 + 1 / 2 e −N such that the curvature ofπ vanishes asymptotically. However, it decays exactly as fast as H so that the curvature of π =π/H becomes constant.
The scale invariance manifests itself in the fact that Eq. (4.11) is invariant under the transformation Ψ 2 → λ 2 Ψ 2 , 2 → λ 2 2 , N → N − 2 ln λ. Since Ψ 0 does not appear in the nonlinear equations for 0 , 2 , it is clear that the nonlinear instability is indeed governed by the curvature of the potential Ψ, i.e. the amplitude of the coefficient Ψ 2 which is also directly proportional to the matter density contrast. The scaling symmetry implies that the redshift at which the divergence occurs is proportional to this coefficient.
The linear relation between the blow-up redshift and the maximal initial density in matter domination agrees with the results from the three-dimensional simulations of Sec. 3, as shown in Fig. 3, indicating that the analytically tractable spherically symmetric case is able to capture the most important features of the divergence dynamics.
While it might appear that the blow-up occurs independent of the value of c 2 s since the solution of 2 does not depend on it, in a more realistic setting the evolution of the linear solution does matter. The correction due to c 2 s acts to increase/decrease the value at the minimum/maximum ofπ, which under more realistic boundary conditions tends to decrease all derivatives. As a crude estimate we could say that this effect becomes important when the time-dependent term of the linear solution becomes ∼ |Ψ 0 |, which occurs roughly at N ∼ −2 ln c s − ln Ψ 2 + ln |Ψ 0 | + 1.6. The numerical constant is not very important and was computed assuming w −1. We may then argue that the critical speed of sound is given by the estimate N b ∼ −2 ln c s −ln Ψ 2 +ln |Ψ 0 |+1.6, where N b is the value of N at which 2 blows up. From Fig. 4 and using the approximate scale invariance we infer that N b + ln Ψ 2 ∼ 0.9 independent of the value of Ψ 2 . This results in an estimate for the critical speed of sound that is very much in agreement with our previous estimate given in Sec. 3. This argument of course only holds as long as the initial assumption that c 2 This means that considering 0 will not change our estimate of the critical speed of sound in a significant way. On the other hand, due to its coupling to 2 , 0 grows without bounds towards N b . This raises the question whether under more realistic boundary conditions, following the thinking of the previous paragraph, the instability could actually be halted. Related to this question it is interesting to note that the PDE (4.1) has a particular exact solution for c 2 s > 0 if we drop the last two terms on the right-hand side (which are often very subdominant). This can be seen by making the ansatzπ = Ce N + 2 / 3 Ψ, for which the corresponding equation becomes (4.14) A regular solution for which the Hopf-Cole transformation remains defined globally is v(r) = D sinh(kr)/(kr), with k 2 = 5 / 4 H 2 0 (1 − 3w) 2 C/c 4 s and C > 0, D > 0. While this is of course a highly fine-tuned solution, it is interesting to note that the radial profile ofπ does not evolve in this potential. Since we can fix C and D in a way to give any desired second-order Taylor expansion around the minimum of Ψ, this shows that for any c 2 s > 0 there exists a local solution where the gradients ofπ do not evolve and therefore do not lead to a blow-up. Under general initial conditions it remains an open question whether the singularity can be avoided in this way.

Planar symmetry
In this subsection we study the blow-up dynamics in more detail, using an even more simplified version of Eq. (4.1) assuming planar symmetry in a non-cosmological setup, This equation is written following our previous study [31] and the mathematical studies [32][33][34][35] where in addition to the non-linear term we consider the linear term c 2 s ∇ 2 π. We therefore consider the two important terms in the dynamics of the scalar field, i.e. the instability part (∇π) 2 and the pressure term which stabilises the system.
In the limit c 2 s → 0 and rescaling π → π/α the equation reads [31] ∂ 2 τ π = (∇π) 2 , (4. 16) and we consider the initial conditions π(0, x) = 0 , (4.17) First we show that in this case the minima and maxima of the scalar field π do not move in space, a property that we also validate numerically, see Fig. 5. We then derive a PDE for the curvature of the scalar field, which at the extrema satisfies an ODE. We explicitly compute the (finite) time at which the curvature of minima becomes infinite.  It is then evident that for any points x s that are locations of extrema of Ψ, D(τ, x s ) = 0 at all times, i.e. these points are also extrema of π (which remain fixed in position).
Taking a further spatial derivative of Eq. (4.18), we obtain a PDE for the curvature of the scalar field, where κ(τ, x) ≡ ∇D(τ, x) is the curvature. In general this equation is not closed, as we need ∇κ(τ, x) to solve the equation and this term is obtained through a higher-order derivative PDE (by taking another spatial derivative of the equation). However, for the extremal points x s where D(τ, x s ) = 0 we can close the equation as the second term vanishes, This is an ODE for the evolution of the scalar field curvature at the extrema. The initial conditions for κ(τ, x s ) are obtained using the initial profile of π and ∂ τ π in Eq. (4.17), With these initial conditions a first integral of Eq. (4.21) yields where we have dropped the argument x s for brevity. Considering the point x s being a minimum, i.e. ∂ τ κ > 0 and integrating the previous equation from τ = 0 to the blow-up time τ b such that at this time the curvature goes to infinity, we obtain Changing the integration variable from κ to s where s 3 = 4 / 3 κ 3 / α∇ 2 Ψ 2 we find To sum up, the minima blow up in a finite time given by Eq. (4.26) which depends on the initial curvature of the potential Ψ at the minimum. It is worth mentioning that if x s is a maximum, it can also become unstable depending on the initial value of ∂ τ κ. Based on the ODE (4.21), ∂ 2 τ κ is always positive as it is sourced by κ 2 . So we roughly expect that the curvature of maxima starts to increase and eventually becomes flat and switches sign after which it blows up in finite time given by Eq. (4.25). In [35] the dynamics of maxima and minima, especially when c 2 s = 0 is studied in more detail. In Fig. 5 we show the numerical results for the scalar field π(τ, x) profile and its first and second time derivative in the top panel and its spatial derivatives in the bottom panel. The scalar field π(τ, x) and its derivatives are obtained by solving the PDE (4.15) numerically on a lattice in 1+1 D assuming c 2 s = 0, periodic boundary conditions and Ψ(x) = cos(4πx/L) 2 . Bottom: The same scalar field and its spatial derivative on the lattice for different times are shown. Due to the numerical noises appearing in the derivatives of the field we only show the results for the derivatives up to τ = 64. It is also interesting to see that ∂ x π behaves similar to a gradient catastrophe that one would see in some situations in fluid dynamics.
Here N grid = 2048 is the number of points and we choose the units such that α = 1 and dx = 1 is the distance between the points on the 1D lattice. Hence L = N grid = 2048 and as a result ∇ 2 Ψ(x s ) = 3.765 × 10 −5 which based on Eq.  Fig. 6. According to the figure our theoretical solution and the numerical results agree very well. Solving the PDE (4.15) for large c 2 s /α changes the behaviour of the system from a divergent to a stable one. For example for the limit α → 0 and c 2 s = 0 we have a stable wave solution.
Similar to the case of spherical symmetry, for c 2 s > 0 one can use a Hopf-Cole transformation to obtain a particular exact solution of Eq. (4.15), where C is an arbitrary constant. In this solution the spatial profile is again constant, as in the spherically symmetric case. However, we find numerically that the spatial profile evolves for more generic initial conditions. But the existence of non-singular solutions for c 2 s > 0 indicates that the limit of small speed of sound needs to be taken with great care, and that further investigations into the physical behaviour in this limit may be warranted 3 .

Conclusions and discussion
Our results show that the EFT description of k-essence DE breaks down for part of the parameter space due to a non-linear instability triggered by a term ∝ ( ∇π) 2 in the equations of motion. Based on our numerical analysis using realistic cosmological simulations we see that this instability happens only when the speed of sound c s is small such that the stabilizing linear term c 2 s ∇ 2 π is suppressed. To gain some analytical insight we investigated the PDE in simplified setups. First, we considered a spherically symmetric scenario and showed that the non-linear term ( ∇π) 2 does indeed lead to a blow-up and a similar behaviour to what we saw in the three-dimensional simulations. We specifically showed that the relation between the blow-up redshift and the initial density of the potential well is similar to what we find in cosmological simulations which implies that the instability is found correctly. We further, in connection with a mathematical discussion, studied a simplified PDE considering only the two important terms. We showed that the system, similar to our simulation results, is unstable for vanishing speed of sound. Moreover, we find numerically that the stability is gained for large values of speed of sound. We derived an ODE for the curvature of the minima and showed that the curvature goes to infinity in a finite time. We compare the numerical 1+1 D solution of the PDE at the minimum point, with the ODE and the blow-up time prediction and find consistent results.
The non-linear instability we found does not appear in the linearized theory, where the evolution is stable for all values of the speed of sound. The presence of such an instability shrinks the w − c 2 s parameter space for the healthy k-essence type theories when treated within the EFT framework. Moreover, similar terms appear in non-linear parametrisations of the EFT framework in the context of more general theories which can therefore suffer from related instabilities. As a result these theories have to be considered more carefully, particularly when ( ∇π) 2 appears in the scalar field equation of motion.
The breakdown of the EFT approach can be either due to the EFT truncation order where higher-order corrections can remove the instability, or it can be a hint of the full theory breakdown. In the case of the latter, it also leads to the breakdown of the weak-field approximation and requires a more careful analysis to decide whether coupling the scalar field to gravity could hide the singularities behind (black hole) horizons without a complete, global breakdown of the evolution. It is difficult to assess whether a tiny but non-vanishing c 2 s is able to prevent a singularity, but even if this is the case the solutions will strongly depend on non-linearities and the truncation of the EFT is still rendered inconsistent.  In particular, we first check the robustness of the time integration of the N -body code, and then discuss extensively how the blow-up redshift depends on the spatial resolution of the simulation. The latter is relevant because increasing the resolution results in higher probability for deeper potential wells and changes the initial conditions for the scalar field accordingly. Finally we discuss how much the critical value of the speed of sound c 2 s ∼ 10 −4.7 depends on our cosmological assumptions and validate that the stability of the system is obtained for large values of c 2 s even if we consider matter domination and go far beyond the present epoch. These tests, in addition to our simplified scenarios (spherical and planar symmetries) in Sec. 4, rule out the possibility of the instability being an artefact.

A.1 Precision of the time integration
In this part we discuss the effect of the precision of the time integration on our results. In Fig. 7, for a specific case where N pcl = N grid = 256 3 , L = 300 Mpc/h and c 2 s = 10 −7 , we show the sensitivity of the blow-up redshift to the time resolution of the simulation. For each precision setting we consider two simulations, one where DE perturbations source other components (stars), and one where this is not the case, i.e. DE is a spectator field (circles). Based on these results, even for the largest dτ (the lowest precision considered) the blow-up redshift does not change significantly. This also shows that the blow-up redshift (for high enough time resolution) does not depend on gravity being sourced by the DE component. In summary, our test indicates that the time precision considered in our cosmological simulations is sufficient to resolve the blow-up phenomenon.

A.2 Spatial resolution
In this subsection we discuss the dependence of the blow-up redshift on the spatial resolution of our simulations. Contrary to the time precision, the spatial resolution is not only a precision parameter of the equation, but also a parameter which changes the initial conditions. Higherresolution simulations probe smaller scales where the amplitude of perturbations is larger, and as a result the scalar field evolution is sourced by deeper potential wells. In Fig. 3 we show the blow-up redshift for different maximum values of the initial density, max(δ ini ), where we obtain a linear relation between the two, i.e. 1 + z b ∝ max(δ ini ) in the matter-dominated era. In Fig. 8 we show the measurement of the blow-up redshift for different spatial resolutions. Our fit over the data in matter domination shows 1 + z b ∝ N 0.22 grid where N grid is the total number of grid points in 3D. Our results in Fig. 3 and Fig. 8, considering only the data in matter domination then suggest In cosmological N -body simulations one can estimate the maximum of the initial matter density contrast, max(δ ini (x)), without measuring it directly through the snapshots by invoking the Fisher-Tippett-Gnedenko extreme value theorem [44]. The extreme value theorem describes the distribution of extrema in a similar fashion to how the central limit theorem concerns the behaviour of averages. It states that for a large number N draws from a normal distribution 4 with average µ and standard deviation σ, the maximum of the sample follows (R = 3dx) 2log(N grid ) N-body 2log(N grid ) N-body direct measurment Figure 9: The maximum of the initial density contrast at z = 100 as a function of the number of grid points is shown. The orange stars correspond to the direct measurement of max(δ ini ) on the initial snapshot of the simulation. The green and blue circles represent max(δ ini ) computed using the formula (A.2) for the mean of the Gumbel distribution. The standard deviation used for the green circles is measured directly from the snapshots while for the blue circles we use Eq. (A.3) and using the linear power spectrum from linear Boltzmann code CLASS [43].
a Gumbel distribution with a mean approximated as In a cosmological simulation we may assume that at early times the density contrast follows a normal distribution N (µ, σ) with vanishing mean, µ = 0, and a standard deviation σ R that depends on the resolution of the simulation, corresponding to a smoothing scale R, and is given by Here W (k, R) is the Fourier transform of a window function with radius R, and P (k) is the power spectrum. We assume a top-hat window function with radius R, with CDM simulation Einstein-de Sitter Figure 10: The blow-up redshift for different speeds of sound is shown. The blue circles represent the ΛCDM scenario, while the red stars show the matter domination case i.e. Ω m = 1 where we choose h = 0.3775 such that w m /h 2 = 1. In the ΛCDM case letting the simulation run to the future (negative redshifts) does not decrease the blow-up threshold for the speed of sound. This happens because in the DE dominated era the potential wells decay and help to stabilise the system. However, in the Einstein-de Sitter scenario, going to negative redshifts can increase the threshold to c 2 s ≈ 10 −4.35 .
N grid will lead to smaller R and thus to higher wavenumbers contributing in the integral, resulting in a larger variance. Based on our numerical measurements, for the choice R = √ 3dx the result of the integral (A.3) agrees with our numerical measurements of the variance. In Fig. 9 we show the numerical results for max(δ ini ) computed in three different ways. The orange stars represent the direct measurement of max(δ ini ) in the initial snapshot of the N -body simulation for a given N grid , the green and blue circles represent the results when we use the Gumbel distribution (A.2) to compute max(δ ini ) where in the latter case we use σ R= √ 3dx computed through the integral (A.3) whereas for the former case we use σ as measured directly from the snapshot of the N -body simulation. In both cases we assume N = N grid , even though the draws are not entirely independent. The figure shows excellent agreement between different approaches for computing max(δ ini ). As a result, we can estimate the maximum of the initial density in an N -body simulation using Eq. (A.2) where the variance is computed by Eq. (A.3).
A.3 Cosmology dependence of c * s In Sec. 2 we showed that there is a critical speed of sound c * s such that the system remains stable for c s > c * s . As the blow-up redshift generally decreases with increasing c s , there is the possibility that the critical value is connected to the onset of DE domination where the stability of the system is guaranteed due to decay of potential wells. In this subsection we rule out this possibility by considering simulations in Einstein-de Sitter cosmology, i.e. where matter always dominates. In order to simulate a Universe with Ω m = 1 while being close to the ΛCDM Universe at early times (including the radiation era) we choose a different Hubble parameter h = 0.3775 such that ω m /h 2 = 1. In Fig. 10 we show the results for c 2 s when we have matter domination compared to the ΛCDM scenario. In the ΛCDM case, due to the late-time DE domination, the PDE does not blow up in the future as the potential wells decay and we obtain the critical value of c 2 s ≈ 10 −4.7 . On the other hand, in Einstein-de Sitter the system could still blow up in the future. However, there is still a critical value for the speed of sound c 2 s ≈ 10 −4.35 where for larger values of the speed of sound the system remains stable and does not blow up even in the far future. This observation is in agreement with the claim that for large c s the stability of the system is restored due to the pressure term in the equation, and is not an effect of late-time DE domination.