Simple numerical implementation of general dark energy models

We present a formalism for the numerical implementation of general theories of dark energy, combining the computational simplicity of the equation of state for perturbations approach with the generality of the effective field theory approach. An effective fluid description is employed, based on a general action describing single-scalar field models. The formalism is developed from first principles, and constructed keeping the goal of a simple implementation into CAMB in mind. Benefits of this approach include its straightforward implementation, the generality of the underlying theory, the fact that the evolved variables are physical quantities, and that model-independent phenomenological descriptions may be straightforwardly investigated. We hope this formulation will provide a powerful tool for the comparison of theoretical models of dark energy with observational data.


Introduction
Recent cosmological observations have concreted the notion of the dark sector: some unknown substance or new gravitational physics is currently dominating the gravitational dynamics of the universe [1,2]. There are no shortage of theories which have been constructed in an attempt to describe these observations (see [3,4] for reviews).
Given the recent deluge of data [5][6][7][8][9][10][11] and upcoming experiments [12][13][14], it is of paramount importance that the understanding of dark sector theories is optimised for confrontation with observations. As such, it is becoming clear that some framework should be constructed which is capable of confronting entire classes of theories with data, and an industry has begun to construct different flavours of such a formalism [15][16][17][18][19][20][21][22][23][24][25][26][27]. Whatever the formalism is, it is important that it allows observations to be transcribed into well defined and meaningful statements about the allowed properties of the dark sector.
The aim of this paper is to present a leap forwards in the development of a model independent framework which can cover as much of the "known theory space" as possible. We do this by bringing together the effective field theory (EFT) approach to linear cosmological perturbations with the formalism for constructing the equations of state for perturbations. The former approach [28][29][30][31] is well suited for incorporating complicated theories, whereas the latter [32][33][34][35] is optimized for comparing theories to data.
The "holy grail" would be to provide a parameterization that covered all of Horndeski's theory [36][37][38], since that is the most general single-scalar field theory with second order field equations. In this paper we provide the penultimate step in such an endeavour: our calculations and results are applicable to almost all of Horndeski's theory, although the formalism can be straightforwardly extended. Our emphasis is torn between theoretical generality and usability -we provide a simple scheme for modifying numerical codes, such as CAMB [39], which will allow observational spectra to be extracted from very complicated models with a minimum of effort. This is entirely due to the way that the equations of state for perturbations work: they modify the fluid equations with terms which are already evolved. We never explicitly have to evolve the scalar fields equation of motion.
We begin this paper by describing a framework in which non-minimally coupled theories may be decomposed into background terms, perturbations, and an effective fluid description in Section 2. Next, the EFT approach is introduced in Section 3. The equations of state for perturbations in the EFT model under consideration are developed in Section 4, and the envisaged implementation of these results in a numerical code is laid out in Section 5. We conclude by discussing the benefits and drawbacks of this approach as well as future prospects in Section 6. Appendix A describes our conventions in detail, while Appendix B provides complete explicit formulae for the equations of state.
for dealing with a dark sector which is non-minimally coupled to gravity. We begin by laying everything out in tensorial notation with very general statements, before specializing to the context of cosmology. In subsequent sections we write down the relevant equations in component form.

General case
We organise the gravitational action into the general form where Ω is the coupling function, L m is the matter Lagrangian, and L d is the dark Lagrangian, containing all dark sector fields. We chose to work in a conformal frame where there is no dark sector coupling to matter (for scalar field theories, this is the Jordan frame). This precludes the description of models that violate the weak equivalence principle, such as models with different couplings to matter and dark matter, within this framework. After varying the action with respect to the metric, the gravitational field equations are The matter and dark sector energy-momentum tensors (EMT) are defined respectively as We now proceed by manipulating these general objects and defining various quantities so as to elucidate the effective fluid nature of the dark sector theory and how nonminimal couplings can be brought into the language of the fluid description. We define a gross dark EMT U µ ν to include the dark sector EMT and the contribution from the derivatives of the coupling function. The definition is (2.4) and the gravitational field equations (2.2) become The components of U µ ν become the "effective" energy density and pressure of the dark sector. After taking into account the Bianchi identity ∇ µ G µ ν = 0 and the conservation of energy for the matter sector ∇ µ T µ ν = 0 (which follows from the choice of frame), it follows from (2.5) that the conservation equation in terms of the gross dark EMT is given by The perturbed gravitational field equations are obtained from (2.5) as As was the case on the background, it is insightful to combine the δU µ ν and G µ ν δΩ contributions in the perturbed field equation (2.7) into the gross perturbed dark EMT, δU µ ν , which we define as This yields for the perturbed field equation. The price we pay for this clean perturbed field equation is that the conservation equation for the gross perturbed dark EMT becomes somewhat complicated. Expanding the conservation equation (2.6) to linear order in perturbations and making use of the perturbed field equation (2.9) yields where δΓ α µν = g αβ (∇ (µ δg ν)β − 1 2 ∇ β δg µν ) is the perturbation to the connection symbols from metric perturbations. We see that the equations of motion for the gross perturbed dark EMT are sourced by not only themselves and metric perturbations, but also by "normal" matter perturbations. Note that the equations of motion for "normal matter" remain unchanged, however.

Cosmology
We wish to apply the above equations to a cosmological context 2 . On a background compatible with Friedmann-Robertson-Walker (FRW) symmetries (isotropy and homogeneity), the background gross dark EMT can be written as where ρ and P are the energy density and pressure of the dark sector fluid, and u µ is the fluid four-velocity (fluid variables without labels are fluid variables of the dark fluid). The components of the gross perturbed dark EMT δU µ ν can be parameterized in a fluids language as We will refer to δρ, v µ , δP and Π µ ν as the perturbed energy, velocity field, perturbed pressure and anisotropic stress of the dark sector fluid: these are the perturbed fluid variables of the dark sector. If the dark sector contains only a single scalar field (without higher time derivatives), then there will only be two scalar components contained in these perturbed fluid variables, corresponding to the field and its time derivative. We later present relations between the perturbed fluid variables and the scalar field, as derived from an effective field theory action.
Using (2.11), the background conservation equation (2.6) yields the coupled fluid equationρ where H ≡ȧ/a is the Hubble expansion, k 0 is the spatial curvature, and overdots represent derivatives with respect to conformal time. For scalar perturbations in the synchronous gauge in momentum space, the coupled perturbed fluid equations (2.10) becomeδ where w = P/ρ is the dark sector "equation of state" and Note that for dark energy with constant equation of state, wΥ is the entropy perturbation wΓ of the dark sector fluid. Also note that these perturbed fluid equations are explicitly sourced by the perturbed matter fluid variables, due to the terms proportional toΩ. Finally, it is evident that this formulation cannot cope with models which cross the "phantom divide" of w = −1, and thus we require w > −1 for this formalism to make sense.

Effective field theory model
The perturbed fluid equations (2.14) are not closed, since we have not yet specified the form of the wΥ or anisotropic stress Π perturbations. For a given single scalar dark energy model, it should be possible to relate both of these perturbations to the scalar field perturbation and its time derivative. In the following sections, we obtain the precise forms of these equations of state for perturbations from an effective field theory (EFT) action. In this section, we briefly summarise the EFT action we employ for coupled scalar field models. Based on the Effective Field Theory of Inflation formalism [40], the Effective Field Theory of Dark Energy [28,29] has been used to describe perturbations in general single-scalar field models of dark energy. The action consists of three terms that determine the background evolution of the cosmology, and further operators that only contribute to the perturbative behavior of the scalar field. It has been demonstrated [28,30,31] that only three perturbative operators are required to completely describe perturbations in the most general single-scalar field theory with second order equations of motion, known as "Horndeski's Theory" [36] and recently rediscovered as "Generalized Galileons" [37,38]. Therefore, it is of interest to describe this action within the current formalism. However, the final term (with coefficientM 2 2 in [31]) significantly complicates the analysis we wish to perform, and thus we leave it for future work. The model without this term is sufficient to describe non-minimally coupled Kinetic Gravity Braiding models [41,42] (corresponding to the first two terms of Horndeski's theory, and a specialized combination of the second two).
The action in unitary gauge and conformal time 3 is given by The non-minimal coupling is implemented by the function Ω(τ ). The functions Λ(τ ) and c(τ ) along with Ω(τ ) are responsible for the dark sector's contribution to the evolution of the background, and the functions M 4 2 (τ ) andM 3 1 (τ ) modify the evolution of linearized perturbations. The Stückelberg trick is used to restore gauge invariance in this action, introducing a scalar field π. The relevant diffeomorphism is τ → τ + π/a, where we introduce the scalefactor for later convenience 4 .
The gravitational field equations are given by (2.2). The Λ(τ ) and c(τ ) terms are the sole contributors to U µ ν on the background. They are then combined with the derivatives of the coupling function Ω to yield the gross energy-momentum tensor U µ ν . Parameterizing the components of U µ ν as (2.11), the effective density and pressure of the dark sector fluid are [28,29] Solving these for c and Λ yield The gravitational field equations (2.5) and conservation equation (2.6) respectively become The perturbed EMT for the action (3.1) can be computed directly, and have been checked against previously derived results transformed from physical to conformal time (see, e.g., the appendices of [28]). We present the components of the perturbed field equations for different operators in the EFT, giving the gross perturbed dark EMT δU µ ν in momentum space (note the inclusion of mode functions, see Appendix A).
• The perturbed EMT for operators Ω, c and Λ is given by • The perturbed EMT for the operator M 4 2 is given by • The perturbed EMT for the operatorM 3 1 is given by Now that we know the form of the dark EMT, we can read off the expressions for the perturbed fluid variables δρ, δP , Π and θ. This yields the following.
These formulae prescribe the explicit way in which the field variables and coefficients from the EFT action combine to modify each component of the perturbed gravitational field equations.

The equations of state for perturbations
We now move on to providing the equations of state for perturbations for these coupled scalar field theories. The key thing to realize is that the perturbed fluid equations (2.14) are not closed until both wΥ and the anisotropic stress Π are specified as functions of fluid and metric components which already have evolution equations (δ, θ,ḣ and η).
In this section we show how this can be done. The first step in the calculation is to write the dark perturbed fluid variables (3.8) as (4.1) The [A IJ ] are the components of what we call the "activation matrix", and can all be calculated in terms of the coefficients in the effective action (they may also depend on wavenumber); they are given explicitly in Appendix B.1. For models that are more general than what we are considering here, there may be more non-zero components, as well as more metric perturbations appearing in the column vector on the left-hand-side. However, note that A 13 and A 23 will always be zero for single-scalar field models with second order equations of motion. In the model under consideration, it transpires that only A 11 , A 31 and A 41 have wavenumber dependence. The next step is to solve for π andπ in terms of δ and θ (and any metric perturbations). We begin with the following subexpression of the activation matrix (4.1): This can be easily inverted to give where we defined The denominator here should not vanish. We give the explicit form of the denominator as a function of terms in the effective action in (B.8).
We now need an expression forπ. Unfortunately, it's not as straightforward as just taking the time derivative of the expression forπ, as this would introduceδ anḋ θ into the expression for δP . Using the equations of motion to eliminateδ andθ would in turn introduce δP (because of the dependence of wΥ on δP ), resulting in a circular definition. Instead, note from (2.14) that the combinationδ + 3H(1 + w)θ is independent of wΥ, and thus δP . Thus, we wish to obtain two expressions forπ, and take the appropriate linear combination of them to precisely give this combination of the fluid derivatives. Taking the time derivative of (4.2) provides (4.5) In order to obtain the desired combination, (4.5) should be contracted with the row vector [1, 3H(1 + w)]. The resulting equation can be solved forπ, and yields where we defined the expression appearing in the denominator as and the expressions appearing in the numerator as We now have expressions for π,π andπ in terms of perturbed fluid and metric variables: these are (4.3a), (4.3b) and (4.6) respectively. Inserting these expressions into the relevant slots in (4.1) yields the following (schematic) expressions.
These are our equations of state for perturbations. The coefficients J i and K i are defined in terms of the A IJ in Appendix B.3. In the following section, we combine these equations with the Einstein equations to demonstrate how these equations of state can be numerically evolved.

Computational steps
We envisage the implementation of these equations in software such as CAMB [39]. In this section, we detail how we think such an implementation might work.
In the scalar sector, the metric perturbation η, the matter sector δ m and θ m , and the dark sector δ and θ are known. The matter sector behavior should be unchanged. We assume that all quantities that depend only on the background are known, or can be computed as necessary.
The perturbed Einstein equations (in synchronous gauge, as used in CAMB) are as follows.
Recall that all fluid variables without labels are those for the dark sector. From the first of these equations, we can computeḣ. The second equation then suppliesη.
In order to proceed, we should now evaluate the denominators D and E, and the numerators F and G. Expressions for these are all given in Appendix B.2. Using these quantities, we can evaluate π andπ from (4.3a) and (4.3b).
Next, we can compute the anisotropic shear stress Π. From the results of the previous section, it is given by (4.9b), although it may be easier to implement in terms of π using the expression from (B.6). Among other things, this provides the driving term in (5.1d) above.
We can now computeδ which is also given in Appendix B.2. The explicit expression for δP can be written in terms of these various quantities as Unfortunately, we do not yet knowḧ. For the moment, it makes sense to compute the quantity We now use (5.1c); some simple rearrangements providë It can be shown that the quantity in the denominator here is always positive. From here, it is now straightforward to construct δP .
We are now in a position to calculate which is the final quantity that needed to be computed. Now we have everything needed to evolve the dark sector fluid equations (2.14).

Discussion
The formalism presented here suggests a way in which the dynamics of perturbations in reasonably complicated scalar field theories can be implemented in software such as CAMB. There are a number of benefits to this approach, as well as some drawbacks. The primary benefit of this approach is that the formalism is both general and self-contained. The action (3.1) that we started from is very general in that it contains most of the generality included in Horndeski's action, and can furthermore be extended to the full theory. At the same time, the implementation of the equations of motion needs to modify only a small handful of equations in CAMB, as it already includes a simple quintessence implementation with constant equation of state. The computation of all the necessary coefficients for the formalism described here is slightly tedious, but not complicated.
A second major benefit to this formalism is that the quantities that are evolved in this formalism are physical quantities that appear in the EMT for dark energy. Thus, the evolution of physical quantities is straightforward to extract from the formalism, as is the meaning of the objects being evolved. Previous efforts to directly implement the equations of motion for the π field ran into difficulties associated with the meaning of the π field as a time displacement for the background scalar field, which is not a particularly meaningful physical quantity. Our formalism is also independent of redefinitions of the π field, and so further removes any ambiguity associated with scalar fields.
An interesting point that may be considered a benefit or a drawback depending on your point of view is that the background evolution of the model must be precomputed in order to construct the various functions of time appearing in the action, as well as their derivatives. While this requires solving the background evolution of specific models, it also allows general phenomenological models to be described, for example, by choosing various functional forms for the coefficients. As a side note, we refer the reader to [31], which relates all of the coefficients appearing in the EFT action to functions appearing in Horndeski's action.
A further interesting point is that we have developed this formalism in conformal time using synchronous gauge. This is primarily motivated by the use of these choices in CAMB, although other choices could of course be made. It may be interesting to investigate what the equations of state look like in a gauge-invariant formulation, in order to make more sense of the structures involved.
One of the drawbacks of this approach is that we cannot consider models which cross (or even touch) the "phantom divide" of w = −1. Doing so causes a number of denominators to become zero, which is clearly problematic. Although it is possible for models not to display instabilities in this regime, the description in terms of an effective fluid breaks down for such cases.
The other drawback to this approach is that we are ignoring the effects of kinetic mixing. In a number of models, including non-minimally coupled models and kinetic gravity braiding models [41,42], the scalar degrees of freedom in the metric are kinetically mixed with the scalar field. To isolate the physical scalar degree of freedom requires the diagonalisation of the kinetic matrix, which allows the speed of sound to be calculated appropriately. The formalism described here simply lumps all of the kinetic mixing terms into the dark sector effective density and velocity field perturbations. As such, it is difficult to extract the speed of sound from this formalism, and one needs to remember that the scalar metric perturbations are not completely nondynamical. Nonetheless, the evolution equations presented here are exactly equivalent to the original equations of motion. We also note that the equations of motion in terms of the original scalar field can be used to obtain such properties in a complementary manner to this approach.
We believe that this formalism will be useful for comparing models to observational data. In particular, we hope to develop a tool by which theorists can rapidly compute CMB and weak lensing spectra (and other perturbative effects) of models, without having to understand the intricacies of modifying CAMB. The next step in this direction will be to expand these results to describe the full Horndeski theory, before diving into the computational realm. and PHY-0968820. JAP is supported by the STFC Consolidated Grant ST/J000426/1.

A.1 Metric
We write the background FRW metric using conformal time as whereg ij is a three-dimensional spatial metric with no time dependence, such that the Riemann tensor of this spatial metric is given bỹ where k 0 is the curvature constant (we reserve k for wavenumber, and use K for the trace of the extrinsic curvature tensor). We use overdots to denote derivatives with respect to conformal time.
We consider perturbations in synchronous gauge. The metric with perturbations is written as The metric perturbation h ij decomposes into two scalar components as where∇ is the covariant derivative associated with the spatial metric. Evidently, g ij h ij = h.

A.2 Momentum space
We follow the conventions of Kodama and Sasaki [43]. For each wavevector k, define Y k (x i ) to be a solution of the equatioñ From now onwards, we suppress the k dependence of Y . Taking derivatives of Y , we define vector and tensor mode functions as We raise and lower indices on Y i and Y ij with the spatial metric only. This ensures that Y , Y i , and Y ij are independent of time, no matter the position of their indices. In synchronous gauge, the metric perturbation decomposes as Following Ma and Bertschinger [44], we define η bỹ Using this definition, the metric perturbation in synchronous gauge becomes In the limit of no spatial curvature, the mode function Y = exp(i k · x) yields Ma and Bertschinger's Eq. (4).

A.3 Energy-momentum tensors
On an FRW background, the background EMT for any sector can be written as where u µ is the velocity vector of the fluid, u µ = (1/a, 0, 0, 0), and P and ρ are the background pressure and energy density. Quantities that belong to the matter sector are given a subscript m , while quantities associated with the dark sector have no subscripts. Again following Kodama and Sasaki, we describe scalar perturbations of the EMT in momentum space using the following general decomposition (we use this decomposition for both the matter and dark sectors).
Note that indices on the modefunctions have been raised only with the spatial metric. We can convert v into the θ variable of Ma and Bertschinger by θ = kv. Similarly, Ma and Bertshinger's σ variable for the anisotropic shear stress is related to Π by σ = 2Πw/3(1 + w). We find it more convenient to work with θ = v/k = θ MB /k 2 , which we use throughout the rest of this document. Note that we use δ = δρ/ρ as the fractional density perturbation.

B Explicit forms of coefficients
In this appendix, we present the explicit expressions for coefficients appearing in our formulae that we kept largely hidden in the main body of the paper.

B.1 Mapping from the effective action to activation matrix
We begin by presenting the components of the activation matrix (4.1) in terms of the quantities derived from the effective action (3.8). Note that the inverse powers of the scalefactor are largely due to using conformal time. Common elements to a number of these quantities are The perturbed density is given by where the activation matrix components are given by The velocity divergence field is given by where the relevant activation matrix components are The perturbed pressure is Here we denoted (M 3 1 ) · ≡ ∂ τM

B.2 Frequently used expressions
We now show how the explicit expressions for the A IJ , which we gave as functions of terms in the effective action in (B.3 -B.6), combine to yield some of the commonly appearing expressions in the coefficients in δP and Π. We begin with explicit expressions that appear in denominators, D and E, defined in (4.4) and (4.7) respectively. We find These expressions appear in the denominators of π,π andπ, and as such, must be nonvanishing. In particular, since it is known from analysing the EFT construction that the stability of the theory requires c + 2M 4 2 > 0, we expect both of these expressions to be positive.
We now move on to the explicit expressions for the terms F and G, defined in (4.8), as well asȦ 14 , which appear in numerators. These are by far the most complicated expressions we have to deal with. Beginning with G, we break it up into

B.3 Complete expressions
For completeness' sake, here we include the full expressions for Π and δP . They can be written schematically as δP = J 1 δ + J 2 θ + J 3ḣ + J 4ḧ + J 5 [δ + 3H(1 + w)θ] (B.13) This is before rearranging the space-space trace equation to obtainḧ. Note that the term with coefficient J 5 can be expressed in terms of δ, θ and Π, but it seems unnecessary to insert the explicit expression (B.12) that displays no dependence on A IJ apart from that implicitly included in Π.
We find that the coefficients J i are given in terms of the A IJ as