# Instabilities Appearing in Effective Field theories: When and How?

Jean-Pierre Eckmann, Farbod Hassani, Hatem Zaag

► To cite this version:

Jean-Pierre Eckmann, Farbod Hassani, Hatem Zaag. Instabilities Appearing in Effective Field theories: When and How?. Nonlinearity, 2023, 36 (9), pp.4844-4861. 10.1088/1361-6544/ace769. hal-03873076

HAL Id: hal-03873076

https://hal.science/hal-03873076v1

Submitted on 16 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
PAPER • OPEN ACCESS

# Instabilities appearing in cosmological effective field theories: when and how?

To cite this article: Jean-Pierre Eckmann et al 2023 *Nonlinearity* **36** 4844

View the [article online](#) for updates and enhancements.

## You may also like

- Stabilization despite pervasive strong cross-degeneracies in a nonlinear diffusion model for migration–consumption interaction

Michael Winkler

- Positive solutions of the Gross-Pitaevskii equation for energy critical and supercritical nonlinearities

Dmitry E Pelinovsky, Juncheng Wei and Yuanze Wu

- Asymptotics of the deformed higher order Airy-kernel determinants and applications

Jun Xia, Yi-Fan Hao, Shuai-Xia Xu et al.
# Instabilities appearing in cosmological effective field theories: when and how?

Jean-Pierre Eckmann¹,*, Farbod Hassani²,³
and Hatem Zaag⁴

¹ Département de Physique Théorique and section de Mathématiques, University of Geneva, Geneva, Switzerland

² Institute of Theoretical Astrophysics, Universitetet i Oslo, 0315 Oslo, Norway

³ Department of Physics, McGill University, 3600 rue University, Montreal QC H3A 2T8, Canada

⁴ Université Sorbonne Paris Nord, LAGA, CNRS (UMR 7539), F-93430 Villetaneuse, France

E-mail: Jean-Pierre.Eckmann@unige.ch

Received 2 May 2022; revised 11 April 2023

Accepted for publication 16 June 2023

Published 2 August 2023

Recommended by Dr Christian Klein

<figure><img src="image_2.png" /></figure>

CrossMark

## Abstract

Nonlinear partial differential equations appear in many domains of physics, and we study here a typical equation which one finds in effective field theories originated from cosmological studies. In particular, we are interested in the equation $\partial_t^2 u(x,t) = \alpha(\partial_x u(x,t))^2 + \beta \partial_x^2 u(x,t)$ in $1+1$ dimensions. It has been known for quite some time that solutions to this equation diverge in finite time, when $\alpha > 0$. We study the nature of this divergence as a function of the parameters $\alpha > 0$ and $\beta \ge 0$. The divergence does not disappear even when $\beta$ is very large contrary to what one might believe (note that since we consider fixed initial data, $\alpha$ and $\beta$ cannot be scaled away). But it will take longer to appear as $\beta$ increases when $\alpha$ is fixed. We note that there are two types of divergence and we discuss the transition between these two as a function of parameter choices.

<figure><img src="image_3.png" /></figure>

<footnote>* Author to whom any correspondence should be addressed.</footnote>
The blowup is unavoidable unless the corresponding equations are modified.
Our results extend to 3 + 1 dimensions.

Supplementary material for this article is available [online](#)

Keywords: blowup solutions, cosmological context,
nonlinear partial differential equations

Mathematics Subject Classification numbers: 35B40

(Some figures may appear in colour only in the online journal)

# 1. Introduction

In physics, effective field theories (EFTs) are employed to describe fundamental theories at the low energy limit in a unified form. This approach is widely used to express different physical phenomena [13, 16–18, 27]. This framework, which is constructed based on perturbative expansion, usually leads to non-linear partial differential equations.

Recently, the EFT approach has become very popular in cosmological studies, especially to study the late time accelerating expansion of the Universe which is driven by the so-called dark energy component [6, 14]. Based on cosmological observations [2, 4, 32] the clustering of the dark energy component is supposed to be small, so the linear approximations are justified [25, 37]. However, in the near future high precision measurements of the Universe will be done by the new cosmological surveys [3, 6, 30, 36]. This motivated cosmologists [15, 20, 21] to study non-linear partial differential equations (PDEs) arising from these EFT approaches to have more accurate predictions of these theories. Thus, cosmological N-body simulations have been developed [1, 23] which describe the evolution of structures in the Universe by solving Einstein field equations as well as a non-linear PDE for the dark energy component. In a nutshell, this non-linear PDE has non-linearities that sometimes are dominated by a $(\partial_x u(x,t))^2$ term [24]. This motivates our current study.

Using extensive numerical simulations with *k*-evolution [23] it was discovered earlier in [24] that the solutions of such equations can form violent singularities at *finite* time. And, depending on the cosmological parameters, these singularities can even appear at a time *before* the current epoch of the Universe. Obviously, this asks for a change of parameters, or for regularising the models with additional smoothing terms. In [24] it is specifically argued that appearance of the Laplace term $\beta \partial_{xx} u(x,t)$ with large enough $\beta$ makes the system stable. It is worth noting that in cosmological context $\sqrt{\beta} = c_s$ where $c_s$ is the ‘speed of sound’ and determines how fast the perturbations of scalar field propagate. In cosmology usually, $0 \le c_s \le 1$. However superluminal cases ($c_s \ge 1$) are also considered in the literature [10, 12]. However, we will show that the Laplacian term does not regularise enough to avoid the finite time blowup: it just shifts the blowup time to a later epoch. Perhaps, not seeing the instability in the realistic cosmological N-body simulations might be due to the short time period, or some other phenomena which are present in cosmological setups.
<figure><img src="image_4.png" /><figcaption>Figure 1. The V-type divergence: the solution profile $u(x,t)$ is shown at different times (different colours), up to the blowup time. The curvature at the minimum diverges in finite time. The initial conditions are: $u(x,0) = 0$ and $\partial_t u(x,0) = -\exp(-30x^2)$ on $[-\pi, \pi]$.</figcaption></figure>

In the paper [24], employing the cosmological N-body code k-evolution, it was found that the PDE for the EFT of dark energy (in particular, k-essence models<fnref n="5" />), for some set of parameters leads to an instability in finite time. It was shown in [33] that the main source of finite-time instability is due to the presence of the non-linear term in

$$ \frac{\partial^2 u(x,t)}{\partial t^2} = \alpha \cdot \left( \frac{\partial u(x,t)}{\partial x} \right)^2, \qquad (1.1) $$

which appears naturally in EFT theories (in cosmological studies $u(x,t)$ is called $\pi(x,t)$). Here, we consider only fixed $\alpha > 0$ but for general theories $\alpha$ can be time dependent. This time dependence appears in k-essence theories through the Hubble parameter $H(t)$ [24] and we have checked that the divergence persists when taking into account the time dependence of $H$. The local in space existence is guaranteed by the Cauchy–Kowalevski theorem, for (locally) analytic initial conditions. A global in space (but of course finite in time) existence can be obtained by considering initial conditions which are analytic in a strip around the real axis, by reducing the widths of the strip in time<fnref n="6" />.

Beyond this setting, it is easy to see that (1.1) has solutions which diverge in finite time [24]: in fact, for some initial conditions we can consider $u$ of the form $u(x,t) = f(t)x^2$. This leads to $f''(t) = 4\alpha \cdot f(t)^2$. When the initial condition is $f(0) = \frac{3}{2\alpha t_0^2}$ and $f'(0) = \frac{9}{\alpha t_0^4}$, then we have

$$ f(t) = \frac{3}{2\alpha(t_0 - t)^2}, \qquad (1.2) $$

<footnote n="5">⁵ These are a general class of theories in which the action contains at most one temporal and one spatial derivative acting on the field [7]. The k-essence models have been proposed as a possible explanation for the late time cosmic acceleration. In the k-evolution code, the k-essence field and other cosmological components, such as dark matter and baryonic matter, are implemented. The 'k' in 'k-essence' stands for 'kinetic,' which refers to the kinetic energy of the theory. The idea behind k-essence is that the negative pressure of its fluid description is caused by its kinetic energy.</footnote>

<footnote n="6">⁶ We thank Pierre Collet and Martin Hairer for helpful remarks on this point.</footnote>
<figure><img src="image_5.png" /><figcaption>Figure 2. The M-type divergence: as time advances, the support of the function spreads, and the function at the edge gets steeper, until the derivative will diverge at some time $T_*$. The simulation is for the equation (1.3), with $\alpha = 0.15$ and $\beta = 0.05$ and initial conditions $\partial_t u(x, 0) = -\exp(-30x^2)$, $u(x, 0) = 0$ on $[-\pi, \pi]$.</figcaption></figure>

which diverges as $t \uparrow t_0$, when $\alpha > 0$. As discussed in [33–35], this divergence is of a local type and the curvature of the minima increases to infinity in a finite time. Here we call this divergence the ‘V’ type as shown in figure 1. We give a general divergence proof in section 2. Consider now the equation

$$ \frac{\partial^2 u(x,t)}{\partial t^2} = \alpha \cdot \left( \frac{\partial u(x,t)}{\partial x} \right)^2 + \beta \frac{\partial^2 u(x,t)}{\partial x^2}. \qquad (1.3) $$

When $\beta = 0$ this is (1.1). We next consider the case when $\beta > 0$ (and $\alpha > 0$ is fixed). The local in time Cauchy problem for this equation can be solved in $W^{1,\infty} \times L^\infty(\mathbb{R})$ (see [8] for a similar strategy). From the finite speed of propagation, for given initial data at $t = 0$, two cases appear for the domain of definition of the solution:

* either it is $\{(x,t) | t \ge 0\}$ and the solution is said to be global;

* or it can be expressed as

$$ \mathcal{D} = \{(x,t); 0 \le t < T(x)\} $$

for some Lipschitz function $x \mapsto T(x)$ (with Lipschitz constant $1/\sqrt{\beta}$), where the solution is said to ‘blow up’ in finite time. (see appendix B for details). In this paper, we deal with blow-up solutions.

As for the blow-up behaviour, there are two scenarios: when $\beta$ is very small, the solution will be of V-type but when $\beta$ is larger, then it will be of a shape we call M-type, as illustrated in figure 2.

It was seen in simulations presented in [24], that the instability seems to vanish when one adds a Laplace term with large enough coefficient $\beta$.

However, this is not the whole story: in fact, whenever $\alpha > 0$ this system is always unstable and even large $\beta$ cannot cure the instability forever, (this holds rigorously for a large class of initial conditions, and seems to hold in numerics, for any initial condition with compact support). However, increasing $\beta$ does increase the blowup time. Actually, the blowup now
happens at the ‘ends’ of the ‘M’ shape in the profile, whose walls get steeper and steeper, until the derivative becomes infinity. This phenomenon was known for some time in the literature [26, 29], and it is quite generic. In simulations, perhaps one does not wait long enough, or the M gets too wide in the numerically allocated spatial direction, before the divergence happens.

In the following sections, we will study in more detail the domains in the $\alpha, \beta$ plane for which ‘V’ and ‘M’ divergences will happen and the time it takes for blowup to occur. This will tell us for which physical parameters the singularity is so far in the future that it can be neglected, or that it is of a form which can easily be damped by adding additional terms to (1.3). However, such terms do not seem to occur naturally in the case of cosmological models for dark energy, especially within the weak-field approximation [20]. For the convenience of the reader, we repeat in two appendices some details about the V type divergence, and we also repeat—with small variations—a proof of the persistence of divergence for all $\beta > 0$ (when $\alpha > 0$). These appendices are based on [34] and [28, 29].

## 2. Blowup time as a function of $\alpha$ when $\beta = 0$

Before we can study the dependence on $\beta$, we need to study the divergence time for the case $\alpha > 0$. The following is a slight adaptation of the results of [33, 34].

We consider the equation $u_{tt} = \alpha(u_x)^2$ on the real line. We start by writing the solution in the form

$$u(x,t) = f(x) + g(x)t + \alpha \int_{0}^{t} d\tau \int_{0}^{\tau} d\tau' (u_x(x, \tau'))^2. \qquad (2.1)$$

This corresponds to the initial conditions

$$u(x,0) = f(x), \quad u_t(x,0) = g(x).$$

We will consider the case where $f'(0) = g'(0) = 0$, and we ask how the solution behaves near $x = 0$. Depending on the curvatures $f''(0)$ and $g''(0)$, the second derivative $u_{xx}(x,t)$ will, or will not diverge at $x = 0$. Of course, if the functions $f$ and $g$ have vanishing derivatives at some other point(s) $x_0$, the same discussion will apply at those points, and there can be one of these points where $u_{xx}(x_0, t)$ diverges before the one at $x = 0$. In the following proposition, we will neglect this aspect.

**Proposition 2.1.** Assume $f'(0) = g'(0) = 0$. Define

$$c = \frac{1}{2} g''(0)^2 - \frac{2}{3} \alpha f''(0)^3. \qquad (2.2)$$

Then the following cases appear:

(i) If $g''(0) > 0$ then $u_{xx}(0, t)$ diverges in finite time $t_+$ given by

$$t_+ = \int_{f''(0)}^{\infty} \frac{db}{\sqrt{\frac{4}{3}\alpha b^3 + 2c}} = \int_{f''(0)}^{\infty} \frac{db}{\sqrt{\alpha(\frac{4}{3}(b^3 - f''(0)^3) + g''(0))^2}}. \qquad (2.3)$$

(ii) If $g''(0) < 0$ then $u_{xx}(0, t)$ will converge to $b_*$ in a finite time $t_-$, where

$$\frac{2}{3}\alpha b_*^3 = -c, \text{ and } t_- = \int_{b_*}^{f''(0)} \frac{db}{\sqrt{\frac{4}{3}\alpha b^3 + 2c}}. \qquad (2.4)$$

At this point in time, we will have $u_t(0, t_-) = 0$ which corresponds to (iii) and the solution will diverge after another finite time $t_+$ (unless $c = 0$).
(iii) If $g''(0) = 0$, and $f''(0) \neq 0$ then $u_{xx}(0, t)$ diverges in finite time $t_+$ given again by (2.3).

(iv) If $g''(0) = 0$ and $f''(0) = 0$ then $u_{xx}(0, t)$ stays constant.

**Remark 2.2.** When $g''(0) > 0$ and $f''(0) = 0$, then the elliptic integral can be evaluated explicitly and one gets

$$t_+(\alpha, g''(0)) = \frac{A}{\alpha^{1/3}g''(0)^{1/3}},$$

with $A \approx 2.5479$.

**Remark 2.3.** Assume that $u(x, 0) = 0$ and $u_t(x, 0)$ is a smooth, bounded function with several well-separated extrema. Such initial conditions are typical for questions in cosmology. In this case, the blowup will happen first in that point $x_0$ for which $C \equiv \alpha u_{txx}(x_0, 0)$ is maximal (and $u_{tx}(x_0, 0) = 0$). In that case, the formula of remark 2.2 leads to $t_+ \sim 2.547/C^{1/3}$.

**Remark 2.4.** Note that $t_+ + t_-$ is the total time for an initial condition $g''(0) < 0$ to diverge (only when $c \neq 0$, with $c$ defined in (2.2)). And then, the divergence time is $t_-(f''(0), g''(0)) + t_+(b^*)$.

The proof of these statements is given in appendix A.

### 3. Crossover

The equation (1.3) can be rescaled in space and time, to eliminate either $\alpha$ or $\beta$. This rescaling will be used in appendix B, when we show that blowup is unavoidable in many cases. Here, we ask a different question as in cosmological studies it is important to know the blowup as a function of $\alpha$ and $\beta$ for a fixed initial perturbation. And then, the scaling cannot be used.

As we have indicated in the introduction, for a fixed $\alpha > 0$ there will be a crossover between the V-type and the M-type divergence as one varies $\beta$. This crossover is important in general, because the nature of the divergence is different in the two cases: for the V-type divergence, we see a localised blow-up of the second derivative in a minimum, as shown in figure 1. In the second, M-type case, the divergence happens at some distance from the critical point of the initial condition (see figure 2). A vertical wall forms, and this wall moves outward from the centre, basically with the propagation speed of the wave (which is $\sqrt{\beta}$ for (1.3)). The rate of divergence for the case when $\beta > 0$ differs from that for $\beta = 0$.

The divergences in the EFT framework are mainly important, because they are either a hint for entering the strong field regime (where the perturbative expansion is not valid anymore) or a breakdown of the underlying fundamental theory [24]. The divergence type should correspond to physical phenomena happening at high energy/small scales. Especially it may help to introduce appropriate mechanisms to remove such instabilities. As an example the V-type divergence is localised, and in cosmological studies (as suggested in [24]) can be used as an origin of super massive black holes. On the other hand the M-type blowup resembles caustics formation in the Universe [9, 11]. These two cases have different signatures in cosmological observables [31]. And therefore it is useful to distinguish them.

We illustrate the various possibilities in figures 3–5. In figure 3 we consider initial conditions with a minimum, namely $-\exp(-30x^2)$. In this case, depending on $\alpha$ and $\beta$ we can see a V-type divergence or an M-type divergence. When $\alpha$ is large and $\beta > 0$ not too small, we detect clearly an M-type divergence. But when $\beta$ is too close to 0, the numerics breaks down, and one sees something like a V-type divergence in the region ‘V’ of figure 3. However, it is also possible that in fact for all $\beta > 0$ one has M-type divergence, with the two walls of the M
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Distance of singularity at divergence time as a function of $\alpha$ and $\beta^{1/2}$: in the region marked ‘V’ the numerics does not allow to distinguish between V-type and M-type divergence, since the distance of the (two) singularities is too small to distinguish M and V. The other white region is the set of parameters where the simulation hits the boundary before divergence (see appendix C for an explanation). The initial condition is $u(x, 0) = 0$ and $u_t(x, 0) = -\exp(-30x^2)$. The colour coded part shows the distance of the maximum of $|u_{xx}|$ at blowup time.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Divergence time as a function of $\alpha$ and $\beta$: the x axis is $\beta$. The initial condition is $u(x, 0) = 0$ and $u_t(x, 0) = +\exp(-30x^2)$. The white region contains those parameters ($\alpha$ and $\beta$) for which the solution hits the boundary before divergence, as in figure 3. Each colour corresponds to a different divergence time. Note that in this case, since the initial condition satisfies the condition of theorem B.1 (we made the support finite), we *know* that the solution must blow up in finite time, for all $\alpha > 0$ and $\beta \ge 0$.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Distance of singularity at divergence time as function of $\alpha$ and $\beta^{1/2}$: because the initial condition is now $u(x, 0) = 0$ and $u_t(x, 0) = +\exp(-30x^2)$, all divergences seem to be of M-type (see theorem B.1, which only asserts divergence in finite type, but the proof suggests divergence at near the advancing front). Note that for all $\alpha$ and $\beta > 10^{-6}$ considered in this graph, the distance of the singularity from $x=0$ at blowup is at least 0.1795 (in the interval $[-\pi, \pi]$). For this initial condition and the $\alpha, \beta$, considered, blowup always happens before the wave hits the boundary.</figcaption></figure>

so close together that the numerics gets unreliable. With our current understanding we conjecture that it is always M-type when $\beta > 0$. Note that, in any case theorem B.1 shows that for all $\beta > 0$ there is finite time divergence (the proof only applies to initial conditions of a certain positivity type—the function $M$ of appendix B, but we have seen the M-type divergence for general non-trivial initial conditions with compact support).

Indeed, when $\alpha$ is large and $\beta > 0$ is small, then the minimum will get more and more pointy as in figure 1, until the second derivative diverges. This happens in the white region 'V' in figure 3. In the coloured parameter regions, the negative initial conditions get more pointy, until the Laplacian term regularises the central part, which then grows, becomes positive and then diverges away from the centre, as in figure 2.

When the initial condition is positive, with a local maximum, then the situation is somewhat simpler, because the transition from negative to positive is absent. No V singularity can form. Note that this is consistent with the discussion of the cases (i) and (ii) in proposition 2.1. The divergence time as function of $\alpha$ and $\beta$ is shown in figure 4 for the divergence time, and in figure 5 for the divergence distance, by which we mean the distance of the singularity from the coordinate origin. We cannot offer a formula for the curves in either of the figures.

**4. Conclusions and discussion**

In this article we discussed the equation $\partial_t^2 u(x,t) = \alpha(\partial_x u(x,t))^2 + \beta \partial_x^2 u(x,t)$. These types of PDEs appear naturally in the EFT descriptions of physical systems, where one approximates the equations assuming weak fields. We specifically discuss the two divergence types arise from these equations, namely the 'V' and the 'M' type. We discuss how and when these
instabilities are generated. In some cosmological studies [22] it is suggested that the term $\beta\partial_x^2 u(x,t)$ stabilises the system so that the instability caused by $\alpha(\partial_x u(x,t))^2$ vanishes. This stabilisation is especially motivated by realistic cosmological studies. While our results are in agreement with [22, 24, 33] regarding the V-type blowup which happens for small $\beta$, our results show that even for large $\beta$ the instability is unavoidable, however it is always of the new M-type. The reasons that this instability is not seen in the realistic cosmological simulations for large $\beta$ (or $c_s^2$ in cosmology) could probably be due to 1) the role of gravity which is neglected in our study 2) the difference in boundary conditions or 3) the instability exists but will appear beyond the times considered in a study. Preliminary studies show that the results carry over to the 3 + 1 dimensional situation<fnref n="7" />.

## Data availability statement

No new data were created or analysed in this study.

## Acknowledgments

We thank Edriss Titi to have brought us together for this paper. We profited from discussions with Julian Adamek, Pierre Collet, Martin Kunz, Sabir Ramazanov, Pan Shi, Peter Wittwer, and Alexander Vikman. This work was partially supported by an ERC advanced Grant (Bridges 290843) and by SNF Swissmap.

## Appendix A. Divergence time when $\beta = 0$

Here we prove proposition 2.1.

**Proof.** Define

$$a(t) = u_x(0, t), \quad b(t) = u_{xx}(0, t).$$

Then (2.1) leads to

$$a(t) = f'(0) + g'(0)t + 2\alpha \int_{0}^{t} d\tau \int_{0}^{\tau} d\tau' u_x(0, \tau') u_{xx}(0, \tau'),$$

and therefore

$$\ddot{a}(t) = 2\alpha a(t) b(t). \qquad (A.1)$$

Similarly,

$$b(t) = f''(0) + g''(0)t + 2\alpha \int_{0}^{t} d\tau \int_{0}^{\tau} d\tau' ((u_{xx}(0, \tau'))^2 + u_x(0, \tau') u_{xxx}(0, \tau')).$$

From this, we deduce

$$\begin{aligned} \ddot{b}(t) &= 2\alpha(u_{xx}(0, t))^2 + 2\alpha u_x(0, t)u_{xxx}(0, t) \\ &= 2\alpha b(t)^2 + 2\alpha a(t)u_{xxx}(0, t). \end{aligned} \qquad (A.2)$$

<footnote n="7">⁷ In the 3 + 1 dimensional case with spherical symmetry, the equation (1.3) is simply replaced by $u_{tt}(r,t) = \alpha(u_r(r,t))^2 + \beta(u_{rr}(r,t) + \frac{2}{r}u_r)$, where r is the radial coordinate. We have done the corresponding numerical experiments for this case.</footnote>
Since we assume $f'(0) = g'(0) = 0$ we find from (A.1) that $a(t) = 0$ for all $t$ for which $b(t)$ is finite. Therefore, (A.2) reduces to

$$ \ddot{b}(t) = 2\alpha(b(t))^2. \qquad (A.3) $$

We will discuss this equation. For computing the divergence time, it is useful to transform the equation as follows: Multiplying by $\dot{b}$ leads to

$$ \frac{1}{2} \frac{d}{dt} (\dot{b}(t))^2 = \frac{2}{3} \alpha \frac{d}{dt} b(t)^3, $$

or, for some $c$,

$$ \frac{1}{2} (\dot{b}(t))^2 = \frac{2}{3} \alpha (b(t))^3 + c. \qquad (A.4) $$

Note that looking at $t=0$ we find

$$ c = \frac{1}{2} \dot{b}(0)^2 - \frac{2}{3} \alpha b(0)^3 = \frac{1}{2} g''(0)^2 - \frac{2}{3} \alpha f''(0)^3, \qquad (A.5) $$

which is the definition in the proposition. Note that

$$ b(0) = f''(0), \quad \text{and} \quad \dot{b}(0) = g''(0). $$

We consider first the case where $g''(0) > 0$. Then $\dot{b}(0) > 0$ and from (A.4) we find that

$$ \dot{b}(t) = \sqrt{\frac{4}{3}\alpha(b(t))^3 + 2c}, \qquad (A.6) $$

which means $b$ is increasing and the quantity below the square root is always positive. Using standard techniques, we get

$$ dt = \frac{db}{\sqrt{\frac{4}{3}\alpha b^3 + 2c}}. $$

From (A.6) we deduce the divergence time $t_+$,

$$ t_+ = \int_{b(0)}^{\infty} \frac{db}{\sqrt{\frac{4}{3}\alpha b^3 + 2c}}. \qquad (A.7) $$

This proves (2.3).

The case $g''(0) < 0$ is handled similarly, but now (A.6) is replaced by

$$ \dot{b}(t) = -\sqrt{\frac{4}{3}\alpha b(t)^3 + 2c}. \qquad (A.8) $$

This means that $b$ is decreasing until the square root in (A.8) vanishes. This defines $b_*$, and then (A.7) is replaced by

$$ t_- = \int_{b_*}^{b(0)} \frac{db}{\sqrt{\frac{4}{3}\alpha b^3 + 2c}}. $$

This leads to (2.4).

The assertions under (iii) are a simple variant of (i) and (ii). The difference is that because $g''(0) = 0$, we find now that $c = -\frac{2}{3}\alpha f''(0)^3$, and $c \neq 0$ by the assumption $f''(0) \neq 0$. Note that in this case, the positivity of $\dot{b}(t)$ follows from the second order ODE (A.3), given that $b(0) = f''(0) \neq 0$ and $\dot{b}(0) = g''(0) = 0$. The only remaining case is (iv), $g''(0) = f''(0) = 0$, which implies that $b(0) = \dot{b}(0) = 0$, hence, directly leads to $b(t) = b(0) = 0$ by (A.3). Note
that the proof says nothing about the divergence of $u_x$ or $u$. Depending on the initial condition, one can see a divergence of $u_x$ or of $u$ (especially if $\beta > 0$).

<figure><img src="image_9.png" /></figure>

## Appendix B. Finite time divergence

We adapt here the proof of [29] to the 1 + 1 dimensional context. The proof actually works in the same way in higher dimensions, for which it was already spelled out in [28, 29]. Since we deal here only with the fact that the equation will diverge in finite time, it suffices to consider instead of $\alpha > 0, \beta > 0$, the simpler form

$$u_{tt} - u_{xx} = u_x^2. \qquad (B.1)$$

Indeed, if $v(x,t)$ solves $v_{tt} = \alpha(v_x)^2 + \beta v_{xx}$, then $u(x,t) = \frac{\beta}{\alpha}v(\beta^{-1/2}x, t)$ satisfies (B.1). Note that we work in $\mathbf{R}$, and not, as happens in some simulations, in periodic boundary conditions. Assume the initial conditions are

$$u(x, 0) = f(x), \quad u_t(x, 0) = g(x),$$

with $f, g$ having support in $|x| < X$. We may assume $(f,g) \in H$ for some functional space, for example $H = W^{1,\infty}(\mathbf{R}) \times L^\infty(\mathbf{R})$.

Using a standard fixed point technique, we can find a local in time solution $u(t) \in C([0, t_0], H)$ for some small $t_0 > 0$. From the finite speed of propagation, using a cut-off technique, this global (in space) existence result extends to some local in time existence result in slices of backward cones of slope 1. This way, we can see that our solution is defined beyond the above-mentioned strip $\mathbf{R} \times [0, t_0]$, and extends to a larger domain of definition, which happens to be a union of backward light cones, with different heights. From elementary considerations, one of the following cases occurs:

(i) Either the union is the half-space $\mathbf{R} \times [0, \infty)$. We say in that case that the solution is 'global' (for forward time).

(ii) Or, the union writes as

$$\{(x,t) | 0 \le t < T(x)\}$$

for some 1-Lipschitz function $T: \mathbf{R} \to \mathbf{R}$. In that case, we say that $u$ 'blows up in finite time.' Note that by construction, we have a local blow-up time for each $x \in \mathbf{R}$, namely $T(x)$. For more details on the construction of the domain of definition, see [5] and also [8].

Define now, see [29], for $x > 0$,

$$M(x) = \frac{1}{2}f(x) + \frac{1}{2} \int_{x}^{x} g(\xi) d\xi.$$

**Theorem B.1.** Assume there is an $X_0 \in (0, X)$ for which $M(x) \ge 0$ for all $x \in (X_0, X)$, and also

$$\int_{X_0}^{X} M(\xi) d\xi \equiv \varepsilon > 0.$$

Then, the solution blows up in finite time, in the sense of the definitions above.

**Remark B.2.** Note that the theorem is shown under the assumption that $M(x) > 0$. This covers cases where $u(x, 0)$ and/or $u_t(x, 0)$ are positive (or positive near the edge of their support). The case of negative $M$ is not covered by the literature, nor by our proof. However,
we have studied many cases with negative initial M. For example, the case $u(x, 0) = 0$ and $u_t(x, 0) = -\exp(-C \cdot x^2)$ of figure 3. In all these cases we seem to see 'M-type' divergence. We can consider the solution $u(x, t)$, $u_t(x, t)$ as a new initial condition for any $t$. Then, we have observed that M starts out negative, decreases, and finally crosses 0. From that point on, we are again in the domain of validity of theorem B.1, and, indeed, the solution diverges. (This is reminiscent of the two cases in appendix A.)

**Proof.** It suffices to consider the situation where $u(x, 0)$ and $u_t(x, 0)$ have support in $|x| \le X$, but we will consider only the side of positive $x$ in the sequel. Clearly, $u(x, t) = 0$ for $x \ge t+X$ due to the finite propagation speed. One estimates now the function

$$H(t) = \int_{X_1}^{t} (t - \tau) \int_{\tau+X_0}^{\tau+X} u(\xi, \tau) d\xi \, d\tau.$$

Here, $X_1 = (X - X_0)/2$. From the definition, we get

$$H''(t) = \int_{t+X_0}^{t+X} u(\xi, t) d\xi. \qquad (B.2)$$

One has the explicit formula

$$u(x,t) = u_0(x,t) + \frac{1}{2} \int_0^t \int_{x-t+\tau}^{x+t-\tau} u_x(\xi, \tau)^2 d\tau \, d\xi, \qquad (B.3)$$

with the 'free evolution'

$$u_0(x,t) = \frac{1}{2}(f(x-t) + f(x+t)) + \frac{1}{2} \int_{x-t}^{x+t} g(\xi) d\xi.$$

When $x \ge t + X_0$ and $X \ge t \ge X_1$, then $x+t \ge X$ and therefore $f(x+t) = 0$ and therefore, in this region,

$$u_0(x,t) = \frac{1}{2}f(x-t) + \frac{1}{2} \int_{x-t}^{x+t} g(\xi) d\xi.$$

We get from (B.3) and (B.2),

$$H''(t) = G_0(t) + G_1(t),$$

with

$$G_0(t) = \int_{t+X_0}^{t+X} u_0(x,t) dx = \int_{t+X_0}^{t+X} M(x-t) dx = \int_{X_0}^{X} M(x) dx = \varepsilon,$$

by the definition of $\varepsilon$ in theorem B.1. The nonlinearity leads to

$$G_1(t) = \int_{t+X_0}^{t+X} dx \int_0^t d\tau \int_{x-t+\tau}^{x+t-\tau} d\xi \, u_x(\xi, \tau)^2.$$

In lemma B.3 below, we show that for $t > X_1 \equiv (X - X_0)/2$ one has

$$G_1(t) \ge \frac{1}{t+X} \int_0^t d\tau \int_{\tau+X_0}^{\tau+X} d\xi (t-\tau) (\xi-\tau-X_0) u_x(\xi, \tau)^2. \qquad (B.4)$$

We use now the Schwarz inequality in the form

$$\int \varphi \psi = \int \varphi^{1/2} (\varphi^{1/2} \psi) \le \left( \int \varphi \psi^2 \right)^{1/2} \left( \int \varphi \right)^{1/2},$$
with $\varphi = (t - \tau)(\xi - \tau - X_0)$ and $\psi = u_x$. This leads to

$$G_1(t) \geq F^2(t)/J(t),$$

with

$$F(t) = \int_{0}^{t} \int_{\tau+X_{0}}^{\tau+X} (t-\tau)(\xi-\tau-X_{0})u_{x}(\xi, \tau)d\xi d\tau$$

and

$$J(t) = \int_{0}^{t} \int_{\tau+X_{0}}^{\tau+X} (t-\tau)(\xi-\tau-X_{0})d\tau d\xi = \frac{(X-X_{0})^{2}t^{2}}{4}.$$

If we integrate the expression for $F$ by parts (in $\xi$), we get

$$F(t) = - \int_{0}^{t} \int_{\tau+X_{0}}^{\tau+X} (t-\tau)u(\xi, \tau)d\xi d\tau = -H(t).$$

Therefore, we find finally

$$H''(t) \geq G_0(t) + \frac{H(t)^2}{J(t)}. \qquad (B.5)$$

Fix now $T$ and we will show that the solution cannot exist for $T > T_*$, where $T_*$ will be computed in the proof: we use here lemma 1 from [28] adapted to the 1d case. The ingredients are that

$$H''(t) \geq G_0(t) = \varepsilon > 0, \qquad (B.6)$$

for all $t \geq 0$ and

$$H''(t) \geq G_1(t) \geq 4 \frac{H(t)^2}{(X - X_0)^2 t^2}, \qquad (B.7)$$

for $t > X_1$ (as long as the solution exists). Furthermore, $H(X_1) = H'(X_1) = 0$.

Fix now $T_1 = 2(X_1 + 1)$. Then, for $t > T_1$, we have $t > \frac{1}{2}(t + 1)$, and we replace from now on (B.7) by the simpler

$$H''(t) \geq K_1 \frac{H(t)^2}{(t+1)^2}, \text{ for } t > T_1, \qquad (B.8)$$

with $K_1 = 16/(X - X_0)^2$.

The idea is now to deduce from (B.6) and (B.2) an inequality of the form

$$H'(t) \geq C H^{1+\delta}(t) \text{ for } t > T_1 \text{ with } \delta > 0. \qquad (B.9)$$

This implies divergence in finite time, when $H(T_0) > 0$. Indeed, if $H(T_0) = c^{-1/\delta} > 0$, then

$$H(t) = \frac{1}{(c - C\delta(t - T_0))^{1/\delta}}. \qquad (B.10)$$

One can reformulate this as follows: if $H(T_0) = A$ and $A \leq 1/e$, the optimising $\delta$ in (B.10) is $\leq 1$ and therefore we find that the divergence time is proportional to $-\log(A)$. Note that, if, for example, the leading edge of the support (at $x = 0$) is like $|x|^2$ for $x < 0$, then this will lead to earlier divergence compared to $|x|^3$.

We now begin the proof proper. If $B > 0$, we will use repeatedly the inequality

$$\frac{x}{x+B} \geq \frac{1}{2}, \text{ for all } x \geq B. \qquad (B.11)$$
From (B.6) we find

$$H(t) \geq K_2 \epsilon t^2, \text{ for all } t > 0, \qquad (B.12)$$

with $K_2 = \frac{1}{2}$.

Substituting (B.12) into (B.2), we get

$$H''(t) \geq K_1 K_2 \epsilon H(t) \frac{t^2}{(t+1)^2} \geq K_3 \epsilon H(t), \text{ when } t > T_2, \qquad (B.13)$$

for some large enough $T_2 = \text{const.} T_1$, not depending on $\epsilon$. Since $H'(t) > 0$, we can multiply (B.13) by $H'$ and write it as

$$\frac{d}{dt}(H'(t)^2) \geq K_3 \epsilon \frac{d}{dt}(H(t)^2) \text{ when } t > T_2.$$

We integrate from $T_2$ to $t$ and obtain

$$H'(t)^2 \geq K_3 \epsilon (H(t)^2 + H'(T_2)^2 - H(T_2)^2) = K_3 \epsilon H(t)^2 + K_4 \epsilon \text{ when } t > T_2,$$

for some $K_4$. From (B.6), we conclude that for large enough $T_3$, one has

$$K_3 \epsilon H(t)^2 + K_4 \epsilon \geq K_5 \epsilon H(t)^2, \text{ when } t > T_3.$$

Combining the last two equations we find

$$H'(t) \geq K_5^{1/2} \epsilon^{1/2} H(t), \text{ when } t > T_3.$$

Integrating from $T_3$ to $t$ leads to

$$H(t) \geq H(T_3) \exp\left(K_6 \epsilon^{1/2}(t-T_3)\right) \geq H(T_3) \exp\left(\frac{1}{2} K_6 \epsilon^{1/2} t\right), \text{ when } t > 2T_3, \qquad (B.14)$$

with $K_6 = K_5^{1/2}$. Substituting again into (B.2), we get

$$H''(t) \geq K_7 H(t)^{1+\delta}, \text{ for any } \delta > 0,$$

since the exponential in (B.14) (to the power $\delta > 0$) will dominate the factor $(t+1)^{-2}$ of (B.2), only if $t$ is sufficiently large. (Note that $K_7$ and this new minimal time $T_4$ will depend on $\delta$.) We now multiply the last equation by $H'$ and we obtain

$$\frac{d}{dt}(H'(t)^2) \geq \frac{2K_7}{2+\delta} \frac{d}{dt}(H(t)^{2+\delta}), \text{ for } t > T_4.$$

Integrating from $T_4$ to $t$ we find

$$H'(t)^2 \geq \frac{2K_1}{2+\delta} (H(t)^{2+\delta} - H(T_4)^{2+\delta}) + H'(T_4)^2.$$

Taking square roots on both sides and choosing $T_*$ sufficiently larger than $T_4$, we finally arrive at (B.9) from which we see that there is a divergence in finite time, as in (B.10). □
We still need to show the inequality (B.4).

**Lemma B.3.** Let

$$G_1(t) = \int_{t+X_0}^{t+X} dx \int_0^t d\tau \int_{x-t+\tau}^{x+t-\tau} d\xi u_x(\xi, \tau)^2 .$$

Let $X > X_0 > 0$, and assume $u_x(x,t) = 0$ for all $|x| \ge t+X$, $0 \le t \le T$. Then one has for all $t \ge X_1 \equiv (X - X_0)/2$ the inequality

$$G_1(t) \ge \frac{1}{t+X} \int_0^t d\tau \int_{\tau+X_0}^{\tau+X} d\xi (t-\tau)(\xi-\tau-X_0)u_x(\xi, \tau)^2 . \quad (B.15)$$

**Proof.** We first decompose the triple integration into 3 pieces: the original integration is over the domain

$$\begin{align}
& X_0 + t \le x \le X + t, \nonumber \\
& 0 \le \tau \le t, \tag{B.16} \\
& x - t + \tau \le \xi \le x + t - \tau. \nonumber
\end{align}$$

Also note that the integrand has support in $\xi \le \tau + X$. The three pieces are

$$\begin{align}
& 0 \le \tau \le t - X_0, \nonumber \\
& \tau + X_0 \le \xi \le \tau + X, \tag{B.17} \\
& t + X_0 \le x \le \xi + t - \tau, \nonumber
\end{align}$$

and

$$\begin{align}
& t - X_0 \le \tau \le t, \nonumber \\
& \tau + X_0 \le \xi \le 2t - \tau + X_0, \tag{B.18} \\
& t + X_0 \le x \le \xi + t - \tau, \nonumber
\end{align}$$

and

$$\begin{align}
& t - X_0 \le \tau \le t, \nonumber \\
& 2t - \tau + X_0 \le \xi \le \tau + x, \tag{B.19} \\
& \xi - t + \tau \le x \le \xi + t - \tau. \nonumber
\end{align}$$

One can show that (B.17)–(B.19) defines a domain which coincides with that of (B.16), and that the three regions are disjoint.

We now give lower bounds for the three regions. For (B.17) we get

$$\begin{align*}
& \int_0^{t-X_1} d\tau \int_{\tau+X_0}^{\tau+X} d\xi u_x(\xi, \tau)^2 \int_{t+X_0}^{\xi+t-\tau} dx, \\
& \int_0^{t-X_1} d\tau \int_{\tau+X_0}^{\tau+X} d\xi u_x(\xi, \tau)^2 (\xi - \tau - X_0).
\end{align*}$$

We bound the last factor from below by

$$(\xi - \tau - X_0) \ge \frac{t - \tau}{t + X} (\xi - \tau - X_0).$$

Similarly, for (B.18), the $x$ integration is bounded from below in exactly the same way. Finally, for (B.19), using the support property $\xi \le \tau + X$, we find (since $t > X_1$ and $X_1 < X$),
$$ \xi - \tau - X_0 \le X - X_0 = 2X_1 < t + X_1 < t + X. $$

This leads to a bound for the $x$ integral of the form

$$ \int_{\xi-t+\tau}^{\xi+t-\tau} dx = 2(t-\tau) \ge (t-\tau) \frac{\xi-\tau-X_0}{t+X}. $$

Collecting terms, we finally find that

$$ G_1(t) \ge \frac{1}{t+X} \int_0^t d\tau \int_{\tau+X_0}^{\tau+X} d\xi (t-\tau) (\xi-\tau-X_0) u_x(\xi, \tau)^2. $$

<figure><img src="image_10.png" /></figure>

## Appendix C. Numerics

We integrate all the equations by using the Dorman–Prince [19] Runge–Kutta integrator. The functions are discretised in 2¹³ equidistant points. Derivatives are computed by using five-point stencils. Divergence is defined by $\max(|u_{xx}|) > 10^6$. Special care has been given to assert the quality of the results: we compute with a tolerance (if achievable) of $10^{-11}$. There are two situations where the integration can fail: the time steps get too short (this happens sometimes when $\alpha$ is large and $\beta$ is small). The other problem is the size of the domain in $x$: we take periodic boundary conditions on $[-\pi, \pi]$, and initial data which vanish at these boundaries. If, during time evolution, the value of $|u(\pm\pi, t)|$ exceeds $10^{-4}$, we consider that the wave-part of the evolution has ‘hit’ the boundary, and we stop the calculation. This happens especially if $\beta$ is large and $\alpha$ is small, because in this case, the wave moves with speed $\sqrt{\beta}$, and may hit the boundary before $\alpha$ can lead to a divergence.

**Conjecture C.1.** We believe that this phenomenon might account for the idea that large $\beta$ regularises the PDE in cosmological simulations as discussed in [22]. But, as we show in section B the mathematical fact is that all solutions diverge in a finite time (unless they are 0). However, this will require a detailed study in a cosmological context.

## ORCID iDs

Jean-Pierre Eckmann https://orcid.org/0000-0001-6118-1077

Farbod Hassani https://orcid.org/0000-0003-2640-4460

Hatem Zaag https://orcid.org/0000-0002-1038-1201

## References

[1] Adamek J, Daverio D, Durrer R and Kunz M 2016 General relativity and cosmic structure formation *Nat. Phys.* **12** 346–9

[2] Ade P et al 2016 Planck 2015 results. XIII. Cosmological parameters *Astron. Astrophys.* **594** A13

[3] Aghamousa A et al 2016 The desi experiment part I: science, targeting, and survey design (arXiv:1611.00036)

[4] Alam S et al 2017 The clustering of galaxies in the completed SDSS-III Baryon oscillation spectroscopic survey: cosmological analysis of the DR12 galaxy sample *Mon. Not. R. Astron. Soc.* **470** 2617–52
[5] Alinhac S 1995 *Blowup for Nonlinear Hyperbolic Equations* (Progress in Nonlinear Differential Equations and Their Applications vol 17) (Birkhäuser Boston, Inc.)

[6] Amendola L et al 2018 Cosmology and fundamental physics with the Euclid satellite *Living Rev. Relativ.* **21** 2

[7] Armendariz-Picon C, Mukhanov V F and Steinhardt P J 2001 Essentials of k-essence Phys. Rev. D **63** 103510

[8] Azaiez A, Masmoudi N and Zaag H 2019 Blow-up rate for a semilinear wave equation with exponential nonlinearity in one space dimension *Partial Differential Equations Arising From Physics and Geometry (London Mathematical Society Lecture Note Series* vol 450) (Cambridge University Press) pp 1–32

[9] Babichev E 2016 Formation of caustics in k-essence and Horndeski theory *J. High Energy Phys.* JHEP04(2016)129

[10] Babichev E, Mukhanov V and Vikman A 2008 k-Essence, superluminal propagation, causality and emergent geometry *J. High Energy Phys.* JHEP02(2008)101

[11] Babichev E and Ramazanov S 2017 Caustic free completion of pressureless perfect fluid and k-essence *J. High Energy Phys.* JHEP08(2017)040

[12] Bonvin C, Caprini C and Durrer R 2006 A no-go theorem for k-essence dark energy Phys. Rev. Lett. **97** 081303

[13] Cheung C, Fitzpatrick A L, Kaplan J, Senatore L and Creminelli P 2008 The effective field theory of inflation *J. High Energy Phys.* JHEP03(2008)014

[14] Clifton T, Ferreira P G, Padilla A and Skordis C 2012 Modified gravity and cosmology Phys. Rep. **513** 1–189

[15] Cusin G, Lewandowski M and Vernizzi F 2018 Nonlinear effective theory of dark energy *J. Cosmol. Astropart. Phys.* JCAP04(2018)061

[16] Endlich S, Nicolis A, Porto R A and Wang J 2013 Dissipation in the effective field theory for hydrodynamics: first order effects Phys. Rev. D **88** 105001

[17] Goldberger W D and Rothstein I Z 2006 An effective field theory of gravity for extended objects Phys. Rev. D **73** 104029

[18] Gubitosi G, Piazza F and Vernizzi F 2013 The effective field theory of dark energy *J. Cosmol. Astropart. Phys.* **2013** 32

[19] Hairer E, Nørsett S P and Wanner G 1987 *Solving Ordinary Differential Equations. I: Nonstiff Problems (Springer Series in Computational Mathematics* vol 8) (Springer)

[20] Hassani F 2020 Characterizing the non-linear evolution of dark energy models *PhD Thesis* Geneva University, Department of Theoretical Physics

[21] Hassani F, Adamek J and Kunz M 2020 Clustering dark energy imprints on cosmological observables of the gravitational field *Mon. Not. R. Astron. Soc.* **500** 4514–29

[22] Hassani F, Adamek J, Kunz M, Shi P and Wittwer P 2022 A new instability in clustering dark energy?

[23] Hassani F, Adamek J, Kunz M and Vernizzi F 2019 *k*-evolution: a relativistic N-body code for clustering dark energy *J. Cosmol. Astropart. Phys.* JCAP1912(2018)011

[24] Hassani F, Shi P, Adamek J, Kunz M and Wittwer P 2022 New nonlinear instability for scalar fields Phys. Rev. D **105** L021304

[25] Hu B, Raveri M, Frusciante N and Silvestri A 2014 Effective field theory of cosmic acceleration: an implementation in CAMB Phys. Rev. D **89** 103530

[26] John F 1981 Blow-up for quasilinear wave equations in three space dimensions *Commun. Pure Appl. Math.* **34** 29–51

[27] Pich A 1998 Effective field theory: course *Les Houches Summer School in Theoretical Physics, Session 68: Probing the Standard Model of Particle Interactions* pp 949–1049

[28] Rammaha M A 1995 Upper bounds for the life span of solutions to systems of nonlinear wave equations in two and three space dimensions *Nonlinear Anal.* **25** 639–54

[29] Rammaha M A 1997 A note on a nonlinear wave equation in two and three space dimensions *Commun. PDE 22* 799–810

[30] Santos M G et al 2015 Cosmology from a SKA HI intensity mapping survey (arXiv:1501.03989)

[31] Sawicki I, Marra V and Valkenburg W 2013 Seeding supermassive black holes with a non-vortical dark-matter subcomponent Phys. Rev. D **88** 083520
[32] Scolnic D M et al 2018 The complete light-curve sample of spectroscopically confirmed SNe Ia from Pan-STARRS1 and cosmological constraints from the combined pantheon sample *Astrophys. J.* **859** 101

[33] Shi P et al On a second order in time Hamilton-Jacobi type equation issued from cosmology (in preparation)

[34] Shi P et al Scale-invariant solutions to a Hamilton-Jacobi type equation issued from cosmology (in preparation)

[35] Shi P et al Stability of a blow-up solution for a second order in time Hamilton-Jacobi type equation (private communication)

[36] Walcher C J 2019 4MOST Scientific Operations *Messenger* **175** 12–16

[37] Zumalacárregui M, Bellini E, Sawicki I, Lesgourgues J and Ferreira P G 2017 hi_class: Horndeski in the cosmic linear anisotropy solving system *J. Cosmol. Astropart. Phys.* JCAP08(2018)019