The Effective Fluid approach for Modified Gravity and its applications

In this review we briefly summarize the so-called effective fluid approach, which is a compact framework that can be used to describe a plethora of different modified gravity models as general relativity (GR) and a dark energy (DE) fluid. This approach, which is complementary to the cosmological effective field theory, has several benefits as it allows for the easier inclusion of most modified gravity models into the state-of-the-art Boltzmann codes, that are typically hard-coded for GR and DE. Furthermore, it can also provide theoretical insights into their behavior, since in linear perturbation theory it is easy to derive physically motivated quantities such as the DE anisotropic stress or the DE sound speed. We also present some explicit applications of the effective fluid approach with $f(R)$, Horndeski and Scalar-Vector-Tensor models, namely how this approach can be used to easily solve the perturbation equations and incorporate the aforementioned modified gravity models into Boltzmann codes so as to obtain cosmological constraints using Monte Carlo analyses.


Introduction
At the end of the previous century statistically significant evidence from observations of Type Ia supernovae (SnIa) revealed that the Universe is currently undergoing a phase of accelerated expansion [1,2]. This accelerated expansion is typically attributed to a cosmological constant Λ, which in addition to the standard Cold Dark Matter (CDM) scenario, can alleviate several deficiencies of the latter [3]. Together these two components form the ΛCDM model, which has been found to be in excellent agreement with recent cosmological measurements [4,5]. Despite that, still the cosmological constant provides a problem for theoretical physics due to the large discrepancy between the predicted and observed values of Λ [6,7].
Since then, the ΛCDM model has become the standard cosmological model as it provides the best description of the observations on cosmological scales [4,5,8]. As a consequence, many alternative explanations have emerged and for the most part there are two main approaches. The first one, which is also more concretely based on high energy physics is the one where Dark Energy (DE) models [9], due to as yet unobserved scalar fields, dominate the energy budget of the Universe at late times, and if their mass is sufficiently light they also lead to an accelerated expansion [10,11]. The second approach is based on the assumption that covariant corrections to the Theory of General Relativity (GR) can alter gravity, usually dubbed modified gravity (MG), at sufficiently large scales [12]. However, several cosmological probes in extra-galactic scales are in good agreement with GR [13,14].
Overall, both approaches with either DE or MG models provide realistic and plausible explanations for the accelerated expansion of the Universe at late times. Furthermore, both kinds of models can also fit the cosmological observations at the background level equally well to the ΛCDM model, as they can always go arbitrarily close to the cosmological constant. Thus, these models are in principle degenerate at the background level, despite laborious efforts to arXiv:2212.12768v3 [astro-ph.CO] 2 Feb 2023 break the degeneracies with model independent approaches [15,16]. Fortuitously, the recent discovery of gravitational waves by the joint LIGO/Virgo collaboration [17] has allowed the community to rule out several MG models [18][19][20][21][22][23][24][25][26][27].
One of the main remaining classes of MG models is the so-called f (R) models [28][29][30][31]. Also in this case the background evolution of the Universe is degenerated with DE models with a redshift dependent equation of state w(z) [32][33][34][35]), albeit the linear theory perturbations are in general vastly different and have a particular time and scale dependence [36]. This is particularly important as in general the DE perturbations have a definite effect on the growthrate of matter perturbations and the so-called growth-index γ [37]. However, tat this point current data do not favor any particular model f (R) model [38,39].
This discussion obviously makes it clear that the perturbations of MG models are of great importance and several approaches exist in the literature, e.g. see Refs. [33,34,36,40,41,[41][42][43][44][45][46][47][48][49][50][51][52][53][54][55][56]. However, even though many authors consider MG models, they still sometimes fix the background to that of the flat ΛCDM model, as for example was done in Ref. [57], but model the evolution of the Newtonian potentials via two functions µ(a, k) and γ(a, k) which take into account possible deviations from GR. These functions have been implemented in modified versions of the codes CAMB [58] and CLASS [59]. Even though these parameterizations are only valid at late times, also new parameterizations which are valid at all times have appeared for MGCAMB, see Ref. [60], and MGCLASS, see Ref. [61], which is build to be compatible with parameter estimation codes. A definite flaw of the approach in early versions of MGCAMB was that the background expansion was fixed to that of the ΛCDM model, even though it is known that for f (R) models it is quite different, as is for example the case for the Hu-Sawicki model [40].
A complimentary approach was carried out in Ref. [62], where the author studied perturbations of f (R), which are degenerate to the ΛCDM at the background level, by utilizing the full set of covariant cosmological perturbation equations and modifying the publicly available code CAMB, in a new version called FRCAMB. Furthermore, an unreleased extended version of the aforementioned code with arbitrary background expansion rates was created by Ref. [63].
A totally different approach to modeling DE and MG models have been proposed in the form of the Effective Field Theory (EFT) [64] and applied to cosmology in Ref. [65] in the form of a code called EFTCAMB. The main aadvantage of this approach is that it does not utilize any approximations, albeit the mapping of specific MG and DE models into the EFT formalism is somewhat complicated in most cases. Some of the aforementioned codes, namely MGCAMB and EFTCAMB, where used by the Planck Collaboration in Ref. [66] to derive cosmological constraints of the MG and DE models. Overall, no conclusive and statistically significant evidence for models beyond ΛCDM was found.
Finally, another interesting approach was followed by Ref. [67] where the authors proposed the so-called Equation of State (EOS) approach for perturbations which maps f (R) models to a DE fluid at both the background and linear perturbation order [42,49], see also [68][69][70]. The EOS approach has been implemented in a modified version of the code CLASS [59] in Ref. [71], albeit the problem with this approach is that the interpretation of the perturbation variables is not clear.
In this work we will also map the MG models as a DE fluid by utilizing the DE equation of state w(a), the sound speed c 2 s (a, k) and the anisotropic stress π(a, k), as these variables are enough to describe any MG fluid at the background and linear order of perturbations [72]. This has the advantage that it makes the comparison with popular DE models such as quintessence (w(a) ≥ −1, c 2 s = 1, π(a, k) = 0) and K-essence (w(a), c 2 s (a), π(a, k) = 0) straightforward. This is clearly important, as in general in MG models the anisotropic stress is non-zero π(a, k) = 0, whereas in standard quintessence π(a, k) = 0, so that any statistically significant deviation of the anisotropic stress from zero, would be a smoking gun for MG models [72,73].
In order to simplify the analysis of the perturbation equations, the quasi-static and subhorizon approximations are frequently utilized. The former is based on the observation that in matter domination the Newtonian potentials are mostly constant, thus terms in the linearized Einstein equations with potentials with time derivatives can for the most part be safely neglected. The latter, is based on the observation that only perturbations with wavelengths shorter than the cosmological horizon are important. Some of the previous codes, i.e. FRCAMB, EFTCAMB, CLASS_EOS_FR did not in fact apply the sub-horizon approximation to the perturbation equations, albeit the quasi-static approximation has been studied extensively and implemented in MGCAMB [45,74].
It has been argued, see for example Ref. [74], that the quasi-static approximation breaks down outside the DE sound-horizon k is the physical Jeans scale, rather than outside the cosmological horizon. However, in the aforementioned analysis the anisotropic stress was neglected and a constant DE c 2 s was utilized, both assumptions being unrealistic in general.
This approach allows us in general to discriminate between traditional DE and MG models, as they both have vastly different predictions for the equation of state w(a), the sound speed c 2 s (a, k), and the anisotropic stress π(a, k). The last two quantities are in particular important as they in general may leave observable traces in the large scale structure (LSS), the Cosmic Microwave Background radiation (CMB) and Galaxy Counts (GC) [36,75]. Furthermore, an important aspect of the anisotropic stress is that in general it can stabilize the growth of matter perturbations in cases where that would be not possible [73,[75][76][77], while the sound speed affects the clustering of matter perturbations [78][79][80]. Both of these effects are crucial, as they can be used to break the parameter degeneracies between the models [81,82].
Even though the ΛCDM model seems to be in good overall agreement with observations [4,5], this might easily change when forthcoming galaxy surveys like Euclid, DESI and stage IV CMB experiments arrive. Furthermore, there also seem to remain some more issues with cosmological data such as direct Hubble constant measurements, weak lensing data, and cluster counts, where different aspects of DE models or MG models could be important [75,[83][84][85][86][87][88][89][90][91], thus the effective fluid approach would be particularly useful.
This review is organized as follows: in Sec. 2 we present the theoretical framework of the effective fluid approach and its application to f (R), Horndeski and Scalar-Vector-Tensor models. Then, in Sec. 3 we present several concrete applications of our approach, namely designer Horndeski models, the numerical solutions of the perturbation equations and the necessary modifications to Boltzmann codes so that comparison with the CMB data and Monte Carlo analyses can be made. Finally, in Sec. 4 we summarize our the effective fluid approach and present our conclusions.

Theoretical framework
Here we now describe the theoretical framework necessary to illustrate the Effective Fluid approach, especially related to linear order of perturbation theory. On large scales the Universe is homogeneous and isotropic, thus it can be described at the background level by a Friedmann-Lemaître-Robertson-Walker (FLRW) metric. In order to describe the large scale structure of the Universe, we need to consider the perturbed FLRW metric, which in the conformal Newtonian gauge is given by: given in terms of the conformal time τ defined via dτ = dt/a(t) and we also follow the notation of Ref. [92]. 1 On large cosmological scales, where the average density of the matter particle species is very low with respect to terrestrial ones, namely on the order of ρ ∼ ρ cr = 1.8788 × 10 −26 h 2 kg m −3 , this means we can assume the matter species can be described as ideal fluids with an energy momentum tensor where ρ, P are the fluid density and pressure, while U µ = dx µ √ −ds 2 is its velocity four-vector given to first order by U µ = 1 a(τ) (1 − Ψ, u), which satisfies U µ U µ = −1 and we have defined u =˙ x andḟ ≡ d f dτ . Then, the elements of the energy momentum tensor are given to linear order of perturbations by: whereρ,P are background quantities, functions of time only. On the other hand, the perturbations of the fluid's density and pressure are given by δρ, δP and are functions of ( x, τ). Finally, is an anisotropic stress tensor. In the context of GR we find that the perturbed Einstein equations are given in the conformal Newtonian gauge by [92]: where the velocity is given by θ ≡ ik j u j and the anisotropic stress by (ρ +P)σ ≡ −(k ikj − 1 3 δ ij )Σ ij , which as mentioned is related to the traceless part of the energy momentum tensor via Σ i j ≡ T i j − δ i j T k k /3 and where k j andk j are the wavenumbers and unit vectors in Fourier/k-space of the perturbations.
To find the evolution equations for the perturbation variables we use the energy-momentum conservation T µν ;ν = 0 conservation given by the Bianchi identities in GR as: 1 In this review our conventions are: (-+++) for the metric signature, the Riemann and Ricci tensors are given by where we have defined the rest-frame sound speed of the fluid c 2 s ≡ δP δρ and its equation of state parameter w ≡P ρ . After eliminating θ from Eqs. (10) and (11) we find a second order equation for δ [75]:δ where the dots (. . .) stand for a complicated expression, while the anisotropic stress of the fluid is defined as π ≡ 3 2 (1 + w)σ. As can be seen, the k 2 term behaves as a source driving the perturbations of the fluid, however as the potential scales as Ψ ∼ 1/k 2 in relevant scales (due to the Poisson equation), the dominant terms are only the sound speed and the anisotropic stress. Thus, we can define a parameter, namely an effective sound speed which controls the stability of the perturbations, as [75]: Clearly, c 2 s,eff not only characterizes the propagation of perturbations, it also defines the clustering properties of the fluid on sub-horizon scales, see Ref. [75]. In general, the sound speed c 2 s can be both time and scale dependent, i.e., c 2 s = c 2 s (τ, k), e.g. as noted in Ref. [93], the sound speed in small scales for a scalar field φ (in the conformal Newtonian gauge) is given by where m φ is the mass of the scalar field.
However, the sound speed is equal to unity only in the scalar field's rest-frame. see for example Chapter 11.2 of Ref. [93]. In f (R) theories, as they in fact can be viewed as a nonminimally coupled scalar field in the Einstein frame, see for example Ref. [74,94], the sound speed is also scale dependent, when we are not in the rest frame of the equivalent DE fluid.
In the end, it is common to use the scalar velocity perturbation V ≡ ik j T j 0 /ρ = (1 + w)θ, instead of the fluid velocity θ as the former can remain finite when the equation of state −1, see for example Ref. [95]. Then, Eqs. (10)-(11) become where a prime means a derivative with respect to the scale factor a, while H(t) = da/dt a is the cosmic-time Hubble parameter.

f (R) models
The simplest application of the effective fluid approach in theories beyond GR, is of course in the context of f (R) models. These can of course be studied directly, as was done in Ref. [36] or as an the effective DE fluid [67]. Specifically, the modified Einstein-Hilbert action is given by: where L m is the matter Lagrangian and κ ≡ 8πG N . Varying the action of Eq. (16) with respect to the metric, we obtain the field equations [36]: where we have defined F ≡ f (R), G µν is the usual Einstein tensor and T (m) µν is the energymomentum tensor of the matter fields.
Bu moving all the modified gravity contributions to the right hand side, we can rewrite the field equations as the Einstein equations being equal to the sum of the energy momentum tensors of the matter fields and that of an effective DE fluid [49]: where we have defined As f (R) theories are also diffeomorphism invariant, one can show that the effective energy momentum tensor given by Eq. (19) also satisfies a conservation equation: Writing the theory in this way implies that the Friedmann equations are the same as in GR [92]: albeit with the addition of an effective DE term on the right hand side described by a density and pressure: where H =˙a a is the conformal Hubble parameter. From Eqs. (23) and (24) we can define an effective DE equation of state for the f (R) models as: which agrees with the one given in Ref. [36].
Using the effective energy momentum tensor of Eq. (19) we can define the effective DE pressure, density and velocity perturbations as: while the difference of the two Newtonian potentials Φ and Ψ is given by Thus, the anisotropic stress is given by [92] 2.1.1. The quasi-static and sub-horizon approximations As can been seen, the expressions for the DE perturbations given by Eqs. (26)-(30) are somewhat cumbersome and can be significantly simplified, without much loss of accuracy, by using the sub-horizon and quasi-static approximations. The former implies that only the modes deep in the Hubble radius (k 2 a 2 H 2 ) are important, while the with the latter we neglect terms with time derivatives. As an example we find that the perturbation of the Ricci scalar is where in the last line we applied the two approximations. Then from the perturbed Einstein equations we find the modified Poisson equations [36]: where G eff and Q eff are both equal to unity in GR and are given by [36]: where we have set Alternatively we can also write the Poisson equation for Φ in the effective fluid approach, where we can to introduce the DE density ρ DE and get: which implies thatρ which allows us to determine the evolution of the DE density perturbation directly.
With these approximations we can also directly derive a second order differential equation, when ignoring neutrinos, for the time evolution of the matter density contrast [36]: where derivatives with respect to the scale factor a are denoted by a prime . Finally, we can also define the DE anisotropic parameters Applying the sub-horizon and quasi-static approximations, and using the Poisson equations, we can estimate the effective density, pressure and velocity perturbations of the DE fluid as: while the DE anisotropic stress parameter π DE is given by As can be seen, the DE anisotropic stress can also be written i nthe more general form as where we have defined the functions f 1 (a) = , which are reminiscent of Model 2 in Ref. [75]. Now using Eqs. (42) and (43), we find that effective DE sound speed of the fluid is given at this level of the approximation by and the DE effective sound speed is As can be seen from the previous expressions, for the ΛCDM model ( In the case of f (R) models, such as the Hu & Sawicki (HS, hereafter), where the DE equation of state w DE crosses w DE (a) = −1, one would expect singularities to appear because of the 1 + w term in the denominator in Eq. (11) [96]. However, we can absorb the 1 + w term by introducing V DE = (1 + w DE )θ DE , as this combination remains finite for well-behaved f (R) models, as seen by inspecting Eq. (44).
As a final remark, it should be noted that the effect DE fluid described here and in Eq. (19), in fact violates the energy conditions of GR [98], which can be expressed via the DE density and pressure: where NEC, WEC, DEC and SEC stand for the null, weak, dominant and strong energy conditions. As the inequalityρ DE ≥ 0 still holds, then the NEC, WEC and DEC conditions can be mapped equivalently to the constraint w DE ≥ −1. However, as can be seen in Fig. 1 in the case of the HS model the NEC, WEC and DEC are violated for redshifts z 1.65.

Horndeski models
The most general Lorentz-invariant extension of GR in four dimensions with a nonminimally coupled scalar field and second order equations of motion is the so-called Horndeski theory [99]. This theory contains several free functions that, in appropriate limits, reduce to several well-known DE and MG models. However, the recent observation of a binary neutron star and its accompanying optical counterpart, has produced an amazingly tight constraint on the speed of propagation of gravitational waves (GWs) [100] −3 · 10 −15 ≤ c g /c − 1 ≤ 7 · 10 −16 .
This constraint implies that the functional forms of two of the free Horndeski functions were then limited to be [20] G 4X ≈ 0, G 5 ≈ const., as seen from the formula for the GW speed of propagation [101] Thus, in the case of the effective fluid approach, we will only focus on the remaining parts of the Horndeski Lagrangian, in particular where we have defined and φ is a scalar field, X ≡ − 1 2 ∂ µ φ∂ µ φ is its kinetic term, and φ ≡ g µν ∇ µ ∇ ν φ; K, G 3 and G 4 are free functions of φ and X. 2 Finally, we also assume L m includes the matter fields.
These specific terms correspond to different dynamics, for example K(φ, X) contains the k-essence and quintessence theory, however does not contribute to the perturbations [24]. On the other hand, the term G 3 (φ, X) contains the so-called kinetic gravity braiding (KGB) models, with G 3X = 0 corresponding to mixing of the kinetic term of the scalar and the metric, with G 3 (φ) only modifying the background as a dynamical DE. The last term, namely G 4 , in in fact the one that contains the non-minimal coupling of the scalar to the Ricci curvature, and contains most scalar-tensor type theories. To give some concrete examples, the action (52) can easily be seen to reduce to the following subclasses: • f(R) theories: These are equivalent to a non-minimally coupled scalar field written as [102] where φ ≡ f ,R √ κ has units of mass and we have set f ,R ≡ d f dR .
• Brans-Dicke theories: These are the archetype of a scalar-tensor theory, with where V(φ) is the potential and ω BD is the well-known Brans-Dicke parameter [103]. • Kinetic gravity braiding: These models contain a mixing of the scalar and tensor kinetic terms [104] and are given by • The non-minimal coupling model: This is given in terms of a coupling constant ζ as In the case of inflation, the Higgs-like inflation model is given by ω(φ) = 1 and V(φ) = λ φ 2 − ν 2 2 /4.
• Cubic Galileon: The well-known models are given by [105] Then, varying the action of Eq. (52) with respect to the metric and the scalar field φ we can obtain the equations of motion. First, doing a variation we find [101]: from which the field equations follow. The gravitational field equations are: where we have set and T

(m)
µν is the energy-momentum tensor of matter. When K = G 3 = 0 and G 4 = 1 2κ we can see that the Eqs. (70) reduce to those of GR. Similarly, the equations of motion of the scalar field are given by where again we have set Here it should be noted that while it would seem that the term ∇ µ J i µ leads to higher than second-order derivatives, in fact it does not as was first noted in Ref. [101]. In fact this is due to that the commutations of higher derivatives can be shown to cancel out as Doing some algebra it is possible to show that the scalar field equation (74) can be reduced to [97] − For the sake or brevity, in what follows we will denote by X the kinetic term of the scalar field evaluated at the background and by δX its linear order perturbation.

Background expansion
At the background level, assuming an unperturbed flat FRLW metric, it is easy to show that the modified Friedmann equations are given by where we have defined the quantities Gathering all the terms we can write down the explicit equations as Again, in the limit of K = G 3 = 0 and G 4 = 1 2κ , these reduced to the standard Friedmann equations as expected. By rearranging and collecting the terms in Eqs. (91)-(92) we can define the density parameter of an effective DE fluid and its effective pressure: thus reducing the modified Friedmann equations Eqs. (91)- (92) to their traditional GR form Using Eqs. (93)-(94) then allows us to also define the effective DE equation of state as . (97) We can also write down the explicit scalar field equation (82) as [106] and by setting it is possible to rewrite the scalar field equation (74) as a conservation equation from which it is easy to deduce the existence of a Noether symmetry under constant shifts of the field φ → φ + c, given by [104] using the fact that that X = 1 2φ 2 , then we find that the conserved quantity is given by so that the scalar field equation reduces to a much simpler forṁ When P φ = 0, then the solution is simply for a constant J c . Here we can also classify some particular subcases: first, when J c = 0 then the scalar field is on the attractor solution, while when J c = 0 then the system is not on the attractor and new dynamics may arise [97].

Linear perturbations
By using the perturbed FLRW metric of Eq. (1) with the field equations (70) we can obtain the linear theory predictions for the perturbations [107,108] Again, in the case when K = G 3 = 0 and G 4 = 1 2κ , we can see that Eqs. (106)-(109) reduce to the GR limit given by Eqs. (6)- (9) with no anisotropic stress. On the other hand, if we consider Eq. (82) we similarly find the perturbed equations of motion for the scalar field where the full expressions for the variables A i , µ, M, ν, B i , C i and D i can be found in Appendix B of Ref. [109].

The effective fluid approach for Horndeski models
Following a similar approach as in the case of the f (R) models, we can now apply the effective fluid approach to the Horndeski models as well. While we already showed this for the background effective DE density and pressure given by Eqs. (93)- (94), so in what follows we also do the same for the linear theory of these models, under the sub-horizon and quasi-static approximations.
Obviously, the first step is to define an effective DE fluid by moving all MG contributions to the right hand side of the field equations and define the DE effective energy-momentum tensor T DE µν . Doing so with the gravitational field equations (70), we find: By considering the decomposition of the T DE µν tensor into its components, given by Eqs. (3)-(5), we can extract the expressions for the DE effective perturbations in the pressure, density, and velocity. Doing so, we find the latter have the general structure: where the dots (. . .) indicates long expressions.
Following the same procedure as before and using the sub-horizon and quasi-static approximations to Eqs. (106), (108), and Eq. (110) we find [109] where the parameters G eff and Q eff are Newton's effective constant and the lensing variable We also define the DE anisotropic stress parameters as which are in agreement with the ones found in Ref. [107]. For these models Eq. (39) is still valid, albeit with G eff given by Eq. (121).

Scalar-Vector-Tensor models
An interesting extension of Horndeski models is Scalar-Vector-Tensor (SVT) models, which also include a vector degree of freedom and also include generalised Proca theories [110]. The most general SVT Lagrangian is given by [111] L ≡ where the Lagrangians L ST i include the scalar-tensor terms, while the Lagrangians L SVT i contain the scalar-vector-tensor terms.
The SVT models in fact contain a scalar field ϕ, a vector field A µ and the gravitational field g µν , which are related to each other via the L SVT i terms, given by where for the sake of brevity we defined g ξ ≡ ∂g ∂ξ , denoting the derivative of a g with respect to the scalar ξ. In the previous equations, we also defined the kinetic terms and couplings between the scalar and vector fields as As usual, the antisymmetric tensor F µν and its dualF µν are related to A µ via where ε µναβ ≡ µναβ √ −g , µναβ is the Levi-Civita symbol. Furthermore, the Lorentz invariant quantities can be constructed from F µν as which vanish when A µ → ∇ µ π, where π is a scalar field. Moreover, the symmetric tensor S µν is related to A µ as while the tensors M and N are given by where Finally, the L µναβ tensor is given by As always, the scalar-tensor interactions are contained in the Horndeski theory discussed in the previous section: where the terms f i ,f i , h 5i ,h 5i , and G i are free functions. It should be noted that even though this theory contains several free and a priori undetermined functions, in practice it has been significantly constrained by the recent GW discovery [17].

The effective fluid approach for SVT theories with non-vanishing anisotropic stress
Following the same methodology as before, the quasi-static and sub-horizon approximations can be applied to the SVT model in order to determine the DE fluid parameters in the effective fluid approach. These where found by Ref. [111] to be where the coefficients Z i (i = 1, . . . , 15) are given in Appendix E of [111]. As in previous cases, the sound speed given in Eq. (144) does not on its own determine the stability of subhorizon perturbations, but in fact the relevant quantity is the effective sound speed given, as in the previous cases, by the difference of the DE sound speed and a term proportional to the anisotropic stress π = 3 2 (1 + w)σ. The effective sound speed again in this case is defined as [75]

Designer Horndeski
One particularly interesting model to demonstrate the effective fluid approach is a class of designer Horndeski parameterization, discovered in Ref. [109] and rediscovered later in Ref. [112]. These models have a background expansion exactly equal to that of the ΛCDM model, but also have different perturbations, and are particularly useful in searches for deviations from ΛCDM [35,97].
For example, in a Horndeski model with just the G 2 and G 3 terms, i.e. of the KGB type, we can use the modified Friedmann equation and the scalar field conservation equation to solve for the unknown functions, while demanding that H(a) corresponds to the of the ΛCDM model. In the previous equations J c is a constant which quantifies our deviation from the attractor, as in the case of the KGB model [106]. Solving Eq. (146) and (147) for (G 3X (X), K(X)) yields [109]: where we have written H = H(X), i.e. in terms of the kinetic term. in fact, Eqs. (148) give us a whole family of designer models that behave as ΛCDM at the background level but have different perturbations. For example, assuming that X = c 0 H(a) n , where c 0 > 0 and n > 0, we find [109]: As can be seen, this designer model, designated as HDES hereafter, has a nice and smooth limit to the ΛCDM model and it also recovers GR when J c ∼ 0.

Numerical solutions of the perturbation equations
Here we now present the numerical solutions of the perturbation equations in the effective fluid approach, using as an example the HDES model, given by Eq. (149). In particular: • First, we consider the numerical solution of the full system of equations given by Eqs. (106)- (109), which we call "Full-DES". • Second, we consider the numerical solution of the effective fluid approach given by Eqs. (14)- (15), which we call "Eff. Fluid". • Third, we consider the numerical solution of the growth factor equation (39) using the appropriate expression for G eff , which we call "ODE-Geff". • Finally, we also consider the ΛCDM model.
As a concrete example, we setc 0 = 1, n = 2, Ω m,0 = 0.3, k = 300H 0 and σ 8,0 = 0.8, unless otherwise specified. Then we show the evolution of the growth-rate parameter f σ 8 (z) for the HDES model on the left panel of Fig. 2. Specifically, we show the "Full-DES" brute-force numerical solution, the effective fluid approach, the ΛCDM model and the numerical solution of the G eff equation and as can be seen, the agreement between all approaches is excellent. On the other hand, in the right panel of Fig. 2 we show the percent difference between the

Modifications to CLASS and the ISW effect.
Finally, we now show how the effective fluid approach can be implemented into a Boltzmann code, like CLASS [59]. Our modifications to the code are denoted as EFCLASS [97,109], while we also compare with the hi_CLASS code [114], which solves the full set of dynamical equations, but at the cost of significantly more complicated modifications.
In our case however, the modifications for the effective fluid approach are much easier, as we only require the DE velocity and the anisotropic stress [97,109]. In the particular case where we consider the HDES model, then there is the further simplification that the anisotropic stress π DE is also zero, as can be seen from Eq. (109), since G 4φ = 0. Thus, to modify CLASS we only require the expression for the DE velocity, which for n = 1 is given by Next, we show the results of applying Eq. (150) to the CLASS code and comparing with hi_CLASS. This is shown in the left panel of Fig. 3, where the low-multipoles of the TT CMB spectrum for a flat universe with Ω m,0 = 0.3, n s = 1, A s = 2.3 · 10 −9 , h = 0.7 and (c 0 ,J c , n) = (1, 2 · 10 −3 , 1) can be seen. Our EFCLASS code is denoted by the green line, hi_CLASS is given by the orange line, while the blue line corresponds to the ΛCDM model. As can be seen, the agreement between EFCLASS and hi_CLASS is remarkable and even with our simple modification is roughly ∼ 0.1% accuracy across all multipoles (shown on the right panel of Fig. 3).
Finally, we also compare our modifications of the effective fluid approach with a direct calculation of the Integrated Sachs-Wolfe (ISW) effect. Specifically, the temperature power spectrum is given by [115]:  where I ISW (k) is a kernel that depends on the line of sight integral of the growth and a bessel function, while P ζ (k) is the primordial power spectrum, given by the primordial power spectrum times a transfer function [97,115] where A s is the primordial amplitude, k 0 is the pivot scale and T(k) is the Bardeen, Bond, Kaiser and Szalay (BBKS) transfer function [116]. We show in Fig. 4 a comparison between CLASS and hi_CLASS for the ΛCDM model (left) and the HDES models (right), for the same parameters as in Fig. 3. Overall, there is good agreement at all multipoles, except for = 2, as the BBKS formula is only 10% accurate on large scales or equivalently, small multipoles.

Conclusions
In this review we briefly presented the so-called effective fluid approach, which is a framework that allows for treating modified gravity models as GR with an ideal DE fluid, described by an equation of state, a sound speed, an anisotropic stress and DE pressure, density and velocity perturbations. While in general MG models have very complicated evolution equations, they are significantly simplified under the effective fluid approach with the addition of the quasi-static and sub-horizon approximations, as they can be written in terms simple DE fluid equations. However, it should noted that in general it is not possible to discriminate modified gravity models from GR plus DE exotic fluids, as one may always move the modifications of gravity to the right hand side of the field equations and rewrite them as the Einstein equations with extra matter/DE contributions.
Thus, the main advantage of this approach, as highlighted in this work, is that it allows us to easily including most MG models in the Boltzmann codes, as the latter are typically finetuned and hard-coded for GR and a DE fluid. Here, we presented the main formalism and the results for the DE fluid quantities for several MG models, including the f (R), Horndeski and Scalar-Vector-Tensor models, in all of which we presented the critical quantities that describe the evolution of the perturbations and the effective DE fluid equations.
We also described some specific applications, i.e. the numerical solutions of the fluid equations for growth factor, the designer Horndeski model and how our approach can be implemented in the Boltzmann codes. In the case of the HDES model we found that with our effective fluid approach and just a simple modification, the agreement between our modified EFCLASS code and the more complicated, but exact, hi_CLASS code, is roughly ∼ 0.1% accuracy across all multipoles, as shown on the right panel of Fig. 3.
Finally, we should also stress that an open point of debate in the community is the proper application of the quasi-static approximation, especially in more complicated models, as is extensively discussed in Ref. [117]. Even though this issue does not affect our effective fluid approach, it is a point that should be addressed prior to the advent of the next-generation surveys, in order to avoid unwanted theoretical errors in the predictions. Still, as demonstrated in this work, the effective fluid approach can provide a powerful framework, covering most viable MG models and allowing for a simple (and educational) way to modify the Boltzmann codes, which are necessary in data analyses. Acknowledgments: In this section you can acknowledge any support given which is not covered by the author contribution or funding sections. This may include administrative and technical support, or donations in kind (e.g., materials used for experiments).

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: