Signatures of Horndeski gravity on the Dark Matter Bispectrum

We present a detailed study of second-order matter perturbations for the general Horn- deski class of models. Being the most general scalar-tensor theory having second-order equations of motion, it includes many known gravity and dark energy theories and General Relativity with a cosmological constant as a specific case. This enables us to estimate the leading order dark matter bispectrum generated at late-times by gravitational instability. We parametrize the evolution of the first and second-order equations of motion as proposed by Bellini and Sawicki (2014), where the free functions of the theory are assumed to be proportional to the dark energy density. We show that it is unnatural to have large 10% ( 1%) deviations of the bispectrum introducing even larger ~ 30% (~ 5%) deviations in the linear growth rate. Considering that measurements of the linear growth rate have much higher signal-to-noise than bispectrum measurements, this indicates that for Horndeski models which reproduce the expansion history and the linear growth rate as predicted by GR the dark matter bispectrum kernel can be effectively modelled as the standard GR one. On the other hand, an observation of a large bispectrum deviation that can not be explained in terms of bias would imply either that the evolution of perturbations is strongly different than the evolution predicted by GR or that the theory of gravity is exotic (e.g., breaks the weak equivalence principle) and/or fine-tuned.


Introduction
One of the most challenging problems in modern cosmology is to understand the nature of the accelerated expansion of the Universe at late-times [1,2]. The simplest model that explains this behaviour is the Λ-cold dark matter one (ΛCDM), where gravity is described by General Relativity (GR) at all scales and the the cosmological constant Λ is used in order to have acceleration on cosmological scales. It is a simple model, with a minimal number of parameters that fits exceptionally well a suite of observations ranging from the early Universe to today and over 4 order of magnitude in scale [3,4]. However, the value of Λ appears to be too small to be explained by fundamental physics (e.g., [5]).
One popular alternative is to add a scalar degree of freedom to the Einstein equations, either by acting as an additional and exotic fluid of the system (as for Dark Energy (DE) models), or by modifying directly the laws of gravity (as for Modified Gravity (MG) models). In this paper we focus on the Horndeski theory introduced in [6][7][8]. The action reads where g µν is the metric, L m is the matter lagrangian, which contains just one single fluid composed by Dark Matter (DM) and baryonic matter, and L i are Here, K, G 3 , G 4 and G 5 are arbitrary functions of the scalar field φ and its canonical kinetic term X = −φ ;µ φ ;µ /2. The subscript X represents a derivative with respect to (w.r.t.) X, while a subscript φ would denote a derivative w.r.t. the scalar field φ. The action, Eq. (1.1), is the most general action for a single scalar field that has second-order equations of motion on any background and satisfies the weak equivalence principle, i.e., all matter species are coupled minimally and universally to the metric g µν . It encompasses many of the classical DE/MG models studied to explain the late-time cosmic acceleration: 1 quintessence [9,10], kinetic gravity braiding [11][12][13], galileons [14,15], f (R) [16]. We must note that with the Horndeski lagrangian it is not possible to describe models which break the weak equivalence principle, such as coupled quintessence [17] and non-universal disformal couplings [18][19][20], or Lorentz-invariance-violating models (e.g. Hořava-Lifshitz gravity [21,22]), or beyond Horndeski theories [23,24]. Measurements of the expansion history of the Universe are in excellent agreement with ΛCDM, e.g. [4]. Then, it is justified to assume it as our fiducial model for the background evolution. However, the same behaviour could be obtained assuming a particular DE/MG model with the requirement that the equation of state of the new degree of freedom is w (t) ≡ p(t) ρ(t) ≈ −1, where ρ (t) and p (t) are the background energy density and pressure of the DE/MG fluid respectively.
It is important to look at the growth of cosmological perturbations via gravitational instability as traced by the Large Scale Structure (LSS) of the Universe to learn about the nature of the cosmic acceleration. Indeed in ΛCDM, and in minimally coupled quintessence models, once the expansion history is given, the growth of structures is fully determined. This one-to-one correspondence is broken by the additional degree of freedom if we use a generic sub-class of the Horndeski action, Eq. (1.1). As explained in greater detail in the following section, it is therefore important to have a description that separates the expansion history from the evolution of the perturbations. This result can be achieved by using an Effective Field Theory approach, which has been developed for Horndeski models in [25][26][27][28][29][30][31].
In this paper we look at the evolution of matter perturbations at first and second orders by using the general Horndeski action Eq. (1.1). The first statistics of interest is the Power Spectrum (PS) P (t, k) defined as is the 3-dimensional Dirac delta, and . . . indicates the ensemble averaging over the possible configurations of the Universe. The PS encodes all the statistical information of a Gaussian random field. However, due to gravitational instability, which amplifies the initial small fluctuations, the late Universe is highly non-Gaussian. The lowest order statistics sensitive to the non-linearities is the bispectrum B t, k 1 , k 2 , k 3 , defined as where the Dirac delta imposes that only closed triangle configurations are to be considered. In GR the expansion history of the Universe determines the growth of the perturbations at all orders. Therefore, linear and non linear growth histories do not add qualitatively new information (beside shrinking the statistical error-bars). In this paper we want to show in which cases the DM bispectrum carries significantly new information about the growth of structures that is not included in the linear-order growth. With the data coming from current and forthcoming Large Scale Structure (LSS) surveys, such as BOSS [32], WiggleZ [33], DES [34] or Euclid [35], it will be possible to measure the PS and the bispectrum with unprecedented statistical precision (∼ 1%). The PS is a well known tried and tested statistical tool, a work-horse of any statistical analysis of LSS, which is extremely high signal to noise, and with well understood statistical properties. The linear growth rate can be measured in several independent ways (e.g., weak lensing, redshift space distortions etc.). On the other hand (mildly) non-linear growth is much harder to measure. The most widespread tool is the bispectrum, but it has much lower signal to noise and the information is spread out over several configurations. It is complicated to model and time consuming even in a standard GR-ΛCDM model (eg. [36,37]). Then, we want to know to what extent the modelling developed for this fiducial case (i.e. GR-ΛCDM) is more general and can be applied more widely.
The general result for the leading order DM bispectrum assuming Gaussian initial conditions reads (e.g., see [38] for a review) B t, k 1 , k 2 , k 3 = F 2 t, k 1 , k 2 P (t, k 1 ) P (t, k 2 ) + cyc., (1.8) where cyc. contains the possible permutations w.r.t k i and F 2 is a model dependent kernel which was derived for Horndeski models in [39]. We shall show the expression for F 2 in the next section.
Here, it is important to note that the functional form of the bispectrum in terms of the kernel and the linear PS is identical to the standard one. As a result of Eq. (1.8), the tree-level DM bispectrum can be separated in terms of linear quantities, i.e. P t, k , and second-order quantities, i.e. F 2 t, k 1 , k 2 . The DM bispectrum has been investigated using second-order perturbation theory or simulations in simple DE/MG theories as Brans-Dicke/f (R) or cubic galileon [40][41][42][43][44]. The result is that, in these Horndeski sub-models, the DM bispectrum kernel appears to be close to the standard one within deviations at the percent level, and then it can be approximated as F 2 F 2GR . The main goal of this paper is to investigate if it is possible to generalise this approximation for a wider class of Horndeski-type DE/MG models. In these cases, it would be possible to avoid the secondorder analysis, and concentrate just on the evolution of linear perturbations by applying directly the bispectrum constraints obtained assuming ΛCDM/GR e.g., [36,37].
The paper is organised as follows. In Section 2 we describe the approach we used and show the main equations, while Appendix A is reserved to useful formulas that are too long to fit in the main text. In Section 3 we present and discuss the main results, and in Section 4 we draw our conclusions.

Approach
We assume that the Universe at large enough scales is well described by small perturbations on top of a flat Friedmann-Robertson-Walker (FRW) metric. We consider scalar perturbations in the Newtonian gauge, and adopt the notation of [45]. Then, up to second order the line element reads where for brevity we have neglected the superscript (1) for linear perturbations. Here Ψ and Φ are the metric perturbations and in this gauge they can be interpreted as the gauge invariant Bardeen's potentials [46]. The perturbation of the scalar field is rewritten as v X ≡ −δφ/φ, where φ and δφ are the background and the perturbation of the scalar in Eq. (1.1) respectively. In the following we will show the equations for the metric perturbations, after having decoupled them from the scalar field perturbations.
In [31], see also [26][27][28] for an equivalent approach, the authors describe the evolution of linear perturbations in general scalar-tensor theories belonging to the Horndeski class of models, Eq. (1.1). The idea is to identify the minimum number of operators that fully specify the linear evolution of those models (see [31] for details). These operators are independent of each other, i.e., they can be parametrized independently. One of the advantages of this description is that it separates the effects of the background expansion from the effects of the perturbations. Then, the total amount of cosmological information can be compressed into five functions of time plus one constant, namely where Ω m0 is the value of the matter density fraction today, H (t) is the Hubble function and α i (t) represent the linear freedom of the Horndeski class of models. α K , the kineticity, is the most standard kinetic term present in simple DE models as quintessence or k-essence. α B is called braiding, it represents a mixing between the kinetic terms of the metric and the scalar. α M , the Planck mass run rate, describes systems in which the Planck mass is varying with time. α T is the tensor speed excess, and it directly modifies the speed of tensors. Both α M and α T produce anisotropic stress between the gravitational potentials Ψ and Φ, and they are signatures of modifications in the evolution of the gravitational waves [47]. Generalizing the results of [31] to second-order in perturbation theory we find four other functions of time which are free to vary in the general Horndeski theory. However, for our purposes we just need two of them, since the other two appear in front of terms that are subdominant on sub-horizon scales (k 2 a 2 H 2 ). Then, the total freedom one can have is represented by where the new functions α 5 and α 4 represent the second-order freedom of this theory and are defined in Appendix A. From their definitions, α 5 and α 4 can be considered as second-order corrections to α M and α T , then we can expect that realistic DE/MG theories will have Models that do not satisfy the second inequality have a strong non-minimal coupling, and they are expected to drastically modify the evolution of first-order scalar and tensor perturbations. Then, they can be excluded by observations that rely on linear theory, without the need of second-order analysis. The first inequality in Eq. (2.4) describes the fact that α 5 and α 4 are built with the same Horndeski functions as α M and α T , i.e., G 4 and G 5 (see Eq. (A.1-A.2) and Ref. [31] for comparison).
In principle it is possible to create a model in which the second-order functions are larger than the first-order ones, but it would require fine tuning. We consider Eq. (2.4) as a naturalness condition, which we will however not impose a priori. While our analysis below is general and does not assume this hierarchy, it is reasonable to concentrate attention on the "natural" regime of Eq. (2.4).
As specified above, with this approach we separate the contribution of the background from the perturbations. Indeed, the α i functions are constructed to be independent form the background in the general Horndeski class of models. Then, we can assume a flat ΛCDM cosmology for the evolution of H (t) and fix Ω m0 = 0.31 [4], this being our fiducial model. The Friedmann and the conservation equations read whereΩ m ≡ ρm 3M 2 * H 2 is the DM fractional density, and M 2 * (t) the effective Planck mass of the theory [31]. Here, we do not need to show explicitly the contribution of the energy density of the DE fluid, since it has been eliminated using the Friedmann constraint, i.e.Ω Λ = 1 −Ω m . As shown in [48], if the sound speed of the DE/MG fluid is sufficiently close to the speed of light ( 0.1c) we can also assume the Quasi-Static approximation, in which time derivatives are considered to be sub-dominant w.r.t. space derivatives (i.e. a 2Φ k 2 Φ). In this regime it is possible to prove that α K (t) does not enter the equations of motion. Finally, we choose a parametrization for the remaining α i , precisely the one suggested in [31] where c i are arbitrary constant. With this choice we ensure the standard evolution of perturbations during radiation and matter domination, while we allow for modifications when DE starts dominating. This parametrization reproduces exactly the imperfect fluid behaviour on its tracking solution introduced in [11], and it can describe general models in which the deviations w.r.t. ΛCDM are smooth.
Note that the naturalness condition Eq. (2.4) then becomes: Before evaluating the linear and second-order evolution, we will adopt a further simplification. It has been noticed in [31] the existence of the braiding scale This represents the only new scale in the linear equations beyond the usual Jeans length. Models that evolve outside this scale have a perfect fluid structure, while clustering DE on linear scales is possible well inside the braiding scale. In the following we focus only on two classes of models, the ones that have sub-braiding evolution (α 2 B α K ), and the ones with super-braiding evolution (α 2 B α K ). As we shall see, in this limit we ensure that the linear-order DM density evolution is scale independent. In the intermediate regime, where α 2 B ∼ α K , we expect the growth of perturbations to be scale dependent. This would modify the shape of the PS, and these models in this regime could be already studied at this level without including the bispectrum. In fact the shape of the (matter) PS at linear and mildly non-linear scales is a high signal to noise and relatively robust quantity to measure. It is much less affected by e.g., non-linear galaxy bias than measurements of e.g., the growth rate of structures.

Linear-order and second-order evolution
The linearized Einstein equations together with the equation of motion for the scalar field, as usual, provide the Poisson equations for the gravitational potentials, e.g. [49][50][51], where δ m ≡ δρ m /ρ m is the density contrast, while G Ψ (t) and G Φ (t) are the effective Newton's constants. In GR we have G Ψ = G Φ = 1, while their definition in our case is given in Appendix A. It is interesting to note that G Ψ and G Φ are scale independent in the limits we are considering, i.e. sub-braiding and super-braiding. Eqs. (A.3-A.4) are presented in Appendix A in the sub-braiding case. Their expression in the opposite limit is identical to Eqs. (A.3-A.4) imposing α B → 0. Note that, in order to obtain the previous equations, we already decoupled the metric from the scalar field perturbations v X . Its effects are fully encoded in G Ψ and G Φ . In addition, the conservation of the stress-energy tensor gives the linear evolution of the DM density contrast, which reads Following the same procedure described in [40], we can get the evolution of the second-order density fluctuations. The result is an equation which is identical to Eq. (2.13), with the addition of a non-vanishing source (S δ ) composed by products of first-order quantities (e.g. δ m δ m ), 14) The exact form ofS δ is not important for the purpose of this paper. More details can be found in [39], where the bispectrum for the Horndeski class of models was first studied. The solution of Eq. (2.14) is given by where the kernel F 2 coincides with the kernel shown in Eq. (1.8) and encodes the second-order modifications w.r.t. the usual Newtonian one [38]. Here we have identified the evolution of the second-order DM perturbations with a particular solution of Eq. (2.14), neglecting the homogeneous part. The homogeneous solution contains both a possible primordial non-Gaussianity and a non-primordial contribution. Since in our analysis the evolution of cosmological perturbations is standard up to redshift z ∼ 1 by construction, the non-primordial contribution remains subdominant on mildly non-linear scales, see e.g. [52,53]. In order to select the leading-order contribution on sub-horizon scales, we assume Gaussian initial conditions and we perform an expansion in terms of ≡ aH ki in F 2 , where k i stands for k, q 1 and q 2 . One could argue that with this method we are neglecting terms proportional to aH qi , which are important in the integral when q i → 0. However, these terms represent wavelengths that are super-horizon, they can be considered as a backreaction and reabsorbed into the background evolution. Then, the kernel F 2 read [39] The full definition for the function C (t) is given in Appendix A in Eqs. (A.5-A.6), and it represents the only effect of Horndeski type of DE/MG at the leading order on mildly non-linear scales. After some algebra, it is possible to demonstrate that our results, Eq. (2.16) and Eqs. (A.5-A.6), reproduce the results of [39]. With our notation, the standard value for this function in an Einstein de Sitter Universe is C (t) = 34/21. The dependence of this value on the cosmology is very weak. Indeed in GR, for ΛCDM and minimally coupled models as quintessence, the corrections to the standard value are negligible [54][55][56][57][58]. Note that the denominator of Eq. (A.6) does not produce divergences. Indeed, after some algebra it is possible to demonstrate that, if we tune the numerator will suppress these divergences with factors O 2 . Even if Eqs. (A.5-A.6) appear long and complicated, it is possible to identify three main objects. The first one depends only on first-order quantities, the second depends on (α 5 ,α 5 ), and the third on (α 4 ,α 4 ). Then, using our parametrization for the Horndeski functions α i , Eq. (2.7), we can rewrite C(t) in the expression for the second order gravitational kernel Eq. (2.16) as

Results
In this section we show the main results of this paper. At linear-order we focus on the effective Newton's constant G Ψ , Eq. (A.3), the growth rate f ≡ d ln δm d ln a and the slip parameterη ≡ 2Ψ Ψ+Φ . The Newton's constant measures the strength of the gravitational force, and modifies directly the matter PS. The growth rate tracks the growth of linear perturbations. The slip parameter quantifies the anisotropic stress between the two gravitational potentials Ψ and Φ, and can be a signature of non-minimal coupling between the metric and the scalar field [47]. It can be measured by comparing weak lensing with redshift space distortions and the estimated errors with data from future surveys such as Euclid are ∼ 10% in case it is not scale dependent [59]. All these quantities are indications of deviations w.r.t. the prediction of a pure ΛCDM, which we consider as our reference model. From Eqs. (2.11-2.12-2.13) it is possible to note that they are independent each other. With present surveys f is measured with ∼ 6% precision [60,61] and G Ψ is measured (although in a model-dependent way) with ∼ 10% precision [62] butη is currently poorly constrained (∼ 100% error) e.g., [63].
Future data will increase this precision to roughly one order of magnitude. Then, for viable models we impose that the most stringent condition is satisfied. Should future surveys not find any deviations from standard cosmology, this condition will become much stronger, with deviations allowed only below the 1%. We explore the parameter space {c B , c M , c T } and we ensure that the models we pick do not suffer from gradient and ghost instabilities, both on the scalar and the tensor sectors (see [31] for details). For understanding the behaviour of the second-order bispectrum we study the quantities {A 0 , A 5 , A 4 } evaluated today and defined in Eq. (2.17) as a function of the free parameters. Indeed they are the only modifications to the kernel of second-order DM matter perturbations. We choose to estimate them today in order to maximise the effects of DE/MG on the bispectrum. With current surveys it is possible to estimate the bispectrum with a ∼ 10% precision, due to statistics and systematic uncertainties [36,37]. With future surveys it will be possible to reduce the statistical uncertainties, and, optimistically reach a ∼ 1% precision and accuracy when averaged over all bispectrum configurations. Then, in the following we consider bispectrum kernels modified by ∼ 1% or less as kernels that are close to the standard result. On the contrary, models that can show interesting signatures of DE/MG on mildly non-linear scales with the next generation of experiments are those models with modifications of the bispectrum kernel 10%.

Small braiding
In this section we consider models with small braiding, i.e. α B → 0. This ensures that the scalar degree of freedom evolves always on super-braiding scales. In this case the scalar field has a perfect fluid structure, with the possibility of having a non-minimal coupling given by {c M , c T }. However, after some algebra one can show that the dependence on c T vanishes both at first and at second orders. This result is a consequence of choosing ΛCDM for the background and α B = 0, it is therefore independent of the parametrization for the α i functions we chose. Then, at first-order we are left with just one parameter to vary, i.e. c M . In Fig. (1) we show the evolution of the effective Newton's constant G Ψ top panel), the growth rate f (central panel) and the slip parameterη (bottom panel) as a function of the scale factor. Increasing the value of |c M | is equivalent to increase the deviations w.r.t. the standard values G Ψ = 1, f = f ΛCDM andη = 1. Note that the extreme models, c M = ±0.5, are already disfavoured from current data as they induce effects on observable quantities that are larger than the limit of Eq. (3.1). Thus c M 0.2 but for forecasted errors from future surveys, should there be no deviations from standard cosmology, then c M 0.1. In Fig. (2) we show the second-order interesting quantities. In particular, in the left panel we plot the value of A 0 (a = 1), Eq. (2.17), normalized by its value in a ΛCDM universe as a function of the parameter c M . One can see that for small values of c M the standard limit is recovered, while one can reach deviations of the order of 1% when c M = 1. However such large value for c M induces uncomfortably large changes on first-order quantities, see discussion around Fig. (1). We can conclude that A 0 in Eq. (2.17) is effectively standard. Therefore the only models that have a chance to significantly change the bispectrum are those where α 5 , α 4 are non-zero.
In the right panel we plot the values of A 0 /A 5 and A 0 /A 4 , Eq. (2.17), calculated today as a function of c M . In the relevant range for c M , A 5 is 10 to 20 times smaller than A 0 and |A 4 | is 20 to 40 times smaller than A 0 . The ratio A 0 /A i therefore indicates how large the parameters {c 5 , c 4 } have to be in order to modify the bispectrum kernel, Eq. (2.16). The result is that, in the limit c M → 0 or when it is negative we need {c 5 , c 4 } {1, 2} in order to have a kernel modified by 10% and {c 5 , c 4 } > {0.1, 0.2} to have a kernel modified by more than 1%. When c M → 1, we need even larger values of {c 5 , c 4 } to have the same deviation.
In conclusion, for models with negligible braiding, it is necessary to require c M 1 in order to have the evolution of the linear perturbations not in tension with current constraints, Eq. (3.1) [60,61]. As shown in Fig. 2, at second-order the kernel, Eq. (2.16), is weakly dependent on the value of c M and its modifications are mostly due to the parameters c 5 , c 4 . It is possible to note that the kernel is modified by less than 1% w.r.t. the standard one if we impose c 5 , c 4 0.1. Differences of order 10% can arise for models with c 5 , c 4 1. Then, an interesting bispectrum kernel modification can happen only if the parameters of the DE/MG model we consider do not satisfy the expected natural hierarchy, Eq. (2.9).

Large braiding
Here we want to study models that evolve on sub-braiding scales. For these models we expect a more interesting physics, since the braiding causes the DE to cluster on linear scales. On the other hand, the equations are longer and there are no reduction on the linear-order free parameters. We then explore the parameter space {c B , c M , c T } using two grids of 30 points in each direction. As we shall see, the dependence of the interesting first and second-order functions on the parameters is smooth, then the choice of using just 30 points is justified. If we require that the linear order quantity f should not be modified by more than 6% w.r.t. standard gravity, Eq. (3.1), we find that the second-order parameters should be at least c 5 , c 4 2(1) in order to have an effect of at least 10%(3%) on the kernel Eq. (2.16).
If we want to keep c 5 , c 4 small ( 1) and have large modifications on the kernel ∼ 10% (3%) we have to allow for deviations of ∼ 30% (10%) at linear-order. However, to find signatures of these kind of models, a higher-order correlations (and thus second-order) analysis is not needed, since any signature of DE/MG could be seen easily using observables that rely on linear theory. There are few exceptions to this result in the region of the parameter space where c T ∼ −0.5. Here we can keep c 5 , c 4 small ( 1), have deviations on the linear growth rate 6% and have large deviations on the bispectrum kernel (∼ 10%). However, such c T values imply a speed of tensors significantly smaller than the speed of light which can be ruled out by observations of ultra-high energy cosmic rays [64][65][66].
In Fig. (3) we show the linear evolution for few representative models. The left panels refer to the effective Newton's constant, in the central panels we plot the linear growth rate and in the right panels we show the time evolution of the slip parameter. In the top panels we present the behaviour of a minimally coupled theory (c M = c T = 0). We vary c B in the region where ghost and gradient instabilities are avoided, i.e. c B > 0. It is possible to note that c B enhances the effective Newton's constant and the linear growth. This is expected, and it is the signal of a stronger gravity, which should be screened on small scales to satisfy solar system constraints. It is also expected the absence of anisotropic stress (top right panel), sinceη = 1 only in non-minimally coupled theories. In the central panels we fix c B = 1 and c T = 0, letting c M free to vary. In this case we can compensate the effect of the braiding by choosing a negative value for c M . By fixing c M −0.4 we can recover nearly the standard predictions for G Ψ , but the growth rate f appears slightly suppressed. This is because in our description c M affects directly the evolution of the matter density through Eq. (2.6). Finally, in the bottom panels we show the evolution of G Ψ , f andη for different values of c T . We set c B = 1 and c M = 0. Note that we have chosen extreme values for c T , and c T = −0.5 could be ruled out by the Cherenkov radiation emission argument [64][65][66]. Despite this fact, the parameter c T seems to affect the evolution of linear perturbations less than the others.
As in the previous section, to understand the bispectrum modifications we focus on the functions {A 0 , A 5 , A 4 }. In Fig. (4) we show the value assumed today by A 0 (left panels) and A 5 and A 4 (right panels) as a function of c B (top panels), c M (central panels) and c T (bottom panels). For all the models under consideration A 0 deviates from its standard value by less than 1%. Therefore also in this case this function cannot be responsible for large deviations in the bispectrum and only models with non-zero α 5 and α 4 can modify it. Also the functions A 5 and A 4 appear rather suppressed compared to A 0 . Indeed, |A 0 /A i | > 10 for every model under consideration. Moreover, for some value of c M (central right panel) and c T (bottom right panel) the curves apparently diverge. For these models |A 5 , A 4 | A 0 , which indicates that these terms can not realistically modify the bispectrum kernel. In the cases where |A 0 /A i | is minimized we still need to choose c 5 , c 4 1 in order to achieve deviations in the kernel of at least 10% . One could argue that these deviation can be reached without imposing c 5 , c 4 1 just by increasing the braiding c B . However, comparing these plots with the ones in Fig. (3) we note that the limit c B 1 produces large modifications of the first-order evolution. Even if we introduce a negative c M to compensate the linear-order deviations we expect for A 5 and A 4 to be further suppressed (see central-right panel of Fig. (4)).
In conclusion, even for general models with large braiding, at second-order in perturbation theory we can not produce large deviations in the bispectrum kernel for models that respect the expected hierarchy, Eq. (2.4) without generating even larger effects at first order.

Bispectrum shapes
In this section we show how the shape of the bispectrum is modified for Horndeski-type DE/MG models. This is useful in order to estimate what are the configurations for which the modifications of the bispectrum kernel, if any, are maximised, i.e., the easiest to detect. In the previous sections we considered the evolution of the function C (t) and we showed that it is unnatural to have large modifications in C (t) evaluated today. This is the only modification that appears in the bispectrum kernel, Eq. (2.16), and, its magnitude at any particular time, is just a number. Therefore here in full generality we do not consider particular models, just the value of C (z obs ), where z obs is the redshift at which we want to measure the bispectrum. We perturb C (z obs ) from it standard value (C ΛCDM 1.62) by fixed percentages (1%, 2%, 5%, 10%). Recall that the deviations in the bispectrum kernel are maximised today and vanish at early times by construction. However, forthcoming and future surveys do not measure the bispectrum today, but at higher redshifts (z ∼ 0.5 − 1). Then, our main result, i.e. the modifications of the bispectrum for natural models are 1% today, can be considered as an upper bound on the modifications of the bispectrum that a generic survey will measure if gravity satisfies our naturalness condition.
, (3.2) which removes most of the cosmology and scale dependence. Here, the PS can be written as where δ + (t) is the growing solution of Eq. (2.13) and encodes all the time dependence of the PS. T (k) is the transfer function, which we take as standard using the fitting formula given in [69]. This is consistent with our parametrization, since we only allow for modifications at late-times while we assume that the Universe behaves as standard during the epochs dominated by radiation and matter. Finally n s is the scalar spectral index, and we fix it at n s = 0.96 [4]. It is important to keep in mind that, since we are looking at models either sub or super-braiding, the growth of perturbations is not scale dependent. As a consequence, one can see using Eqs.
In the left panels we show A0/AΛCDM evaluated today, while in the right panels we show A0/A5 (red line) and A0/A4 (blue line) evaluated today. In the top panels we fix cM = cT = 0 and we vary cB. In the central panels we fix cT = 0, cB = 1 and we vary cM. In the bottom panels we fix cM = 0, cB = 1 and we vary cT. In the left panels it is possible to note that the standard term of the kernel, Eq. (2.16), have deviations w.r.t. the GR one 1% for every case under consideration. The right panels show that c5 and c4 should be 2 in order to have an interesting effect on the bispectrum kernel (∼ 10%). In the right central and bottom panels it is possible to note divergencies in the curves for cM = −0.5 and cT = 0.6 respectively. Indeed for these models |A5, A4| A0, which is a sign that these terms can not realistically modify the bispectrum kernel.
i.e. ∆C ≡ C−CΛCDM CΛCDM . It is important to note that in the left panel the limit θ 12 → π involve extremely large scales (k i → 0). Then, since we are assuming to be well inside or outside the braiding scale, Eq. (2.10), the behaviour of the curves in these limits could probably be corrected for real models.
In every panel the differences in the reduced bispectrum are maximised at θ 12 0.6π, where they are three times bigger than the differences in C, i.e. ∆Q ∼ 3∆C. For all the other configurations the deviations from the standard bispectrum are suppressed. This can be connected with the findings of the previous sections, where we considered modifications in C of ∼ 3% as large ones. Indeed, the observable is Q in particular configurations, and not C. Then, we have been conservative in considering the maximum deviations as the deviations that could be detected by forthcoming and

Conclusions
We have presented a systematic analysis of the DM bispectrum generated at late-times by gravitational instability in the general Horndeski class of models. We assumed a ΛCDM evolution for the Hubble function H (t) and we parametrized the perturbations. In this way we were able to separate the effects of the expansion history of the Universe from the first and second-order effects.
We focused on two class of models, the ones with large braiding (α 2 B α K ) and the ones with small braiding (α 2 B α K ). In the intermediate regime, where α 2 B ∼ α K , we expect the linear evolution of perturbations to be scale-dependent. Then, any deviation from standard gravity can be seen looking, e.g., at the PS without the need of a bispectrum analysis.
In the small braiding case we noted that the only parameter that is left free to vary at linearorder is the Planck mass run-rate (c M ). Then, we showed that, even if we allow for large deviations at linear-order ∼ 10% (∼ 1%) we can not have large modifications at second-order ∼ 3% (∼ 1%) if c 5 , c 4 c M , c T . Since the parameters c 5 and c 4 represent the second-order contribution of non minimal coupling between the metric g µν and the scalar field φ, they can be seen as corrections of the linear parameters c M and c T . This argument leads to Eq. (2.4), which we consider to be the natural behaviour of realistic DE/MG models. The large braiding case is more complicated. We then explored the parameter space for realistic values of {c B , c M , c T }. We found that, assuming that the linear growth f is measured with a 6% precision (or better), there can not be significant deviations in the bispectrum if c 5 , c 4 c M , c T . The only way of having large modifications in the bispectrum kernel (∼ 3%) with c 5 , c 4 c M , c T , is to allow for unrealistic deviations on the growth rate f (∼ 30%). These models are expected to be excluded by linear observations, without the need of the analysis of the bispectrum.
Large and observable modifications of the DM bispectrum kernel are still possible in the Horndeski class of models. However, they can not come from the function A 0 in Eq. (2.17), since it can be modified by 1% for all the models under consideration. Indeed, only models with with large α 5 and α 4 , i.e. α 5 , α 4 α M , α T , have a chance to modify it. Note that in popular MG models such as as Brans-Dicke/f(R), kinetic gravity braiding, cubic galileons these functions are zero. It is therefore not surprising that no significant signatures of these models in the bispectrum were found in these cases. Large α 5 , α 4 imply a large non-minimal coupling between the curvature and the derivatives of the scalar field. 3 This can be achieved in, e.g., quartic and quintic galileons at the price of having large α M and α T .
In conclusion, we have shown that it is not possible to have large modifications in the DM bispectrum kernel for models following our parameterisation of the free Horndeski functions and that satisfy our naturalness condition, Eq. (2.4), without introducing even larger changes in firstorder quantities. Even if in these models there is no extra qualitatively new information in the DM bispectrum, we are certainly not advocating not to use it. Our findings imply that the kernel, Eq. (2.16), can be modelled as the standard GR one and applied more generally. Such bispectrum analysis is still useful to: 1) have new information on real word effects, e.g., bias parameters, to remove degeneracies; 2) decrease statistical errors; 3) offer a powerful consistency check.
Eq. (2.4) can be seen as a prescription, under which we can consider the DM bispectrum kernel as standard. Then, the observation of a large bispectrum deviation that can not be explained in terms of bias would imply either that the linear evolution of perturbations is strongly different than the evolution predicted by GR or that the theory of gravity is exotic and/or fine-tuned.
The modified Newton's constants in Eqs. (2.11-2.12) are At second-order, the function C (t) that modifies the bispectrum kernel, Eq. (2.16), reads where δ g and δ d are the growing and the decaying modes of Eq. (2.13) respectively, while W (t) ≡ δ g (t) δ d (t) − δ g δ d (t) (t) is the Wronskian. Finally c (t) is defined as