The phenomenology of beyond Horndeski gravity
Dina Traykova Emilio Bellini Pedro G. Ferreira





































a]Dina Traykova,
a]Emilio Bellini,
a]Pedro G. Ferreira
[a]Astrophysics, University of Oxford, DWB, Keble Road, Oxford OX1 3RH, UK
dina.traykova@physics.ox.ac.uk
emilio.bellini@physics.ox.ac.uk
pedro.ferreira@physics.ox.ac.uk


We study the phenomenology of the beyond Horndeski class of scalar-tensor theories of gravity, which on cosmological scales can be characterised in terms of one extra function of time, ,  as well as the usual four Horndeski set of free functions.
We show that  can be directly related to the damping of the matter power spectrum on both large and small scales. 
We also find that the temperature power spectrum of the cosmic microwave background (CMB) is enhanced  at low multipoles and the lensing potential is decreased, as a function of .  
For a particular choice of time dependence we find constraints on  of order O(1) using measurements of the temperature and polarisation of the CMB, as well as the lensing potential derived from it, combined with large scale structure data.
We find that redshift space distortion measurements can play a significant role in constraining these theories.
Finally, we comment on the implications of the recent observation of an electromagnetic counterpart to a gravitational wave signal;
we find that these measurements reduce the number of free parameters of the model but do not significantly change the constraints on the remaining parameters.





The phenomenology of beyond Horndeski gravity
    [
    Received January 31, 2019; accepted February 25, 2019
=========================================================







§ INTRODUCTION



In the last two decades, there has been tremendous progress in testing gravity on cosmological scales.
We can now measure the expansion history, the Cosmic Microwave Background (CMB) and the Large Scale Structure (LSS) of the universe with incredible precision. 
Distance measurements from supernovae <cit.>, combined with measurements of LSS clustering <cit.> and the CMB <cit.> have confirmed that, at recent times, the universe is undergoing a phase of accelerated expansion. 
The simplest model that can explain these measurements  is ΛCDM, in which the energy density content of the Universe is dominated by an additional stress-energy component in the form of a cosmological constant (Λ) that causes the observed acceleration. 

The standard model of cosmology assumes that the theory of gravity is General Relativity (GR); indeed the current assumption is that GR  works incredibly well over 15 orders of magnitude in length scale <cit.>.
However, the large amount of fine-tuning required to explain the value of the cosmological constant has led to the exploration of alternative models, from Dark Energy (DE) models, where an additional dynamical degree of freedom (d.o.f.) replaces Λ, to theories that directly modify the laws of gravity, i.e. Modified Gravity (MG) models <cit.>. 
A particular class of DE/MG models are the  scalar-tensor theories, where the additional d.o.f. is represented by a scalar field.

The Horndeski class of theories <cit.>  has received a lot of attention as it is given by a general scalar-tensor action that leads to at most second-order equations of motion on any background and also satisfies the weak equivalence principle, i.e. all matter (except from the scalar, ϕ) is coupled minimally and universally to the same metric.
In these theories, there is one extra propagating degree of freedom, ϕ, in addition to the standard two tensor modes of GR.
It includes, as special cases, a wide variety of dark energy and modified gravity models, such as quintessence <cit.>, f(R), Brans-Dicke <cit.>, kinetic gravity braiding <cit.>, covariant galileons <cit.>, Chameleons <cit.> and others (for a review see <cit.>).

It has been shown that Horndeski is not as general as previously believed and can be further extended to the so-called beyond Horndeski theories <cit.>.
While in general beyond Horndeski have higher (than second) order equations of motion, this does not necessarily lead to Ostrogradski instabilities <cit.>.
Indeed, it is possible to show that the time-derivatives remain second-order in the equations of motion for the propagating degrees of freedom for any ADM decomposition <cit.>.
As expected, beyond Horndeski theories lead to new interesting phenomenological properties that are not described by the Horndeski class <cit.>.
For example, a  recent  study <cit.> shows that the scalar field in beyond Horndeski theories can significantly modify the speed of sound in the atmosphere of the Earth.
The question of whether this set of theories is viable for explaining dark energy has been studied extensively in <cit.>;
self-acceleration solutions exist that can lead to a de Sitter universe at late times without the need of a cosmological constant.

The aim of this work is to explore the phenomenology of beyond Horndeski models and constrain the free parameters of the theory using the latest CMB and large scale structure data.
Such fits to data have been performed previously on specific modified gravity models <cit.> or specific data sets <cit.>.
More generally this was done for the Horndeski models considering a collection of latest CMB and large scale structure data in <cit.>.

We use the formalism introduced in <cit.> for Horndeski and used in <cit.> and extended to beyond Horndeski in <cit.> (see also <cit.> for a different approach to the linear perturbation theory in the Horndeski and beyond Horndeski Lagrangian).
In this approach the evolution of linear perturbations is described fully by the fraction of matter density today, the Hubble parameter H(t) and five other functions of time.
We have chosen the parametrisation for these five functions proposed in <cit.>, in which these are proportional to the fractional energy density of dark energy, Ω_ DE(t).
In this way the effective freedom of the model is reduced to just five free parameters.

The structure of this work is as follows. In Section <ref> we review the theoretical properties of the beyond Horndeski action and the equations of motion that it leads to.
Then in Section <ref> we discuss the phenomenology of this model.
We present our results for the constraints on the parameters of this model with data in Section <ref>. Finally we discuss our findings in Section <ref>.




§ THEORETICAL FRAMEWORK

In this section we discuss the general properties of the beyond Horndeski class of theories, present its Lagrangian, background equations and stability.
We then briefly discuss the recent restrictions arising from GW170817 and finally present the linear perturbation equations for the model. The action for beyond Horndeski reads <cit.>,

    S[g_μν,ϕ]=∫d^4x √(-g)[∑_i=2^51/8π G_N L_i(g_μν,ϕ) +ℒ_m(g_μν,ψ_M)] ,

where g_μν is the metric, g is its determinant and  ℒ_m generically describes all matter species, which we consider to be minimally and universally coupled to the metric g_μν. The sum is over four Lagrangians defined as,

    L_2 =      G_2(ϕ,  X) ,
     L_3 =      -G_3(ϕ,  X)ϕ ,
     L_4 =      G_4(ϕ,  X)R+G_4X(ϕ,  X)[(ϕ)^2-ϕ_;μνϕ^;μν]
        + F_4(ϕ,  X)ϵ^μνρ_       σϵ^μ^'ν^'ρ^'σϕ_;μϕ_;μ^'ϕ_;νν^'ϕ_;ρρ^' ,
     L_5 =      G_5(ϕ,  X)G_μνϕ^;μν
        - 1/6G_5X(ϕ,  X)[(ϕ)^3+2ϕ_;μ^νϕ_;ν^αϕ_;α^μ-3ϕ_;μνϕ^;μνϕ]
        + F_5(ϕ,  X)ϵ^μνρσϵ^μ^'ν^'ρ^'σ^'ϕ_;μϕ_;μ^'ϕ_;νν^'ϕ_;ρρ^'ϕ_;σσ^' ,

where ϕ_;μ≡∇_μϕ, ϕ_;μν≡∇_μ∇_νϕ, G_i(ϕ,  X) and F_i(ϕ,  X) are arbitrary functions of the scalar field ϕ and its canonical kinetic term 2X≡ - g^μνϕ_μϕ_ν.
The subscripts ϕ and X denote partial derivatives (G_iX≡∂ G_i/∂ X), and ϵ^μνρσ is the antisymmetric Levi-Civita tensor.

As suggested by <cit.> and shown in <cit.>, this action as it is written can lead to an extra propagating degree of freedom.
However, this issue can be avoided by imposing the following condition <cit.> ,

    2XG_5XF_4 = -3F_5[G_4+4XG_4X+XG_5ϕ] .

This is visible in the context of beyond Horndeski when non-linear terms are considered. However, in this work we only consider linear perturbations and do not need to impose this condition in obtaining our results.

Starting from this action, one has to choose the functional form of the G_i (and F_i), initial conditions for the scalar field, and get predictions solving for the background and the perturbations. 
However, this can be expensive and inefficient if the objective is to test the robustness of the standard cosmological model just by detecting deviations from it.
A more practical approach for this purpose is to use parametrized frameworks which were built to compress the information of a general class of models into few functions of time.
In particular we are going to use the framework introduced in <cit.>, where it was shown that the expansion history of the universe and the linear perturbations can be described completely by six functions of time: one responsible for the background evolution – H(t) or w(t) – and the remaining affecting only the perturbations.
It is possible to show that for a general theory that can be described as a sub-class of beyond Horndeski these functions are completely independent of the background and each other, which means that they represent the minimal set of functions that can describe all models within this set of theories.
Here we describe their physical meaning and definitions in terms of the G_i and F_i functions:


  * The kineticity, , is a generalisation of the standard kinetic term of simple DE models.
On its own this parameter describes some of the simpler DE models, with the only modification to ΛCDM being the addition of an extra fluid (e.g. quintessence).
It is defined as,

    = 
    4 H ϕ̇ X[3 G_5X + X(7 G_5X X - 120 F_5 + 2X (G_5X X X - 66 F_5X - 12 X F_5X X))]
       + 2 X/H^2[K_X - 2 G_3ϕ + + 2 X (K_X X - G_3ϕ X)]
       + 12 ϕ̇ X/H[G_3X - 3 G_4ϕ X + X(G_3X X - 2 G_4ϕ X X)]
       + 12X [G_4X - G_5ϕ + X( 8 G_4X X - 5 G_5ϕ X + 2X(2 G_4X X X - G_5ϕ X X))]
       + 12X[X(24 F_4 + 2X( 18 F_4X + 4 X  F_4X X))] .



  * The braiding, ⟨$|, describes the mixing of the scalar field and metric kinetic terms.
It contributes to the kinetic energy of the scalar, partially sourced through its mixing with the metric.
In terms ofG_iandF_i,

    ⟨|   = 2 H ϕ̇ X (3 G_5X + 2 X G_5X X - 60 X F_5 - 24 X^2 F_5X)
       + 2 ϕ̇/H[ - G_4ϕ + X( G_3X - 2 G_4ϕ X)]
       + 8 X [G_4X - G_5ϕ + X (2 G_4X X - G_5ϕ X + 8 F_4  + 4 X F_4X)] .

  *  The effective Planck mass,

    = 2 G_4 - 2 X[ 2 G_4X +  H ϕ̇ G_5X - G_5ϕ + 2 X ( 2 F_4 - 6 H ϕ̇ F_5)] ,

can be absorbed in the densities and pressures of matter and the scalar and be effectively hidden in the equations of motion.
But this is not the case for the Planck mass run rate,

    ≡ dln/ dln a = 
    2 Ḣϕ̇ X/H(-G_5X + 12 X F_5) 
        + 2 H ϕ̇ X[3G_5X + 2X (G_5X X - 30 F_5 - 12 X F_5X)]
        - 2 ϕ̈ X [3 G_5X + 2X(G_5X X - 30 F_5 - 12 X F_5X)] 
        - 2 ϕ̇ϕ̈/H[G_4X - G_5ϕ + X (2 G_4X X - G_5ϕ X + 8 F_4 + 4 X F_4X) ]
        + 2 ϕ̇/H[G_4ϕ + X ( G_5ϕϕ - 2 G_4ϕ X - 4 X F_4ϕ)]
        + 4X [ G_4X - G_5ϕ + 2X (G_4X X - G_5ϕ X + 4 F_4 + 2 X (F_4X + 3 F_5ϕ)) ] ,

which measures the variation in time of the Planck mass. This function is non-zero in models where in the action the scalar field couples directly to curvature, and it produces anisotropic stress in the gravitational potentials.


  *  The tensor speed excess,, signifies deviations of the speed of gravitational waves from the speed of light and has an effect on the scalar perturbations too.
As, it is present in non-minimally coupled theories and can be written as,

    = X[2 G_5X(2 H ϕ̇ - ϕ̈) +
    4(G_4X - G_5ϕ + 2 X (F_4 - 3 F_5 H ϕ̇))] .

  *is non-zero for beyond Horndeski theories.
It has been shown that even when ordinary mater is minimally coupled to the metric, in the theories beyond Horndeski, the higher order derivatives in the conformal transformation of the metric can lead to coupling between the sound speeds of dark energy and matter <cit.>.
This mixing between the scalar and matter can be measured by theparameter,

    = 8 X^2 (F_4 - 3 F_5 H ϕ̇) .


Here, the dots denote derivatives with respect to proper time.
As a shorthand notation in the rest of this work we may refer to these functions asα_i.

It is possible to derive the linear perturbation equations and an evolution equation for the perturbation of the scalar in term of these functions.
These equation, combined with the evolution equations for the matter species, provide the full system that one has to solve to obtain the linear matter and CMB power spectra and any other prediction within linear theory.
Given the constraining power of current datasets, any practical use of theα_iimplies the choice of a parametric form to describe their time evolution.
There are many different parametrisations that can be used.
Although they cannot fully represent the beyond Horndeski class of theories, we have to make a choice in order to evolve the equations of motion in time.
For example one could choose to keep theα_iconstant, split them in redshift bins, or evolve them as proportional to a power of the scale factora^nor the the fractional density of the DE/MG componentΩ_s^n, withn=const.
Here we pick the parametrisation proposed in <cit.>,

    α_i = Ω_sα̂_i ,

whereα̂_iare constants. We choose this time evolution because we are interested in DE/MG models that aim at explaining the late-time acceleration of the universe, so it is reasonable to assume that theα_iwill become more important as the DE density grows and the universe starts to accelerate.
If all theα̂_iturn out to be negligibly small, that would suggest that the accelerated expansion is mainly due to a cosmological constant.
Furthermore, this parametrisation provides more freedom than the constant or binnedα_i, but has a smaller parameter space to sample from than forα_i = Ω_s^n_iα̂_i. 
While the reasons we have invoked are practical, we note, however, that the study of specific beyond Horndeski models is in its infancy; specifically, there is, as yet, little understanding of the type of behaviour one might expect if constructing a model from the ground up (i.e. via an action). 
It would be important to explore the range of possible behaviours (as has been done for the case of the equation of state and the CLP parametrization in the case of dark energy and, to a lesser extent, as has been done for "normal" Horndeski models) – only then can we be assured that such a parametrisation is representative of a model space.



 §.§ Background


The Friedmann equations for a flat Friedmann-Robertson-Walker (FRW) metric read

    3H^2= +  ,
       2Ḣ=-3H^2-- ,

where the subscriptmstands for all matter species andsfor DE/MG.
Here we have absorbed the Planck mass into the definitions of the densities and the pressures of matter and the scalar, e.g.ρ̃≡ρ/M_*^2(see <cit.> for notation). As a consequence, the matter and DE densities are not conserved anymore

    ρ̇̃̇_ m + 3H( + ) = - H ,
       ρ̇̃̇_ s +3H(+)= H .

This can be interpreted as energy density being exchanged between matter and the scalar whenevervaries.
Thus in this framework we need to integrate the densities of all matter species appropriately according to the equations above.
With different definitions of DE/MG energy density and pressure we could have assumed the standard evolution for these quantities at the price of keeping track the value of.
Clearly the two approaches are equivalent and one uses whichever is more convenient. 
The expressions forandin terms of the HorndeskiG_iandF_ifunctions are given by,

    ≡ - G_2 +2X(G_2 - G_3ϕ) + 6ϕ̇H(X G_3X - G_4ϕ - 2X G_4ϕ X)
       + 12H^2X(G_4X + 2X G_4X X - G_5ϕ - X G_5ϕ X + 8X F_4 + 4 X^2 F_4X)
       + 4ϕ̇H^3 X(G_5X + X G_5X X - 24 X F_5 - 12 X^2 F_5X) ,

    ≡ G_2 -2X(G_3ϕ - 2G_4ϕϕ) + 4ϕ̇H(G_4ϕ -2X G_4ϕ X + X G_5ϕϕ - 4 X^2 F_4ϕ)
        -  ⟨H|ϕ̈/ϕ̇ 
            - 4 H^2 X^2 (G_5ϕ X - 12X F_5ϕ)
            + 2ϕ̇ H^3 X (G_5X - 12F_5 ) .


 §.§ Stability conditions


The price of introducing new degrees of freedom is that the resulting theory may be unstable on a given background.
Indeed, it is possible for specific choices of theα_ithat the perturbations grow exponentially.
In this section we show the conditions to avoid two types of instabilities, i.e. ghost and gradient.
Ghost instabilities happen when we choose the wrong sign of the kinetic term of a d.o.f., while gradient instabilities can occur when the sound speed of the fluid is imaginary.

After decoupling from auxiliary variables the quadratic action for scalar and tensor modes in beyond Horndeski with no matter reads <cit.>,

    S^(2)=∫ d^4xa^3[Q_ S(ζ̇^2-/a^2(∂_iζ)) + Q_ T(ḣ_ij ^2-c_ T^2/a^2(∂_kh_ij)^2)] ,

whereh_ijare the tensors modes (gravitational waves), andζthe scalar.
To ensure that the propagating degrees of freedom are not ghost-like, we require that their kinetic terms are positive, i.e.

    Q_ S≡2 D/(2-⟨)|^2>0 ,     D≡+3/2⟨^|2 ,
       Q_ T≡/8>0 .

To avoid gradient instabilities we need to require that the speed of sound of the scalar and tensor degrees of freedom are positive.
This means that,

    c_ s^2 =(2-⟨)|^2/2D{-c_ T^2 + 4/ad/dt[a(1+)/H(2-⟨)|]} - +/D H^2 (1+)^2>0 ,
       c_ T^2 = 1+  > 0 .

Note that the only condition affected by the beyond Horndeski terms is the gradient condition for the scalar sector, while the remaining are the same as for Horndeski models <cit.>.

There have been studies of the impact that stability conditions have on the constraints obtained for parametrised models <cit.> and  gives you the freedom to relax or disable the restrictions imposed by these conditions.
From a brief inspection of the results that we obtain with and without these restrictions we find that the combinations ofα_ithat get rejected by the stability conditions are also in regions of extremely low likelihood; if the tests are disabled they lead to exponential growth of the perturbations (and hence are an extremely bad fit to data).
While it is safe to ignore the restrictions,  we have decided to keep them in our analysis, as this speeds up sampling significantly.



 §.§ Recent implications from GW170817 and GRB 170817A

After the detection of an electromagnetic counterpart (GRB 170817A) to the gravitational wave signal (GW170817) from a binary neutron star merger <cit.>, it has been shown that the speed of gravitational waves (GW) has to be incredibly close to the speed of light and thus||<10^-15(we note that
Current constraints onusing cosmological data areO(1)).
This implies that, as pointed out in <cit.>, this result imposes strict constraints on general scalar-tensor and vector-tensor theories that are much tighter than current cosmological constraints and therefore it is safe to take 

    =0 ,   c_ T^2=1 .

It has been suggested that this  constraint can be avoided for some scalar-tensor theories if one takes into account the dynamics of the scalar field directly, when it is not coupled directly to matter <cit.>. 














Furthermore, in Ref. <cit.> the authors claim that any model with≠0predicts a copious decay of GW into DE scalar fluctuations. Given that that we do observe GW, suggests that we must rule out any beyond Horndeski model. The only way to avoid decay of GW and have≠0is to impose the propagation speed of scalar perturbations,, to be equal to the speed of light. While the authors mention that there could be power-law divergent terms even in this case, here our objective is to give a purely phenomenological idea of the behaviour of beyond Horndeski theories. Our main analysis is then performed lettingfree to vary and without any additional condition on(beside the stability condition mentioned in the previous section). However, for completeness we compare this case with the one where we fix=1and show the results in Section <ref>.

It is possible to get models with= 1by expressing one of theα_iin terms of the others. In this paper we choose to fix, which obeys this relation

    = -3/2⟨^|2 - 1/2(2-⟨)|^2{c_ T^2 - 4/ad/dt[a(1+)/H(2-⟨)|]} - +/H^2 (1+)^2 .

This ensures that the speed of sound of the scalar is unity for all the models we consider at all times.
As shown in <cit.>, the value ofis generally unconstrained by data and varying this parameter has negligible effect on the background.
In Section <ref>, we can see thatis not present in the equations of motion in the quasi-static approximation, Eqs. (<ref>) and (<ref>), so it also has no effect on the observables we are looking at in this limit.
In Section <ref> we present the results foras given by Eq. (<ref>) and compare those to the case where= 1.



 §.§ Perturbations

We assume that the universe is spatially flat and is well described by small perturbations to the Friedmann-Robertson-Walker (FRW) metric. 
Taking into account only scalar perturbations, the line element in Newtonian gauge (with the notation of <cit.>) reads,

    ds^2=-(1+2Ψ)dt^2+a^2(t)(1-2Φ)d𝐱^2 .

Following <cit.> we redefine the perturbations of the scalar field as

    v_X≡-δϕ/ϕ̇ .


Here we present the linear perturbation equations for the Fourier components of the action given in Eq. (<ref>) in Newtonian gauge and physical time. 
In Appendix <ref> we present these in synchronous gauge, as implemented in .
These equations have been derived previously in <cit.>.

The modified Einstein (00) equations has the form,

    3(2 - ⟨|) H Φ̇ + (6 -  - 6 ⟨|) H^2 Ψ + (1 + )2 k^2/a^2Φ
        - ( + 3⟨|)H^2v̇_̇Ẋ - H[(⟨+| 2 )k^2/a^2 - 3 ⟨ ̇|H + 3( 2Ḣ +  + ) ]v_X  = - δ_m .

the Einstein (0i),

    2 Φ̇ + (2 -  ⟨|) H Ψ - ⟨H|v̇_̇Ẋ - (2 Ḣ +  + )v_X  = - ( + )v_ m ,

the traceless part of the Einstein (ij) equation,

    (1 + ) Ψ - (1 +  ) Φ + v̇_̇Ẋ - (  - )H v_X = - 3/2( + ) σ_m ,

and the trace part,

    2 Φ̈   -  ⟨H|v̈_̈Ẍ + 2(3 + ) H Φ̇ + (2 - ⟨|)HΨ̇
       + [(2 - ⟨|)(3 + )H^2 -  (⟨H|)  - (2Ḣ +  + )]Ψ
       - [⟨|(3 + ) H^2 + (⟨H|)  +  (2Ḣ +  + )]v̇_̇Ẋ
       - [2 Ḧ + 2(3 + ) Ḣ H +  ṗ̃̇_ m +  H ] v_X  = - k^2/a^2σ_m ( + ) + δ ,

and the equation of motion for the scalar velocity potentialv_X,

    3 ⟨ ̈|Φ H    +  H^2 v̈_̈Ẍ + 
        3[⟨|(3 + ) H^2 + (⟨H|)  - (2Ḣ +  + )] Φ̇ + ( + 3 ⟨|)H^2Ψ̇
        -  2k^2/a^2Φ̇ - 2 [ + ((1 + )+  - ) H ]k^2/a^2Φ  - (⟨+| 2 )H k^2/a^2Ψ
       + [(2 + 9 ⟨|) Ḣ + ( + 3 ⟨̇|̇) H + (3 + )(3⟨+|)H^2 -3( 2Ḣ +  +)]HΨ
       -3 {3Ḣ[(⟨H|)  - (2Ḣ +  + )] + 3 ⟨|[Ḧ +  (3 + )ḢH]H} v_X
       + {[(⟨+|2)H]  + (⟨+| 2 )(1 + )H^2 + 2 ( - )H^2 - 2 (Ḣ +  + )}k^2/a^2v_X
       +[2 Ḣ +  H + (3 + )] Hv̇_̇Ẋ = 0 .

These equations, combined with the evolution equations of the matter perturbations,δ_m,v_m,δ, andσ_m, form the complete set of equations that one has to solve to follow the linear dynamics of beyond Horndeski theories.
We have implemented these equations into <cit.> and used the solutions to study the phenomenology of this class of models and constrain its parameters with current CMB, large-scale structure data in the following two sections.



§ PHENOMENOLOGY OF BEYOND HORNDESKI THEORIES

The purpose of this work is to provide constraints on parametrized beyond Horndeski models, i.e. using.
In this section we briefly discuss the observables we use to calculate our parameter confidence regions, namely the matter power spectrum and the CMB temperature, polarisation and lensing power spectra.
We present the results obtained with  for the evolution of these quantities for different values ofand derive the analytic expressions for these in the quasi-static approximation and observe what effectshas on them.



 §.§ The quasi-static approximation

The key assumption of the quasi-static approximation (QSA) is,

    |Ẋ|≲ H|X| ,

whereXstands for any metric or scalar perturbation.
This assumption states that the time evolution of all parameters and perturbations of the metric and the scalar is comparable to the evolution of the Hubble parameter.
However, this can be only true on sub-horizon scales, where space-derivatives become important

    k^2/a^2|X|≫ H^2|X| .


We can apply this limit to the Hamiltonian constraint equation,  Eq. (<ref>), and the equation of motion for the scalar perturbation, Eq. (<ref>),

    2 k^2/a^2(1 + )Φ - k^2/a^2(⟨+| 2 )H v_X = - δ_ m ,

    2   k^2/a^2Φ̇ +  (⟨+| 2 )k^2/a^2 H Ψ + 2[ + ( (1 + ) +  - )H]k^2/a^2Φ
       + {[(⟨+|2)H]  + [(⟨+| 2 )(1 + ) + 2 ( - )]H^2 - (2Ḣ +  + )}k^2/a^2v_X = 0 .

In this approximation, we neglect any contribution from the matter anisotropic stress,σ_m, and pressure,.
The reason is that the QSA can give predictions only at late-times, where the dominant matter component is DM, which is modelled as a pressureless perfect fluid.

Combining these with the rest of the constraint equations Eqs. (<ref> - <ref>) and their time derivatives, we get the constraint equations for the perturbationsΦandΨin this limit,

    k^2/a^2Φ = -3/2H^2Ω_ mμ_Φδ_ m + λ_Φ Hk^2/a^2v_ m ,
       k^2/a^2Ψ = -3/2H^2Ω_ mμ_Ψδ_ m + λ_Ψ Hk^2/a^2v_ m ,

where

    μ_Φ    = 1/1+(1+⟨+|2/B) ,
    λ_Φ    = -1/2(⟨+|2)/ H^2 ,
    μ_Ψ    = 1/(1+^2A)(1+)^2[c_ T^2 + 2(B+/H)B/+(1+)ad/dt(2B/aH)] ,
    λ_Ψ    =Ȧ+2A/1+^2A ,

    A≡2λ_Φ/(⟨+|2)H ,    
    B≡1/2(⟨+|⟨|+ 2 - 2 - 2 - 2/H) .

These can be combined with the Euler and continuity equations for DM, i.e. 

    δ̇_ m=k^2/a^2v_ m ,    v̇_ m=-Ψ .

Note that the notation we have used here is different from the commonly used effective gravitational constant,μ≡G_eff/Gand gravitational slip parameter,γ≡Φ/Ψ.
In beyond Horndeski we have additional terms involving matter velocity that can not be neglected on sub-horizon scales, i.e. the ones proportional toλ_Φandλ_Ψ.
To avoid confusion we chose to useμ_Φandμ_Ψinstead.
In Horndeski (=0), we would be able to relate the variables used here to the usualμandγas

    μ_ϕ≡μ     μ_Ψ≡μ/γ .

We know that in theΛCDM model these are both equal to one.
So the strength of gravity in other models can be expressed as a the gravitational constant inΛCDM plus a small correction <cit.>, 

    μ = 1 + δμ = -2k^2Φ/3H^2Ω_ mδ_ m .

However, here due to the presence of the beyond Horndeski parameter,, this is not as straightforward. 
The effect of this additional term is non-negligible and we have taken it into account in all the calculations here.
To show the significance of these terms we look at the simplified case where⟨=|==0and= Ω_DE, i.e. a model with a standard kinetic termplus beyond Horndeski. The functionsλ_Φandμ_Φsimplify to

    λ_Φ= -3Ω_ m/2-3(-1)Ω_ m ,    μ_Φ= 1 - λ_Φ ,

and forλ_Ψandμ_Ψ, we have

    λ_Ψ=-9Ω_ m(2-4Ω_ m-3Ω_ m^2)/(2+3Ω_ m)[2 + 3(1-)Ω_ m] ,    μ_Ψ=1 - λ_Ψ .μ_Φandμ_Ψare plotted on Figure <ref> as functions of the scale factor,a. 
In order to isolate its contribution, allα_iexcept forhave been set to zero.
Given that we chose the evolution ofto be proportional toΩ_DE, we see that the the deviations fromμ_Φ= μ_Ψ= 1are very small before the onset of DE. 
On the left panel we have the evolution ofμ_Ψ= 1 -  λ_Ψ, which begins to decrease below one asbecomes non negligible at late times.
On the right isμ_Φ= 1 -  λ_Φ. At late times λ_Φbecomes negative and enhancesμ_Φ.


In this case we can writeδμ_Φ≡-λ_Φand so the modifications to GR coming from the two terms on the RHS of Eq. (<ref>) are3/2H^2Ω_mλ_Φδ_mandλ_ΦHk^2/a^2v_m, and equivalently for Eq. (<ref>).
The ratio of these two contributions is plotted on Figure <ref> for four different values of.
We can see that the these are of the same order and hence have roughly equal contributions to the equations. 
This means that thev_mterm cannot be neglected and must be taken into account in the definition of the QS approximation equations.
It is possible to redefine the effective gravitational constant and hence the constraint equations above using the common notation and taking into account the additional terms we have here.
As usual, we define the growth factor as the logarithmic derivative of the matter overdensity

    f = d lnδ_ m/d lna .

Using,

    k^2/a^2v_ m = δ̇_ m = f Hδ_ m ,

we can re-write the constraint equation forΦin its usual form

    k^2/a^2Φ = -3/2H^2Ω_ mμδ_ m ,

where

    μ = μ_Φ - 2/3Ω_ mλ_Φ f .

We show in Section <ref> that the approximate solution for the growth factor in the quasi-static limit is very close to the exact solution, which makes this a good expression for the effective gravitational constant.



 §.§ Matter power spectrum


Combining Eqs. (<ref>) and (<ref>), we can derive the evolution equation for the matter overdensity,δ_m, in the QS limit

    δ̈_ m + (2+λ_Ψ)Hδ̇_ m-3/2H^2Ω_ mμ_Ψδ_ m=0 .

As we have seen above the presence ofdecreases the effective gravitational constantμ_Ψbyλ_Ψ, which in our parametrisation becomes important only at late times (see Figure <ref>).
Then, increasing the value of, we expect a suppression of the matter power spectrum at small scales.
This has two contributions acting in the same direction: (i) a decreased value of the effective gravitational constantμ_Ψ, and (ii) an enhanced friction term proportional to.
In the matter-dominated era, wherea∝t^2/3,H=2/(3t)andλ_Ψ≪1, an approximate solution to Eq. (<ref>) up to linear order inis given by

    δ_ m∝ t^2/3(1-λ_Ψ) ,

where in the derivation we assumed constant coefficients. 
This confirms that, during this period, whileλ_Ψ>0, the growth of overdensities and hence the matter power spectrum have to be suppressed.

In Figure <ref> we show the matter power spectrum as a function ofkat redshiftz=0as outputted from . 
We have set the other parameters,=1,===0, in order to observe the effect ofalone.
We set the background cosmological parameters to their best fit values form Planck 2015 <cit.> and assumeΛCDM for the background expansion history.
In this case only positive values ofare allowed, as negativewould cause a gradient instability.
AssumingΛCDM background expansion we fixw_s=p_s/ρ_s=-1, sop_s+ρ_s=0andρ̇_s = 0, and given that we chose=Ω_s, the gradient instability, Eq. (<ref>) simplifies to,

    [(ρ_m+p_m)(1-) + 2H^2]>0 ,

which is only satisfied when

    >0     and     < 1 + 2H^2/ρ_m+p_m .

For the time evolution of theα_iwe assumed parametrisation in Eq. (<ref>).
On the left panel of Figure <ref> we show the matter power spectrum generated by  for five different values ofand on the right the relative difference between those and our fiducialΛCDM.
This is consistent with the results of <cit.>.
As expected and discussed in detail in this section,has a damping effect on the matter power spectrum and the deviation fromΛCDM on small scales indeed seems proportional to.



 §.§ The growth rate of structure
 
Using the definition of the growth rate of structure, Eq. (<ref>), and the evolution equation for the matter overdensity, Eq. (<ref>), we have,

    f^' + f^2 + (2 + λ_Ψ + Ḣ/H^2)f - 3/2Ω_ mμ_Ψ = 0 ,

where primes denote a derivative with respect tolna.
We can find an approximate solution to this equation in the matter-dominated era, making the following assumptions.
The growth rate of structure can be expressed as the growth rate of matter in matter domination plus a small correction that represents the DE/MG component

    f = f_ m + δ f .

Similarly, as argued above, we can express the effective gravitational constantμ_Ψ, as its standard value, i.e.μ_Ψ=1, plus a small correction,δμ_Ψ≡-λ_Ψ.
Expanding the growth equation Eq. (<ref>) to first order in the perturbation parameters, we get,
   0^ th order    f^'_ m + f_ m^2 + (2 + Ḣ/H^2)f_ m - 3/2Ω_ m = 0 ,
   1^ st orderδ f^' + 2f_ mδ f + -δμ_Ψ f_ m + (2 + Ḣ/H^2)δ f - 3/2Ω_ mδμ_Ψ = 0 .      
We know that inΛCDM matter-dominated universeΩ_m∼1,H=2/(3t)andf_m^'=0, sof_m=1.
Hence we can find an approximate solution for the deviation of the growth rate fromΛCDM,

    δ f ≃   5/2 e^-5/2ln a∫e^5/2ln aδμ_Ψ dln a
       ∝δμ_Ψ ,

where to obtain the second line we have assumed thatδμ_Ψ∼const., which is reasonable in matter domination where in most realistic models theα_iare negligible.

We plot this result in Figure <ref> and compare it to the one obtained with .
We can see that in the matter-dominated era, the two solutions agree very well and the approximation fails at later times, where the assumptions we have made do not hold.



 §.§ CMB power spectra



















The fluctuations in the CMB temperature and polarisation can be used to constrain very tightly the standard cosmological parameters in theΛCDM model.
Large scale structure gravitationally lenses the temperature and polarisation anisotropies, which makes the CMB also a good probe of the late universe where the effects that we are interested in become important. 
Therefore, beyond Horndeski theories are expected to modify the CMB spectra. In particular we expect to see significant deviations fromΛCDM in the lensing potential and the temperature power spectrum of the CMB.
On the left panel of Figure <ref> we show the relative difference between the CMB lensing potential,C_l^ϕϕ, of a set of MG models and our fiducialΛCDM.
As before we are looking at models where all===0and=1and consider five different values of.
We can see thatC_l^ϕϕdecreases at a nearly linear rate as we increase, which was also shown in <cit.>.
Again, in order to understand this behaviour, we turn to the equations describing this effect in the quasi-static approximation.
The strength of the lensing potential is related to the Weyl potential, which can be obtained by summing the constraint Eqs. (<ref>) and (<ref>),

    k^2/a^2(Φ+Ψ)=-3/2H^2Ω_ m(μ_Φ+μ_Ψ)δ_ m +  H (λ_Φ + λ_Ψ)δ̇_ m .

Similarly to the case of the power spectrum, Section <ref>, we look at the simplest case (where⟨=|==0). Here we haveμ_Φandμ_Ψgiven by Eqs. (<ref>) and (<ref>) respectively.
As in the case with the effective gravitational constant, the weak lensing parameter can not be defined as a function ofμ_Ψandμ_Φas usual.
Here, instead we have

    Σ≡ -2k^2(Φ+Ψ)/3a^2H^2Ω_ mδ_ m   =μ_Φ + μ_Ψ - 2(λ_Φ + λ_Ψ)/3Ω_ mf ,

where again we have usedδ̇ = fHδ.
For the case of===0and=1the above expression reduces to

    Σ = 2 - (λ_Φ + λ_Ψ)(1+2f/3Ω_ m) .

This expression is plotted on the right panel of Figure <ref>.
InΛCDM the lensing potential isΣ= 2, so what we are showing is the modification to this quantity due to the beyond Horndeski parameter,.
It is possible to notice that this is negative during most of the matter dominated era and decreases asincreases, which confirms the damping effect ofon the lensing potential, see left panel of Figure <ref>.

The effects of beyond Horndeski can be seen also on the top panel of Figure <ref>, where we plot the CMB temperature power spectrum as obtained by .
This result is again consistent with <cit.>.
We can see that increasingenhances the CMB temperature power spectrum at lowl.
This effect is well known for DE/MG models, and it is due to the integrated Sachs-Wolfe effect, which affects the power spectrum at lowl, but has no significant effect at largel. The contribution from the combination of the late and early integrated Sachs-Wolfe effect to the temperature anisotropies is plotted on the bottom of Figure <ref>.

















































§ RESULTS

In this section we present the constraints on our specific choice of beyond Horndeski parameters with data. We note that different time dependences may lead to large changes in the uncertainties and that a more systematic analysis including a broader range (or variety) of priors would be desirable. Nevertheless, this can be seen as a first step in assessing constraints on par with what has been done in the case of Horndeski theories.


 §.§ Method

We use the Boltzmann code <cit.> to solve the equations shown in Appendix <ref>.
This code extends the Cosmic Linear Anisotropy Solver Software (CLASS) <cit.> by including the Horndeski class of theories, and it has been tested intensively against other Einstein-Boltzmann solvers in <cit.>.
On top of the public version of , we have added the extra beyond Horndeski terms to the perturbation equations, the definitions of the DE density, pressure and.
We then interfaced it with the modular cosmological parameter estimation code CosmoSIS <cit.>, which runs Monte Carlo Markov Chains (MCMC).
We used the Metropolis-Hastings sampler to obtain our chains and consider those converged when the Gelman and Rubin parameterR-1<0.01<cit.>.

The full set of beyond Horndeski parameters that describes the evolution of the perturbations is{,,,,}.
Assuming the constraints derived from the GW detection that we discussed in Section <ref>, we set=0for all runs and run two sets of MCMC chains, one fixing=1and another withgiven by Eq. (<ref>) to ensure=1for all possible sets of values of{,,}. 
We fix the initial effective Planck mass,=1.
We also assume two massless and one massive neutrino withm_ν=0.06eV and vary the background cosmological parameters and the relevant set ofα̂_i.
Given that we use a parametrisation for the evolution of theα_i, Eq. (<ref>), here we present the constraints for the coefficients that fix this proportionality,α̂_i.



 §.§ Datasets


To constrain the parameters of the theory we include CMB data, measurements of baryon acoustic oscillations (BAO) and redshift space distortions (RSD) from large-scale structure surveys:

  *  CMB: We use the temperature and polarisation power spectra and the lensing potential of the Cosmic Microwave Background (CMB) from Planck 2015 <cit.>. 
    The Planck 2015 likelihood is discussed in <cit.>. We use the high l TT likelihood with l=30-2508, along with the joint TT, EE, BB and TE likelihood in the range l=2-29 and the lensing likelihood, using both temperature and polarisation lensing reconstruction in the multipole range l=40-400.

    
  *  BAO: 
    We use BAO measurements from 6dFGS <cit.> and the Sloan Digital Sky Survey (SDSS) Data Release 7 Main Galaxy Sample (MGS) <cit.>, which probes the supplied expansion history through the redshift-distance and redshift-Hubble relations combined through the relation, D_ V = (D_ A^2(1+z)^2z/H)^1/3. 
Here D_ V is the angle-averaged distance, D_ A - the angular diameter distance and z is the redshift.
We also used measurements from the Baryon Oscillation Spectroscopic Survey (BOSS) Data Release 12 <cit.>, which are lower signal-to-noise and constrain both the angular diameter distance, D_ A and the Hubble parameter, H.
There are also BAO measurements from WiggleZ survey <cit.>, but these partially overlap with the BOSS CMASS volume. We have chosen to use the BOSS measurements here.
These measurements do not constrain the matter power spectrum itself, but just the positions of the BAO peaks.
They are presented in Table <ref> and plotted on the left panel of Figure <ref>.

    
  *  RSD: Finally, we use Redshift Space Distortions (RSD) data. This effect probes the growth of structure and therefore is sensitive to late-time effects and can be used to probe MG.
We use measurements derived form from 6dFGS <cit.>, SDSS MGS <cit.> and Luminous Red Galaxy (LRG) <cit.> samples, WiggleZ <cit.>, and VIPERS <cit.>.
We also use the RSD measurements from BOSS DR12 <cit.> with the full covariance between the 3 f(z)σ_8(z) measurements at different redshifts and the BAO measurements of H(z) and D_ A(z).
All RSD measurements are quoted on Table <ref> and on the right panel of Figure <ref>.

In <cit.> the authors also included data for the full matter power spectrum,P(k), from a number of large scale structure surveys.
However, this dataset does not improve the constraints significantly and requires to choose a galaxy bias factor. So we have decided not to include it here.


 §.§ Constraints on the beyond Horndeski parameters

Here we present the constraints on the beyond Horndeski parameters derived using the method and datasets described above.
As mentioned before, the value ofhas little to no effect on sub-horizon scales and in particular on the observables we are looking at.
We then chose either to fix=1or to get it from Eq. (<ref>) and present the constraints on the other parameters.

The grid of parameters varied and datasets used is shown on Table <ref>, with the 95% confidence regions forand marginalised means.
Our baseline model is the standardΛCDM, where we varied the Hubble constantH_0, the baryon and dark matter densities today,Ω_bandΩ_cdm, the curvature fluctuation amplitude,A_s, the scalar spectral index,n_sand the optical depth at reionization,τ.
On top of that, in the other runs we varied: the beyond Horndeski parameterkeeping the otherα̂_ifixed (ΛCDM+);and the density of massive neutrinos (ΛCDM+Ω_ν+);and the Horndeskiα_i(ΛCDM+α_BMH);and we replaced the cosmological constantΛwith a time dependent equation of state parametrized asw(a)= w_0 + w_a(1 - a)(wCDM+α_H); and finally a time dependent equation of state plus the Horndeskiα̂_iand(wCDM+α_BMH).

The confidence regions ofseem to exclude0in most of these cases and are mainly positive.
The lower limit at0is exact when considering a cosmological constant and==0, and can be inferred by inspecting Eq. (<ref>).
For the remaining cases the contours are smooth around0since the other MG parameters modify the stability conditions, but Eq. (<ref>) still plays a central role in excluding negative values of. 
Both the upper and lower limits seem to be affected largely by the combination of datasets used and the freedom in some of the other parameters.
We see that varyingΩ_νhas little effect on the constraints ofand the largest difference seems to be due to the other MG parametersand.
It is possible to notice from Table <ref>, and Figures <ref> and <ref>, that the contours of the MG parameters are greatly improved by adding the RSD dataset.
This can be explained by the fact that RSD is the only dataset that probes the growth rate of structure,fσ_8, at high significance.
And given that we chose for theα_ito evolve as the dark energy density parameter,Ω_DE, we would expect it to have an effect at late times.


In Table <ref> we show the absolute log likelihood values for the models considered, measured at their best fit values.
It is possible to notice that the total likelihood generally decreases as we increase the number of free parameters.
A better fit to data is expected when enlarging the parameter space, but in this case the improvement does not appear enough to justify the price of adding new parameters.
The best fit models are also plotted, together with the measurements we used, on  the left panel of Figure <ref> for the BAO datasets and on the right for the RSD dataset.
We show how these measurements are fitted by the models we have considered here (green lines) andΛCDM (black dashed lines). 
The beyond Horndeski models, and particularly when combined with the DE equation of state clearly provide more flexibility to fit these measurements than the best likelihood model ofΛCDM.
However, given that most of these models fit within the error bars of these measurements, we cannot conclude that these differences point towards either of the models being a better fit.


On Figure <ref> we have plotted the confidence regions for the relevant parameters of theΛCDM+Ω_ν+run.
The lighter and darker contours represent the 95% and 68% confidence regions respectively.
In this case the only MG parameter we vary is, while the remainingα̂_iare fixed to= 1,===0. As mentioned the gradient stability condition, Eq. (<ref>), requires>0, which can be clearly seen in the contours. Interestingly,does not seem to have any degeneracy with the density of neutrinosΩ_ν, which was the main motivation to study this extension to the baselineΛCDM model.



Figure <ref> shows the contours for thewCDM+α_BMHmodel. In particular, we focus on the posterior distributions of the (beyond) Horndeski parameters,{,,}, and the DE EoS parameters,w_0andw_a. As mentioned, by comparing the two dataset combinations, CMB+BAO and CMB+BAO+RSD, it is possible to notice the effect of RSD in tightening the contours.
The most significant differences seem to be forand.
Unlike in the previous case, here there is no hard bound on0on the value ofcoming from the stability conditions. This allows it to go negative, since the combination with the other parameters may preserve the stability of the theory.
In the case where we had== 0, the data seem to peak the distribution ofaround zero and no large values were allowed. On the contrary here, the degeneracy between,and, spreads out the distribution ofin the positive direction relaxing significantly the contours. We see thatis not degenerate with thew_DE.


On Figure <ref> we show the different constraints on the parameters for the two cases discussed in Section <ref>. 
The combination of datasets used for both contours here is CMB+BAO+ RSD and the model considered isΛCDM+α_BMH.
The red contour represents the case where we set=1and let the speed of sound of the scalar degree of freedom vary with the otherα̂_i.
The blue is for the case where the speed of sound of the scalar=1andis determined by the relation in Eq. (<ref>).
We can see that there is no significant difference in the contours for those two cases, which confirms that the choice forhas little effect on current constraints, and that the=1branch of the results of <cit.> may not have an effect on the estimated values of the rest of theα_i.



§ DISCUSSION


In this paper, we studied the phenomenology of the beyond Horndeski class of theories.
This class of theories includes the most popular DE and MG models that aim at explaining the accelerated expansion of the Universe and has been shown to have self-accelerating solutions.
The minimal set of functions that can describe the the evolution of linear perturbations in this model is{,⟨,|,,}.
To find the solutions of these equations a parametric form for these functions has to be assumed.
Here we chose the time evolution to be given by the fractional density of dark energy, so thatα_i = α̂_i Ω_DE.
Taking into account the implications from the observation of the binary neutron star merger (GW170817) and its electromagnetic counterpart (GRB170817A), we have set the tensor speed excess,= 0.

Constraints on,,andfrom cosmological data are𝒪(1).
In this work our aim was to constrain the parameter, which is peculiar of the theories beyond Horndeski, using a combination of CMB, BAO and RSD datasets.

Studying the phenomenology of this model and in particular the case where=1,===0, we found thatsuppresses the matter power spectrum and the CMB lensing potential, and enhances the temperature power spectrum of the CMB at lowl.
These findings are consistent with other results in the literature.
The quasi-static approximation was used to derive analytic solutions to evolution equations of the matter overdensity,δ_mthe growth rate parameter,fand the weak lensing parameter,Σ, which agreed with the exact solutions obtained using .
We discuss the commonly used notation for the quasi-static approximation and suggest how that may be amended to account for the extra terms that come from the presence of the beyond Horndeski parameter,.

We perform parameter estimation using CMB temperature, polarization and lensing data from the Planck 2015 dataset and measurements of the BAO and RSD from a number of different large scale structure surveys. 
From the MG parameters we set=1, as varying this parameter has been shown to have no effect on the other parameters, as mentioned above we also set=0, and varied,and.
We performed different runs where we variedand the six standard cosmological parameters,h_0,Ω_b,Ω_cdm,A_s,n_sandτ, combined with eitherΩ_ν,and,w_0andw_a, or a combination of all.
We present the confidence contours for two of these runs and show that the beyond Horndeski parameter,, is degenerate with the other MG parameters,and, but not with either theΛCDM parameters,Ω_ν,w_0orw_a.
The constraints onthat we get from these runs are𝒪(1).
We also show that setting the speed of sound of the scalar=1does not have a substantial effect on the results.
Comparing the likelihood of the best fit models for the different runs we have done, we find that the MG models seem to fit the data slightly better thanΛCDM.
However, the difference in the likelihood is too small to comment on the validity of either model over the others and most probably is simply due to the fact that there are more free parameters in the MG models than there are inΛCDM.




§ ACKNOWLEDGEMENTS

We are very grateful to Joe Zuntz for helping us in the implementation of  in . We thank Erminia Calabrese for useful comments and discussions. We are also grateful to Marco Crisostomi and Guillem Domenech for useful comments on this manuscript. DT, EB and PGF are supported by ERC H2020 693024 GravityLS project, the Beecroft Trust and STFC.

vancouver

§ IMPLEMENTING THE EQUATIONS INTO 

In this section we present the equations we implemented into . In CLASS (and ) the units are chosen so that the Friedmann equations are written as,

    H^2 =  +  ,
       H^'/a = -3/2(H^2 +  + ) ,

where the primes denote a derivative with respect to conformal time, however the Hubble parameterHis the physical one,andare the density and pressure of matter in the universe (excluding dark energy), andandare those of the scalar field.

    ≡    - 1/3 K + 2/3 X (K_X - G_3ϕ) -  2 H^3 ϕ^' X/3 a[7G_5X + 4 X (G_5X X - 33F_5 - 12 X F_5X)]
       +H^2[1 - (1-⟨|) - 4X(G_4X - G_5ϕ) - 4X^2(2G_4XX - G_5ϕ X + 8 F_4 + 4X F_4X)] ,

    =   1/3 K - 2/3 X (G_3ϕ - 2G_4ϕϕ) + 4 H ϕ^'/3 a(G_4ϕ - 2X G_4ϕ X + X G_5ϕϕ - 4X^2 F_4ϕ) -  (ϕ^'' + a H ϕ^')/3 ϕ^' a  H ⟨
    |    - 4/3 H^2 X^2 (G_5ϕ X - 12 X F_5ϕ) - (H^2 + 2 H^'/3 a)(1-) + 2 H^3 ϕ^' X/3 a(G_5X - 12 X F_5) .

This form of the densities and pressures is analogous to the way they were defined in <cit.>, so that we get the usual conservation equations,

    ^' = - 3  H a ( + ) ,
       ^' = - 3  H a ( + )  .

In CLASS it is possible to solve the linear perturbation equations in both Newtonian and synchronous gauge.
However, the current version of  only has the synchronous gauge, so here we present the equations of Section <ref> in this gauge as they are implemented into the code.

The line element to first order in the perturbations in synchronous gauge is <cit.>,

    ds^2=a^2[-dτ^2+(δ_ij+h̃_ij)dx^idx^j] ,

with

    h̃_ij(τ,k )=k̂_ik̂_jh+6(k̂_ik̂_j-1/3δ_ij)η+h_ij ,

in Fourier space, wherehandηare the scalar perturbations andh_ijis the tensor perturbation. 

The scalar field perturbationV_Xhere is defined as 

    V_ X≡ aδϕ/ϕ^'= a v_ X .

In this gauge the Einstein equations take the form:

Einstein (0,0)

    h^' =    4 k^2/H a(1 + /2 - ⟨|) η + 6 δ_ mc a/H (2 - ⟨|) - 2 H a ( 3 ⟨+|/2 - ⟨|) V_X^'
        - 2 [3 a H^' + ( + 3 ⟨|/2 - ⟨|)a^2 H^2 + 9 a^2/( + /2 - ⟨|) + (⟨+| 2 /2 - ⟨|)k^2]V_X .


Einstein (0,i)

    η^' = 3 a^2 θ_ m/2 k^2  + 1/2⟨H| a V_X^' + [ H^' a + 1/2⟨H|^2 a^2 + 3 a^2/2 ( + )] V_X .

where 

    θ_ m≡ - ( + ) k^2/a v_m .


Einstein (i,j) traceless

    ξ^' = (1 + ) η - a H (2 +  ) ξ + a H (  -   - ) V_X  -  V_X^'  - 9 a^2 σ_ m/2 k^2  .

To simplify this equation a new perturbationξ(τ,k⃗)has been introduced that relateshandηthrough

    ξ=h^'+6η^'/2k^2 .


Einstein (i,j) trace

    D h^'' =     2 (λ_1 - λ_9) k^2 η - 6 ⟨|k^2 η^'/H a + 2 a H λ_3 h^' -  9 δ a^2/
         + (3 a^2 H^2 λ_4 -2  k^2) V_X^' + 2 a H k^2 (λ_10 + λ_5 + 3 H^2 a^2/k^2λ_6)  V_X ,


And the equation of the evolution of the scalar perturbations is given by

    D(2 -  ⟨|) V_X^'' + [⟨|k^2/H a(2-⟨|) + 8 H a λ_7] V_X^' + [2 k^2 λ_11 + 2  k^2 - 8 H^2 a^2 λ_8] V_X 
       = 2 k^2/H a( + 3 (+)/H^2 M_*^2(1+)) η + 2  k^2/H^2 a^2(2 - ⟨|) η^' + 3 a /2 H M_*^2[ 2 λ_2 δ - 3 ⟨|(2-⟨|) δ]  .

where 

    D =  + 3/2⟨^|2 ,
       λ_1 = (1+) - 3⟨(|-) ,
       λ_2 = -3( + )/H^2 - (2 - ⟨)|H^'/aH^2 + ⟨^|'/aH ,
       λ_3 = -1/2(2+)D-3/4⟨λ|_2 ,
       λ_4 = λ_2 - 2⟨^|' - ⟨|^'/aH ,
       λ_5 = 3/2⟨^|2(1+)+(D+3⟨)|( - ) + 3/2⟨λ|_2 ,
       λ_6 = (1-3⟨H|^'/ a H^2)λ_2/2-D H^'/aH^2(2+ + H^''/a H H^') - 2⟨^|'-⟨|^'/2aH - 3^'/2aH^3 ,
       λ_7 = D/8(2-⟨|)(4++2H^'/aH^2+D^'/aHD) + D/8λ_2 ,
       λ_8 = -λ_2/8(D-3λ_2 + 3⟨^|'/aH) + 1/8(2-⟨)|[(3λ_2-D)H^'/aH^2-9⟨|^'/2aH^3]
           - D/8(2-⟨|)[4++2H^'/aH^2+D^'/aHD] ,
       λ_9 =  3[⟨|( 1 + ) +  ⟨|^'/a H] ,
       λ_10 = λ_9 + 3 ⟨|H^'/aH^2 -  ,
       λ_11 =[1/2(2-⟨)|(⟨+|2H^'/a H^2) + 3 (+)/H^2 M_*^2 (1+)] ,
       c_ sN^2 = 1/2(2-⟨)|[⟨+| 2 +2(1+) - (2-⟨)| + 2^'/a H - 2(1+) H^'/a H^2]
           + (1+)⟨^|'/a H - 3 (+)/H^2 M_*^2 (1+)^2 ,

wherec_sN^2is the numerator of the speed of sound of the scalar,

    c_s^2=c_ sN^2/D .

To include the effects of beyond Horndeski theories, we added the terms involving,λ_9,λ_10andλ_11to the speed of sound of the scalar and perturbation equations and theF_iand their derivatives to the definitions of the scalar density and pressure and theα_iin .

For now we have only included two parametrisations for this model in , one where the evolution ofα_iis given by the evolution of the DE density parameterΩ_DEand and one where they are proportional to the scale factor,

    α_i=(1-Ω_ m)α̂_i ,
       α_i= a α̂_i .

For both of these you need to set or vary the initial set{,,,,}.
For the first one we have also added the option to have the value ofcomputed using Eq. (<ref>).
In this case the initial set of parameters you need to specify is{,,,}.