On validity of the quasi-static approximation in scalar-tensor theories

The discovery of cosmic acceleration motivated extensive studies of dynamical dark energy and modified gravity models. Of particular interest are the scalar-tensor theories, with a scalar field dark energy non-minimally coupled to matter. Cosmological constraints on these models often employ the quasi-static approximation (QSA), in which the dynamics of the scalar field perturbations is proportional to the perturbation in the matter density. Using the QSA simplifies the physical interpretation of the phenomenology of scalar-tensor theories, and results in substantial savings of computing time when deriving parameter constraints. Focusing on the symmetron model, which is a well-motivated scalar-tensor theory with a screening mechanism, we compare the exact solution of the linearly perturbed field equations to those obtained under the QSA and identify the range of the model parameters for which the QSA is valid. We find that the evolution of background scalar field is most important, namely, whether it is dominated by the Hubble friction or the scalar field potential. This helps us derive a criterion for the symmetron model, but same argument can be applied to other scalar-tensor theories of generalized Brans-Dicke type. We consider two scenarios, one where the scalar field is only coupled to dark matter and where it couples to all of the matter.


I. INTRODUCTION
The Λ Cold Dark Matter (ΛCDM) model provides an acceptable fit to the data, yet many other dark energy (DE) and modified gravity (MG) models have been proposed and studied in the literature [1][2][3][4] with the aim to explain the observed cosmic acceleration [5,6]. The primary motivation for these studies comes from the extreme fine-tuning required to reconcile the small observed value of Λ with the large vacuum energy density predicted by particle physics [7,8]. The necessity to postulate yet another dark component, CDM, the evidence for which is purely gravitational, further adds to the interest in possible extensions of General Relativity (GR).
Modifications of GR typically involve new degrees of freedom [3], such as a scalar field which can, along with being a source of dynamical DE, mediate a force between matter particles often referred to as the "fifth force". There are very tight constraints on the presence of such a fifth force acting on matter that is part of the standard model of particle physics [9]. Hence, any modified gravity model with a scalar field non-minimally coupled to all matter that is designed to have cosmological implications must also include a screening mechanism that would suppress the fifth force inside the Solar System [4,[10][11][12][13]. On the other hand, if the scalar field coupled only to DM, then the Solar System tests of gravity would not apply and the screening would not be necessary. While it is unnatural for a light scalar field playing the role of DE not to couple to all matter [14], once one allows for the existence of a dark sector, one should allow for the possibility that DM and DE have their own special interaction -a possibility extensively studied in the literature [15][16][17][18]. Ultimately, we would like to use observational data to differentiate between the dark-matter-only and all-matter coupled DE, although this can be challenging in practice [19].
Deriving cosmological constraints on scalar-tensor theories often rely on the quasi-static approximation (QSA), where the scalar field is assumed to be tracking the evolution of the matter fluid to which it is coupled. Using the QSA not only helps to simplify the numerical implementation, potentially saving substantial computing time required to resolve rapid scalar field oscillations, but also makes it possible to relate parameters of the scalar-tensor theory to the phenomenological parameters, such as the effective Newton's constant and the gravitational slip, used to constrain modifications of GR on cosmological scales with the help of publicly available codes such as MGCAMB and MGCLASS [20][21][22][23]. Understanding the range of applicability of the QSA is, therefore, important for the practical implementation and meaningful interpretation of model-agnostic constraints on scalar-tensor theories.
In this paper, we consider the symmetron model [13,24], which is a well-motivated scalar-tensor theory with a screening mechanism and a renormalizable potential. The potential of the symmetron field acts as the effective cosmological constant which causes the accelerated expansion of the universe at late times. The matter particles feel the presence of a fifth force at late time, due to the non-minimal coupling with the scalar field. When matter density is higher than a certain critical value, the coupling vanishes and one recovers GR. This ensures that GR is recovered at early times and in dense regions such as our galaxy.
In what follows, we provide a comprehensive analysis of the evolution of the background cosmology and linear perturbations in the symmetron model, with and without using the QSA. We also separately consider the cases of the symmetron coupled only to DM and to all mater. We find that in both cases the QSA provides an accurate solution when the combination of the Compton wavelength and the redshift of the symmetry breaking transition is relatively small, but breaks down when this combination is large. Generally, the QSA holds when the background solution for the scalar field evolves sufficiently slowly over a period of oscillation of the scalar field perturbation, set by the inverse of the Compton wavelength. This paves the way to using MGCAMB and MGCLASS for deriving observational constraints on the parameters of the symmetron model and other scalar-tensor theories.

II. THE QUASI-STATIC APPROXIMATION FOR A SCALAR FIELD COUPLED TO DARK MATTER AND ALL MATTER
In this Section, we present the equations of motion of a cosmological theory that includes a scalar field conformally coupled to matter, starting with the case where the scalar field is coupled only to DM. Throughout this paper, g µν denotes the baryon frame metric whose geodesics are followed by the Standard Model (SM) particles, or "baryons", since the observable quantities are measured in this frame. We use the metric signature (−, +, +, +).
Evolving the scalar field equations can be time consuming due to rapid oscillations, making it impractical for applications requiring many repeated runs, such as MCMC parameter constraints. Instead, in some circumstances, one can use the QSA in which the scalar field perturbation is algebraically related to the matter density contrast. The QSA involves two assumptions: (1) that the Fourier modes of interest are sub-horizon, k/H ≫ 1, and (2) that the time derivatives of the metric and the scalar field are much smaller than their spatial derivatives. In ΛCDM, the sub-horizon condition automatically implies the smallness of the time derivatives of the linear metric perturbations, as their growth rate becomes comparable to spatial derivatives only for perturbations on near-Hubble scales. In scalar-tensor theories, however, the scalar field has its own dynamics, and the first assumption no longer implies the second. Moreover, the time-derivatives of the scalar field and the metric are not necessarily small. However, the QSA can still be a good approximation if neglecting the time-derivatives has no impact on observables. For instance, rapid time-oscillations of the scalar field around a central position may simply average out in their contribution to the observed matter power spectrum [25].
In what follows, we present the exact and the QSA form of the relevant perturbation equations with more details provided in the Appendix.
A. Scalar field coupled to dark matter Let us consider the following action: with the scalar field ϕ and the potential V measured in M Pl and M 2 Pl Mpc −2 units, respectively, where M Pl ≡ (8πG) −1/2 is the reduced Plank mass (Mpc denotes megaparsecs). In the above, L SM is the Lagrangian density of the standard model of particle physics, describing the known baryons, leptons and gauge bosons, and L dm is the Lagrangian of DM particles. The metric g µν determines the curvature scalar R and defines the baryon frame in which the SM particles follow the geodesics of g µν . The DM particles follow the geodesics ofg µν , related to g µν through a conformal transformationg where A(ϕ) is a coupling function. When the coupling function is non-constant, A, ϕ ̸ = 0, the scalar field mediates a fifth force between DM particles. In this DM-only coupled case, the form of the Einstein-Hilbert action is unmodified in the baryon frame. Next, we consider the perturbed flat Friedmann-Lemaitre-Robertson-Walker (FLRW) metric in the Conformal Newtonian gauge, where a is the scale factor and τ is the conformal time. The two scalar potentials, Ψ and Φ, characterize the Newtonian gravitational potential and the curvature perturbations, respectively. We write the scalar field as a sum of the homogeneous (background) part and a perturbation which depends on space and time, and, similarly, split the stress-energy components into their background and perturbation parts. The complete set of equations for the background quantities and linear metric and matter perturbations, before and after applying the QSA, is given by where H =ȧ/a, with the dot denoting the derivative with respect to τ , ρ dm , ρ b , ρ γ , and ρ ν are the DM, baryon, photon, and neutrino energy densities, respectively, while δ i , θ i , and σ i are the energy density fluctuations, divergence of fluid velocity, and anisotropic stress, respectively, following the notation of [26]. The coupling function β = A, ϕ /A determines the strength of the DM-scalar interaction, c s is the baryon sound speed, andτ c is the differential optical depth which is important before and around the epoch of recombination, where the photons and baryons are tightly coupled. The effective potential V eff is defined via its derivative appearing on the right hand side of Eq. (5.2), and the effective scalar field mass appearing in Eq. (6.10) is m 2 = V ′′ eff (ϕ min ), where ϕ min is the scalar field at the minimum of V eff . The continuity equation for DM, Eq. Under the QSA, where the Fourier modes of interest are sub-horizon, the time derivatives of both the metric and the scalar field perturbations are assumed to be negligible compared to their spatial derivatives and removed from Eqs. (5.4), (5.6), (5.8), and (5.10). At the background level, the QSA implies that the scalar field evolves slowly, remaining at the minimum of the slowly changing effective potential with no oscillations around the minimum. In this approach, one eliminates the derivative of the background scalar field in equations where it appears, and obtains an algebraic relation between the background scalar field and the matter density. This is justified as long asφ ≪ H. One should, however, keep in mind that this is not a good approximation in models in which βφ ∼ H [27]. If the background field is at the minimum of the effective potential, we have V ′ eff = 0, and the DM density perturbation is then also algebraically related to the scalar field perturbation δϕ via Eq. (6.10). Also, the contribution of the scalar field density perturbations, V, ϕ δϕ, to the Poisson equation is negligible compared to the matter density, hence we have omitted it in Eq. (6.4). The complete set of equations that we numerically evolve to compare the exact solution to those in the QSA are given in Appendix A 1, along with additional details on how they are derived.

B. Scalar field coupled to all matter
The action in the case of a scalar field conformally coupled to all matter, when written in the Einstein frame (in which the gravitational part of action has its usual form) is similar to the one in Eq. (1), except now both normal matter and DM follow the geodesics ofg µν . In the Einstein frame, the baryon continuity and Euler equations acquire the explicit ϕ-dependent terms previously present only in the DM equations. Hence, the baryon and DM perturbations have similar equations of motion, with the additional coupling to photons for the baryons. In our numerical implementation, we evolve the equations in the Einstein frame and use the conformal transformation to convert to the baryon frame g µν when interpreting the results, since it is the frame in which the observational measurements are made. The action in the baryon frame takes the form and the baryon frame stress-energy is The baryon frame QSA equations that are modified compared to their counterparts in the DM-only coupled case of Eqs. (6) are . The full set of equations and the conversion between baryon and Einstein frames for all variables are provided in Appendix A 2.
The last term on the right-hand side of (10.2) arises from the conformal transformation of the metric from the Einstein to Jordan frame and describes the effect of the fifth force mediated by the scalar field. In the Einstein frame, the equivalent contribution of the fifth force manifests as an additional force term on the right-hand side of the Euler equation (6.9). It is important to note that both (6.4) and (10.2) are intended to encompass the scalar field density perturbation on the right-hand side in the sum over all species. However, their effect is negligible on sub-horizon scales due to the fact that the sound speed of the propagating degree of freedom is equal to the speed of light (c s = c), hence the scalar field density contribution is dropped under the QSA.
There are two notable differences between Eqs. (6) and (10). One is that δϕ, which is the fifth force potential, is now sourced by both DM and baryon inhomogeneities, implying a larger enhancement of the growth rate. The other important distinction is that the two gravitational potentials, Φ and Ψ, are different, developing the so-called "gravitational slip" [28]. As a result, the effective gravitational couplings affecting relativistic particles, which follow the geodesics of Φ + Ψ, and non-relativistic matter, which follows Ψ, become different, making it possible, in principle, to distinguish between the DM-only and all-matter coupling. However, confirming this in practice is challenging, as it is difficult to find observables that measure the true Ψ [19].
A commonly used practical way to model cosmological perturbations in modified gravity is to introduce the phenomenological functions µ, γ, and Σ, which parameterize the deviations of Eqs. (10.2) and (10.3) from their ΛCDM forms. These functions are identified by where ρ∆ ≡ ρ i ∆ i , with ∆ i = δ i +3Hv i /k being the comoving density contrast, and v i is the irrotational component of peculiar velocity. Under the QSA, ∆ i = δ i . These phenomenological functions are also known as G matter = µ(a, k)G and G light = Σ(a, k)G. Note that the three functions are related via Σ = µ(γ + 1)/2. Combining Eqs. (10.2), (10.3) and (10.5), and considering the fact that V ′ eff = 0 in QSA implies V, ϕ = −β(ρ b + ρ dm ), one finds the expressions for the phenomenological functions to be (see also [29,30])

III. THE SYMMETRON EXAMPLE
In what follows, we introduce the symmetron model in the CDM-coupled case, with all main features being the same for the all-matter coupled scalar field. In this model, the potential and the coupling function are given by [13] where V 0 is set by the present value of the DE density, and (µ, M, λ) are three parameters whose values are discussed later in this section. Using the DM energy density given by Eq. (6.3) in Eq. (7) and integrating allows us to write the effective potential as where A 0 = 1 + ϕ 2 0 /2M 2 is the present value of the coupling function. Putting together Eqs. (13), (14) and (15) allows us to write the effective potential in the symmetron model as The effective potential changes its shape depending on the DM background density. At low matter densities, the effective potential has two minima with ϕ ̸ = 0, hence A ̸ = 1 and DE and DM are non-minimally coupled to each other. However, at high densities corresponding to the early universe or highly clustered regions such as the Solar System, V eff has a minimum at ϕ = 0, implying A(ϕ) = 1 and no direct coupling between DE and DM. This is an example of screening, allowing to suppress the fifth force inside the Solar System and stay in compliance with the existing tests of GR [24]. Under the QSA, the scalar field always remains at the minimum of the effective potential. This condition leads to a useful expression for ϕ at densities below the point of spontaneous symmetry breaking (SSB), where is the scale factor at the SSB. The three parameters (µ, λ, M ) appearing in the symmetron Lagrangian can be expressed in terms of three more phenomenologically transparent quantities: the redshift at symmetry breaking, z SSB , the present value of the coupling strength, β 0 , and the Compton wavelength, λ C , which sets the range of the DM-DE interaction. The relations are given by [31] (see Appendix C for more details) The effective mass parameter appearing in Eq. (6.10) can be written as [32]

IV. TESTING VALIDITY OF THE QSA
To test the validity of the QSA, we compare the numerical solutions of the exact equations describing the evolution of linear perturbations in the symmetron model to their solutions obtained under the QSA, considering the cases of the scalar field coupled to DM and to all matter. We are especially interested in the evolution of perturbations after the point of spontaneous symmetry breaking, when the symmetron scalar field develops a vacuum expectation value and starts mediating a fifth force between matter particles.
We fixed the value of the coupling parameter at β 0 = 1 and used the set of equations provided in Appendices A 1 and A 2 to solve for the evolution of cosmological perturbations for a wide range of the other two symmetron model parameters: 1 ≤ z SSB ≤ 500 and 0.1Mpc ≤ λ C ≤ 100Mpc. For the all-matter coupled case, we evolved the equations in the Einstein frame, given by (A31), and transformed the solutions to the baryon frame using (A28). We set the initial conditions deep in the radiation era, as detailed in Appendix A 3. When working under the QSA, we evolve ΛCDM equations up to the time of SSB, and switch to Eqs. (6) and (B4) after the that.   1, and 10. We plot the matter density contrast, δ m = (ρ b δ b + ρ dm δ dm )/(ρ b + ρ dm ), the lensing potential Φ + = (Φ + Ψ)/2, and the product δ mΦ+ that represents the strength of the galaxy-CMB correlation induced by the Integrated Sachs-Wolfe (ISW) effect. All the results are shown relative to the corresponding solution in the ΛCDM model, and we consider both the DM-only and the all-matter coupled cases. In the left block of panels, where the range of interaction between the matter and DE is smaller (λ C = 1Mpc), there is a good agreement between the exact and the QSA solutions for all the values of k-modes. We also see that, as expected, the modifications relative to ΛCDM are stronger in the case of the all-matter coupling. On the other hand, in the right block of panels, corresponding to λ C = 50Mpc, we clearly see that the QSA approximation breaks down.
To separate the effects of assuming the QSA for the evolution of the background from the QSA for the perturbations, we also evolved the exact equations for linear perturbations on the QSA background, and the QSA perturbations on the exact background. This test revealed that, when using the QSA for the background, it makes no difference whether one uses the exact or the QSA perturbations. On the other hand, when using the exact background, the QSA breaks down in the same way whether one uses the exact or the QSA equations for the perturbations. This is illustrated in Fig. 1 which includes the results for the combination of the QSA background and exact perturbations, along with the fully exact and fully QSA solutions.
To assess the accuracy of the QSA, we calculated the root-mean-square (rms) of the fractional difference in the matter density contrast δ m between the exact solution and the QSA for two symmetron models depicted in Fig. 1. The corresponding results are presented in Fig. 2. The top two plots demonstrate the effectiveness of the QSA for λ C = 1Mpc, with differences of less than 1% observed for the DM-only coupled case and less than 1.5% for the all-matter coupled case when comparing the exact solution with the QSA. Conversely, the bottom two plots, corresponding to λ C = 50Mpc, reveal a significant discrepancy. The discrepancy reaches 20% for the CDM-coupled case and extends to 80% for the all-matter coupled case, indicating a failure of the QSA in these scenarios. Examining the solutions for the two models shown in Fig. 1, as well as our numerical results for other values of symmetron parameters, we find that the agreement between the exact and the QSA cases depends primarily on the dynamics of the background value of the scalar field after the point of symmetry breaking. Fig. 3 clearly shows the difference between the two cases. A larger Compton wavelength parameter, λ C ∼ m −1 , causes a longer delay in response of ϕ to the change in the shape of the effective potential at the SSB. It also reduces the frequency of oscillations of the background field around the average value predicted by the QSA solution, and allows for more freedom, i.e. a larger amplitude, for the scalar field to oscillate around that average due to a flatter shape of the effective potential at the minimum, since V ′′ eff = m 2 . We note that Eq. (5.9) suggests that QSA may not be a good approximation when βφ ∼ H. However, our analysis shows that, for the symmetron model, the Hubble friction is always dominant in this equation even for the large values of λ C . Our numerical results show that the QSA holds when 1 ≤ z SSB ≤ 10 and λ C < 10 Mpc, with the best consistency for λ C = 1 Mpc. This range of validity of the QSA also follows from a simple analytical consideration presented below.
Prior to specializing to the symmetron example, it helps to note that the QSA assumes that the background scalar field follows the minimum of the effective potential V eff (ϕ), which includes a term proportional to the decreasing background matter density. As a result, the minimum of V eff is continuously shifting and the QSA assumes that the scalar field manages to keep up with the shift. In addition, the scalar field oscillates around the minimum, which is not described by the QSA. However, in most cases of cosmological relevance, these oscillations do not affect the validity of the QSA when it comes to evaluating observable quantities, as the oscillations tend to average out to zero over cosmological timescales. When the scalar field evolution is dominated by the Hubble friction, it is unable to keep up with the shift in the minimum, and its average value becomes inconsistent with the prediction of the QSA. Hence, the criterion for the validity of the QSA is set by the balance between the time-scale m −1 , set by the curvature of the effective potential V ′′ eff , and the Hubble scale H −1 . These considerations are generally applicable to all scalar-tensor theories of generalized Brans-Dicke type [24,33].
Let us now consider V eff (ϕ) of the symmetron model near ϕ = 0 right after the point of symmetry breaking. Since the value of the scalar field was close to zero at that point, we can approximately write V eff (ϕ) as , and we have omitted the λϕ 4 term since it is small compared to the ϕ 2 terms, and also ignored the matter density term, ρ (0) dm /(A 0 M 2 a 3 ), as it decreases rapidly with a and quickly becomes much smaller than µ 2 . The equation of motion for the scalar field right after the SSB can be written as We consider two regimes for the transition from a zero expectation value of ϕ to a none-zero value at the SSB: a fast transition, whereφ ≫ Hφ and the Hubble friction can be neglected, and a slow transition dominated by the Hubble friction, withφ ≪ Hφ. We can estimate the time-scale of the transition in the slow and the fast cases, T s and T f , as We expect the QSA to break down when the slow and fast transition time-scales are of the same order, i.e. T s ∼ T f or 2H SSB ∼ µa SSB . From this condition and Eq. (19a), which relates µ to λ C and a SSB , we can derive an expression for the critical value of the Compton length above which the QSA fails: where H 0 is the Hubble constant, and we used H SSB ∼ H 0 √ 1 + z SSB assuming the symmetry breaking occurs during matter domination. The plot of λ cr vs z SSB in Fig. 4 shows that the analytically determined range of validity of the QSA, λ C < λ cr , is in good agreement with our numerically obtained results, such as the two examples shown in Fig. 1.

V. SUMMARY
Scalar-tensor theories of gravity provide a well motivated framework to study the universe beyond the ΛCDM model. Among such theories, the symmetron model is a compelling example with a screening mechanism and a renormalizable potential, making it an excellent prototype for late dark energy.
Although many studies of observational constraints on scalar-tensor theories assumed the QSA to be valid in the late universe, there is no guarantee that it holds for all ranges of the parameters in specific modified gravity model. In this work, we focused on the symmetron model to examine in detail the range of validity of the QSA. We found that it is the evolution of the background scalar field that is most important. Namely, whether its dynamics is determined by the Hubble friction or the scalar field potential. In the friction dominated case, m ≪ H, the field is unable to keep up with the shifting minimum of the effective potential, and the QSA is broken. In the opposite limit, the QSA works. We demonstrated this numerically by comparing the exact solutions for the background and perturbations in the symmetron model with the QSA prediction. We also derived an expression (24) for the critical Compton wavelength, given the redshift of SSB, above which the QSA breaks down. Our arguments can be applied to derive similar criteria in other scalar-tensor theories of generalized Brans-Dicke type.
Our analysis paves the way for using MGCAMB and MGCLASS, that assume the QSA, for deriving observational constraints on the parameters of the symmetron model and other scalar-tensor theories.

Appendix A: Background and perturbations in scalar-tensor theories
Here we present the complete set of equations that we use to evolve the background and perturbations for the exact case.

CDM-only coupled case
Varying Eq. (1) with respect to g µν yields the Einstein's field equations, where T SM µν andT dm µν are the stress-energy tensors of normal matter and DM in baryon and DM frames, respectively, given by with Eq. (A2b) relating the DM stress-energy tensors in two frames. The scalar field stress-energy tensor is Varying the action with respect to ϕ yields the scalar field equation of motion (EoM), whereT dm is the trace ofT dm µν that transforms between two frames through We will assume that, in the DM frame, the DM acts as dust, i.e. a perfect fluid with no pressure. Since we are writing our equations in the baryon frame, which is the co-moving frame for normal matter, the stress-energy tensor of baryons, photons, and neutrinos is conserved, The Bianchi Identity for the Einstein tensor G µν implies that the sum of the scalar field and DM stress-energies should be conserved, namely, which, after using Eqs. (A3) and (A4), gives The effect of the fifth force on DM particles enters through the right hand side of Eq. (A8) and vanishes when A, ϕ = 0. The background expansion is given by the Friedmann equation, The background evolution of the scalar field is given bÿ and the DM density evolves according toρ For numerical implementation in codes such as CAMB [34], it helps to write the above equations as where Ω (0) i is the density parameter for an individual component at the present epoch [1], and Ω dm ≡ ρ dm /(3M 2 Pl H 2 0 ). Here, ϕ has units of M pl and the conformal time is in units of megaparsec [Mpc]. Eqs. (A12) form a complete set needed to solve for the background evolution.
Eqs. (A1), (3), and (4) imply four equations for the linear metric perturbations in Fourier space, but one only needs two out of four Einstein equations given by with the choice being a matter of convenience. We select (A13a) and the derivative of (A13b), since bothΨ andΦ appear in Eq. (A14). Using Eqs. (3) and (4) in (A4), yields the EoM for scalar field perturbations: Perturbing Eq. (A8) gives the continuity and Euler equations for the DM perturbations: For the baryons, the conservation equations are the same as in ΛCDM, since they do not couple to the scalar field: The quantity c s is the baryon sound speed given by where R = 3ρ b /4ρ γ is the baryon-photon density ratio.τ c = an e σ T is the differential optical depth, n e is the free electron density and σ T is the Thomson scattering cross section. The Thomson opacity is important before and around the epoch of recombination, when the photons and baryons are tightly coupled [26]. In [35], Hu and Sugiyama provided a fitting formula forτ cτ which can be used for quick estimates once the ionization history is provided. In the above, Y p ≈ 0.23 is the primordial helium mass fraction and h is the dimensionless Hubble parameter. The photon and neutrino phase-space distributions evolve according to the perturbed Boltzmann equations and can be expressed in terms of perturbations in their corresponding black-body temperatures. Expanding the temperature anisotropies in Legendre multiple moments gives an infinite set of coupled differential equations known as the Boltzmann hierarchy [36]. The hierarchy can be truncated at some maximum multiple moment l max chosen according to the desired numerical accuracy. For our purposes, it will be sufficient to work with l max = 3.
A full list for evolution of perturbations is given bẏ where following the truncation scheme provided by Ma & Bertschinger in [26], the equations for the photon and neutrino anisotropies can be written as (A19.10) to (A19.13). Here, F γl and F νl denote the photon and neutrino multiple moments of order l, with F 0 = δ, F 1 = 4θ/3k, and F 2 = 2σ. Note that the evolution of DM is modified via coupling to the scalar field as (A19.4), (A19.5) while baryons are only coupled to photons in (A19.6), (A19.7).

The all-matter coupled case
For the all-matter coupled case, we evolve the equations in the Einstein frame, and transform the solutions to the baryon frame. A conformal transformation transforms the Einstein frame metricg µν to baryon frame metric g µν via Hence, the line element transforms as Since d⃗ x is a comoving coordinate, it is not subject to change of scale in conformal transformation. This immediately gives two expressions relating the scale factor and conformal time between the baryon and Einstein frames as i.e. the conformal coordinate τ is the same in the two frames. The transformation between the cosmological variables in the two frames can be derived from requiring the action to remain invariant, i.e. S =S. Furthermore, we require the scalar field to have the same form of the kinetic and potential energy terms, giving the transformation for the scalar field and its potential [37]: Using Eq. (9) and the invariance of the action allows one to derive the transformation for the matter stress-energy: giving the transformations for the background pressure and energy density: Since the stress-energy tensor of relativistic particles is traceless, the conservation equation for radiation density has the same form in the Einstein and baryon frames,ρ One can use this to derive a useful relation between the conformal Hubble parameter in the two frames: The transformations for the linear perturbations can be found by perturbing Eqs. (A23a) and (A25), and from Eq. (A21). The complete set of relations between the Einstein and baryon frames is given bỹ The divergence of the fluid velocity θ and the shear stress σ in k-space are defined as [26] (ρ + P )θ ≡ ik j δT 0 j , It then follows from Eqs. (A24) and (A29) that these two variables remain invariant under conformal transformations. Eqs. (A23), (13) directly relate the parameters in symmetron model between two frames by The equations of motion when scalar field is coupled to all matter in the Einstein frame are given bỹ is the derivative of effective potential with respect toφ. The conversion of all the variables to the baryon frame is given by (A28), (A30).

Initial conditions
The system of coupled differential equations presented in earlier sections can be soled once the initial conditions are specified. We set them deep in the radiation era where the relevant Fourier modes are still outside the horizon. This satisfies the condition kτ ≪ 1. At the early time, the baryons and CDM have a negligible contribution to the energy density of the universe, so we have ρ total ≈ ρ ν + ρ γ . The Hubble parameter is also inversely proportional to the conformal time and determined by H = τ −1 . We choose our initial conditions so that only the fastest-growing physical modes are present. This is appropriate for the perturbations that are created in early universe. We can then analytically find the time dependence of all metric and density perturbations by directly solving the differential equations similar to the ones in [26]. Using Eqs. (A9), (A19) for the isentropic perturbations in the conformal Newtonian gauge, we derive the following relations for the initial conditions: where β = A, ϕ /A, and R ν = ρ ν /(ρ ν + ρ γ ) is the neutrino-photon density ratio and C is a normalization constant which sets the initial value of Ψ (we choose C = 0.5 in accordance with its value in CAMB [34]). One can check that with the absence of scalar field, these equations will reduce to the corresponding initial conditions in ΛCDM model. Eqs. (A33) apply to a general modified gravity model with the scalar field non-minimally coupled to DM. What is missing is the choice for the initial value of the scalar field and its derivative, which need to be chosen based on the properties of a given model. Since we are interested in the symmetron model, in which the VEV of the scalar field is zero prior to symmetry breaking, we set ϕ to be a small nonzero number, ϕ ini ≪ M pl , andφ = δϕ = δφ = 0 at the initial time.

Appendix B: The quasi-static approximation
In what follows, we present the QSA form of the relevant perturbation equations.

QSA for the scalar field coupled to dark matter
We revisit the equations of the previous section and apply the QSA. We take k ≫ H and assume that Under these approximations, the perturbation of the scalar field δϕ is algebraically related to the DM matter density perturbation, and so are the gravitational potentials Φ and Ψ: where V ′ eff is given by Eq. (7). Since QSA implies that the scalar field remains at the minimum of the effective potential, we have V ′ eff = 0 . Also, the contribution of the scalar field density perturbations, V, ϕ δϕ, to the Poisson equation is negligible compared to the matter density, hence we will omit it in Eq. (B2b). With this, we can rewrite Eqs. (B2) as where with times and distances are now in Mpc units and m 2 = V ′′ eff (ϕ min ), with ϕ min being the minimum of V eff .

QSA for the scalar field coupled to all matter
Applying the QSA in the case of the scalar field coupled to all matter, the affected equations in the Einstein frame areH ϕ =φ min , (B4.2) Transforming these equations into the baryon frame gives Eqs. (10), where we have omitted theṼ ,φ δφ term in (B4.5) due to its negligible effect.