Unveiling dark fifth forces with linear cosmology

We initiate the exploration of the cosmology of dark fifth forces: new forces acting solely on Dark Matter. We focus on long range interactions which lead to an effective violation of the Equivalence Principle on cosmological scales today. At the microscopic level, the dark fifth force can be realized by a light scalar with mass smaller than the Hubble constant today ($\lesssim 10^{-33}\,\text{eV}$) coupled to Dark Matter. We study the behavior of the background cosmology and linear perturbations in such a Universe. At the background level, the new force modifies the evolution of the Dark Matter energy density and thus the Hubble flow. At linear order, it modifies the growth of matter perturbations and generates relative density and velocity perturbations between Dark Matter and baryons that grow over time. We derive constraints from current CMB and BAO data, bounding the strength of the dark fifth force to be less than a percent of gravity. These are the strongest constraints to date. We present potential implications of this scenario for the Hubble tension and discuss how our results are modified if the light scalar mediator accounts for the observed density of the Dark Energy. Finally, we comment on the interplay between our constraints and searches for violations of the Equivalence Principle in the visible sector.


Introduction
The existence of Dark Matter (DM) is supported by a large variety of observations. This strong evidence, however, concerns only its gravitational interactions, providing no insights on its microscopic properties. Current data leave wide open the landscape of possible DM masses and interactions, and besides the development of a large number of experimental strategies, we are just scratching the surface of the theoretically plausible microscopic theories of DM.
Facing this plenitude of theoretical possibilities, it is interesting to take a more datadriven methodology and focus on what we can learn about the nature of DM using present and forthcoming observations. Following this approach, we initiate the exploration of DM scenarios that can be uniquely tested with cosmological measurements at an exquisite level of precision. These are dark sectors where DM is self interacting through a dark fifth force which has a range comparable to the size of the observable Universe.
Cosmology offers the unique possibility of studying the very long range dynamics of DM, where the presence of any force other than gravity results in an effective violation of the weak Equivalence Principle (EP). The idea is to use the Universe as a "scale" and try to measure by how much the inertial and gravitational masses of DM can differ in light of the available data. In this sense, cosmological observations can be viewed as dark fifth force experiments, in analogy with the many tests of long range forces in the visible sector [1][2][3].
In this paper we begin a systematic study of the cosmological signatures of dark fifth forces. As the first step we focus on the constraints from linear cosmology, looking at the current Cosmic Microwave Background (CMB) data from Planck [4] and Baryon Acoustic Oscillations (BAO) data [5][6][7][8]. For concreteness, we consider the case where the dark fifth force is mediated by trilinear interactions of the DM with a light scalar mediator, whose mass we take to be lighter than or equal to the size of the Universe today (m φ ≲ H 0 ≃ 10 −33 eV). We leave the exploration of constraints from non-linear cosmology, as well as different mass regimes and non-minimal realizations of a dark fifth force, for future works.
The presence of a long range force between DM particles severely affects both the cosmological background and the linear regime of perturbations, resulting in a strong upper bound on its strength. We will show this is constrained to be roughly two orders of magnitude weaker than gravity, as long as the mass of the scalar and its initial displacement from the origin are small enough to make its contribution to the vacuum energy negligible. This scenario, which we call a "pure fifth force" or 5F for short, was first studied in Refs. [9][10][11][12] (see also later discussions of possible observational imprints [13][14][15][16][17], including on galactic scales [18,19]).
In the opposite regime, the light scalar vacuum energy can account for the observed amount of Dark Energy (DE) today. In this scenario, which we call Coupled Dark Energy (CDE), the scalar is a quintessence field coupled to the DM through trilinear interactions and minimally coupled to gravity. Such interaction between DM and DE was first studied in Refs. [20][21][22] and has been subject of intense activity in subsequent literature [23][24][25][26][27][28][29][30][31][32][33][34][35][36]. In our simple setup and with natural values for the scalar potential, the tuning of the initial displacement of the scalar field required to match the value of the Cosmological Constant (CC) does not allow us to take the mass of the scalar arbitrarily small and at the same time retain an observable coupling with DM. Saturating this condition (i.e. taking the highest possible mass m φ = H 0 ) yields constraints on the strength of the coupling similar to the ones for a pure fifth force.
Finally, we take a look at the generic implications of the presence of a dark fifth force on a consistent relativistic theory of DM. Requiring naturalness of the scalar potential forces the DM mass to be below about 10 −2 eV. This bound could be overcome in nonminimal models (for instance imposing supersymmetry in the dark sector), but it indicates a generic correlation between a dark sector endowed with a fifth force and ultralight DM produced in the early Universe through the misalignment mechanism [37][38][39]. We explore the parameter space of simple axion DM scenarios and discuss the interplay between the dark fifth force constraints derived here and fifth force constraints in the visible sector.
This paper is organized as follows. In Sec. 2 we illustrate our framework, deriving the equations that govern the dynamics of DM in the presence of a fifth force and of the fifth force field itself. We then study the behavior of the background in Sec. 3 and move to the linear perturbations in Sec. 4. The constraints derived from CMB and BAO data are discussed in Sec. 5. In Sec. 6 we comment on the interplay between our constraints and searches for EP violation in the visible sector. Section 7 contains our conclusions. In a triptych of appendices (A, B, C) we give further technical details about the equations governing the cosmological dynamics.

Framework
In this section we show how to describe the cosmological imprints of a dark fifth force. We start in Sec. 2.1 by studying the behavior of a DM particle in the presence of an external gravitational field plus a scalar fifth force. We then derive the equations governing a fluid of DM particles in the non-relativistic limit. This limit is the one relevant for DM, although in principle our derivation could apply to relativistic fluids like neutrinos [40,41]. The relativistic expressions are provided in Appendix A. In Sec. 2.2 we give explicit examples of full relativistic quantum field theories realizing this fluid behavior and outline their properties.
We assume the Universe is described by a flat FLRW metric and consider only scalar perturbations around it, which can be written in Newtonian gauge as where Ψ and Φ coincide with the two gauge invariant Bardeen potentials [42]. For some applications we will make extensive use of other gauges, in particular the synchronous gauge, to which coordinate transformations are well known [43]. 1 The equations governing the perturbations in synchronous gauge are provided in Appendix B. 1 The synchronous gauge is a full gauge fixing of the metric only if the perturbations are defined in a frame comoving with one (or more) pressureless species that follows the standard General Relativity geodesics. This is not the case for DM interacting with new light degrees of freedom, hence strictly speaking the synchronous gauge cannot be defined in our setup. We then defined the synchronous gauge with respect to a small fraction of regular cold Dark Matter (CDM) ρc, not subject to the new long range force. In

Particles and fluids
In the classical limit, we can derive the worldline of a DM particle (independently of its spin) in an external background with a non trivial metric g µν and a scalar fifth force s [22] as for some affine parameter λ. We take m χ (s) to be a function of s only, whose explicit form depends on the microscopic interactions between the DM and the scalar fifth force. By varying the action we find the geodesic equation where dσ = −g µν dx µ dλ dx ν dλ dλ and the Γ's are the usual Christoffel symbols. The fact that the DM mass does not drop out from the integral in Eq. (2.2) indicates a violation of the EP in the dark sector, which results in a modification of the DM geodesic equation controlled by the space-time profile of the scalar field s.
We define the 4-momentum of the χ particles as with the mass shell condition g µν P µ P ν = −m 2 χ (s) providing the dispersion relation E = m 2 χ (s) + p 2 . This allows us to rewrite the geodesic equation as dP µ dλ + Γ µ νρ P ν P ρ + 1 2 ∂ m 2 χ (s) ∂s ∂s ∂x ν g µν = 0 , (2.5) which in the non-relativistic limit and in flat space-time reduces to the formẍ i = −∇ i Ψ − (∂ log m χ (s)/∂s)∇ i s. We clearly see that the particle acceleration feels both the gravitational force (the gradient of the gravitational potential) and the scalar fifth force. The DM single particle phase space density f χ (x, t, p) satisfies the Vlasov equation [44] ∂f χ ∂t where we treat the DM as a collisionless fluid and encode its interaction with the fifth force in the geodesic equation, very much as we do with gravity. This treatment is consistent as long as the fifth force mediator s has a large occupation number and can be viewed as a coherent field acting on test particles. We then proceed and compute the first two moments of the Vlasov equation, assuming χ to be pressureless (i.e. neglecting terms of order p 2 /E 2 or higher) and expanding to the particular, we set Ωc = 10 −5 today and we checked that this value has no observational impact given the current uncertainties of cosmological data. linear order in the metric perturbations. The resulting continuity and Euler equations reaḋ where ρ χ = d 3 pEf χ /(2π) 3 is the energy density, v i χ = ρ −1 χ d 3 pp i f χ /(2π) 3 is the fluid velocity, H is the Hubble parameter and () · indicates derivatives with respect to the coordinate time t. The equations are fully non-linear in the χ density and velocity perturbations. The piece of Eq. (2.8) containing the second moment of the phase space distribution, Σ ij χ = d 3 pp i p j f χ /[(2π) 3 E], needs to be included for a perturbative treatment beyond linear theory [45,46].
The presence of the new long range force modifies the Euler equation in Eq. (2.8), introducing a new acceleration term similar to the one induced by the gravitational potential, but controlled by the field dependence of the DM mass. This was expected, since the classical non-relativistic 1/r potential generated by the exchange of the fifth force adds to the gravitational potential in flat space-time. Interestingly, the effect of the new force can be rewritten as the gradient of the inertial mass, ρ χ (∂ log m χ (s)/∂s)∇ i s = n χ ∇ i m χ (s), which shows once again how the presence of a scalar force induces a violation of the EP.
The continuity equation in Eq. (2.7) is also modified, reflecting the fact that only the total energy density of χ and s is conserved. This introduces an explicit dependence of the DM density on the time evolution of the fifth force field s , which will impact the behavior of the cosmological background (discussed in Sec. 3) and the growth of density and velocity perturbations (discussed in Sec. 4). These effects cannot be easily decoupled parametrically from the modification of the Euler equation and need to be included in a complete description of the fifth force dynamics. As a result, the imprints of a fifth force in the cosmological observables cannot be fully captured by looking purely at the non-relativistic limit of the Vlasov equation, as done for example in Refs. [12,47].
A completely equivalent way of deriving the fluid equations above would have been to start from the energy-momentum tensor for the microscopic degrees of freedom, and then rewrite it in terms of the macroscopic fluid variables. We explicitly show this derivation in Appendix A.

Fields
The continuity and Euler equations for DM should be complemented with the equations describing the dynamical evolution of the scalar fifth force and gravity. In general, the microscopic Lagrangian controlling the fifth force and DM dynamics can be parametrized as where K χ is a canonical kinetic term for the DM and V χ is its potential energy. 2 On the other hand, K s and V s are the kinetic and potential energy of the fifth force field, and V int parametrizes the interactions between the two species. To strengthen the analogy with the gravitational potential we define the 5th force kinetic term as K s ≡ − ∂ µ s∂ µ s/(2G s ), where s is dimensionless, like metric perturbations are, and G s is a dimensionful constant analogous to G N . We thus define the dimensionless ratio which determines the strength of the 5th force with respect to gravity. The Einstein and Klein-Gordon (KG) equations are written as is the KG operator. Energy-momentum conservation for the sum of the DM and fifth force fluids requires ∇ µ (T µν χ + T µν s ) = 0 , but the individual stress tensors are not conserved. We thus have the freedom to choose which stress tensor the interaction term V int is contributing to. We find it convenient to package the interaction in the DM stress tensor, which makes transparent the appearance of the field dependent mass m χ (s) and its relation with the interaction term V int . With this choice we can where the DM energy density is given by ρ χ = 2(V χ + V int ) assuming the DM to be pressureless (see Appendix A for a detailed discussion).
At this point, we pause briefly to discuss the decoupling limit of the 5th force, which is not entirely trivial due to the use of the dimensionless field variable s. The rewriting K s = − ∂ µ s∂ µ s/(8πG N β) makes it apparent that in the decoupling limit β → 0, s ∼ β 1/2 is required in order to keep fixed the kinetic term of the 5th force field. Therefore the last piece on the left-hand side of the continuity and Euler equations, (2.7) and (2.8), vanishes. The decoupling is also manifest in the KG equation (2.12), where the right-hand side vanishes because V int is at least linear in s. To inspect the decoupling in the Einstein equation (2.11) it is useful to recall the structure of the DM stress tensor, whose last term vanishes for β → 0. In summary, in the decoupling limit we recover two non-interacting species (and fluids) with separately conserved stress tensors, as expected. So far our discussion has been general, but we now focus on the minimal scenarios of interest in this work, where both V χ and V int are quadratic in χ. We assume V χ to contain only a mass term for DM while the DM interaction with the fifth force takes the schematic form V int ∼ χ 2 F (s), with F (s) some function of the 5th force field. As discussed in Appendix A, in this case the relation ρ χ = 2(V χ + V int ) can be written as which provides an explicit map between the microscopic interactions in the field theory picture (the left hand side of the equation) and the fluid description (the right hand side of the equation). We now write two explicit particle physics models as examples of the microscopic interpretation to the parameters G s and m χ (s). The simplest way to realize a dark fifth force scenario is through a trilinear interaction between the Dark Matter particle χ and the scalar mediator φ (taken to have canonical field dimension). While in this paper we focus on this simple case, we emphasize that our formalism applies to any interaction of the form V int ∼ χ 2 F (s), provided F (s) can be expanded as a power series n ≥ 1 c n s n where c n are dimensionless coefficients.
Explicitly we can write where the DM χ is a Dirac fermion in Eq. (2.15) or a real scalar in Eq. (2.16). The non relativistic potential generated in Minkowski space-time by the exchange of the scalar fifth force with the DM particles as external sources can be written as V (r) = − g 2 D e −mφr /(4πr) independently of the DM spin. As a consequence, at distances r ≪ m −1 φ the long range force is equivalent to a shift in Newton's constant, G N → G N (1 + β), where β was given in Eq. (2.10) and we defined (2.17) The field dependent DM mass is m χ (s) = m χ (1 + s) for fermionic DM and m χ (s) = m χ √ 1 + 2s for real scalar DM. A non-polynomial dependence of the DM mass on the fifth force field can in principle be realized if φ is a non-minimally coupled scalar in Brans-Dicke theories [48], a string theory dilaton [49], or the dilaton of spontaneously broken conformal symmetry [50]. The cosmological implications of m χ (φ) = m χ e −βφ/M Pl for an ultralight scalar were considered already in the seminal papers on coupled quintessence [20,21]. In this work we do not consider these scenarios, leaving a detailed analysis for future study.

Fifth force potential and naturalness
As our goal is to study the phenomenology of a new long range attractive force operating in the dark sector, the force mediator φ is assumed to have physical mass smaller than or equal to the Hubble scale today and to interact with DM through trilinear interactions like the ones in Eq. (2.15) or Eq. (2.16). In general, we can consider the renormalizable potential for the scalar fifth force mediator where naturalness requires its parameters to be as large as the size of the radiative corrections controlled by the interaction strength with DM. Explicitly we find for both scalar and fermionic DM where M Pl = (8πG N ) −1/2 = 2.4 × 10 18 GeV is the reduced Planck mass and we used the definitions in Eq. (2.10) and Eq. (2.17) in the estimates above. For simplicity, in this paper we take V s to contain just a mass term m 2 φ φ 2 /2. Nonetheless, we have checked that a very similar cosmology is obtained by considering for instance a purely quartic potential with effective coupling λ φ = β m 2 φ /(s 2 M 2 Pl ), wheres is the typical field value during matter domination (see for example Fig. 1). Taking a linear combination of mass, cubic and quartic would also leave our results qualitatively unchanged.
We recall here that for a scalar field with physical mass much larger than H 0 , the purely quadratic and purely quartic potentials result in very different late-time evolutions, in particular for the cycle-averaged background energy density of the scalar, as thoroughly studied in the decoupled case [51][52][53]. An extension to include a coupling to DM has been made in Ref. [54], where a scalar field with purely quartic potential playing the role of Early DE was considered. 3 Requiring m φ ≲ H 0 and the naturalness of the fifth force potential in Eq. (2.19) can be read as an upper bound on the DM mass (see also Ref. [34]), where we set as reference for β the ballpark value of our cosmological constraints, as derived later in Sec. 5. Therefore, a relativistic field-theoretical description of the dark fifth force mediator without fine-tuning requires the DM to be an ultralight boson. In this mass region the DM may for instance be an axion, produced via misalignment or other non-thermal mechanisms. While this observation does not affect our analysis of the cosmology, we will return to its implications for particle physics in Sec. 6. One should also keep in mind that the naturalness bound in Eq. (2.19) can in principle be circumvented in non minimal DM models, for instance enforcing supersymmetry in the dark sector. In such non minimal scenarios the DM could be heavier and possibly fermionic. Having introduced our theoretical framework, we are now ready to present the background and first order cosmological dynamics in the next two sections. For definiteness, when providing numerical results we always refer to a representative model with scalar DM and purely quadratic fifth force potential, namely −L = (∂ µ χ∂ µ χ + m 2 χ (s)χ 2 )/2 + M 2 Pl (∂ µ s∂ µ s + m 2 φ s 2 )/β with m χ (s) = √ 1 + 2s . However, we emphasize again that our findings have little dependence on the DM spin and the precise form of the φ potential.

Cosmological Background
At the background level, the momentum of the χ particles redshifts with the expansion of the Universe like the one of all other particles,ṗ i = −Hp i . This guarantees the existence of a frame where the spatial part of the fluid 4-velocity of each species is zero, hence homogeneity and isotropy are preserved. We thus split each variable in a background component plus a perturbation around it, with the former depending only on time. For example ρ χ =ρ χ (t) + δρ χ (x, t) and s =s(t) + δs(x, t), and so on for the other species.
The equation of motion for the background DM density is directly obtained from the continuity equation (2.7). Switching to conformal time τ we havē where () ′ denotes derivatives with respect to τ and H = a ′ /a is the conformal Hubble parameter. The field dependent mass and its derivatives are, here and throughout the text, evaluated in the background. Eq. (3.1) can also be derived by writing in the nonrelativistic limitρ χ = m χ (s)n χ and imposing conservation of the comoving number density, n ′ χ + 3Hn χ = 0. We see that while the number density of DM scales proportionally to the volume (i.e.n χ ∼ a −3 ), the energy density does not, unlesss is constant. As we show below, the KG equation forcess ′ < 0 throughout the cosmological history so that the righthand side of Eq. (3.1) is negative and non-zero, corresponding to a net energy transfer from the DM to the 5th force field which modifies the redshift behavior of the DM energy density.
The KG equation for the scalar field is where V s,s ≡ ∂V s /∂s, and for the remainder of this paper by V s we always indicate the self interaction of s in the background, and similarly for its derivatives. The above equation is derived starting from the Lagrangian for the fields, Eq. (2.12), and rewriting the source term as a function of the χ fluid variables according to Eq. (2.14). At the background level the energy density and pressure of the 5th force mediator arē satisfying the equation of motion where w s ≡ P s /ρ s is the equation of state parameter of the 5th force fluid. The sum of this equation with (3.1) is identically zero on the right hand side, reflecting the energy conservation for the total fluid.

New cosmological parameters
The background evolution depends on three new parameters beyond the ones of ΛCDM: i) the dimensionless quantity β as defined in Eq. (2.10) determines the coupling strength of the 5th force in units of G N ; ii) the mass of the 5th force mediator m φ is assumed to be smaller than H 0 in order for the interaction to be long range till today; iii) the dimensionless initial background value of the 5th force fields ini controls the size of its contribution to the vacuum energy, together with m φ and β. 4 We consider two physically distinct scenarios where the scalar field s plays different roles in the cosmological history: • Pure Fifth Force (5F): the scalar field s leaves observable imprints only by mediating a long range force between DM particles, whilst its energy densityρ s remains negligible throughout the cosmological history. In this case we sets ini to a very small value (for definiteness we chooses ini = 10 −4 ) and determine the value of the CC density parameter ω Λ that satisfies the closure condition today, i ω 0 i = h 2 . In the 5F scenario the evolution of the s background undergoes the phases A and B described in Sec. 3.2.
• Coupled Dark Energy (CDE): in addition to mediating a 5th force, s also accounts for the DE. In this case we set the CC to zero, ω Λ = 0, and fix the initial field values ini by imposing the closure condition i ω 0 i = h 2 . The fine-tuning of the CC is thus replaced by the fine-tuning of the initial conditions for s. The required initial field value is always much larger than in the case of a pure 5th force: assuming a quadratic potential V s , for a scalar mass m φ ∼ H 0 one findss ini ∼ O(1) whereas for smaller masses the initial condition scales ass ini ∼ H 0 /m φ . 5 In the CDE scenario the s background evolution consists of phases A, B, and C as detailed in Sec. 3.2.
Finally, the energy density of the χ particles, ω χ , replaces the one of standard CDM. For practical reasons, Boltzmann codes like CLASS [55,56] take as input the present day value of the energy density of each species, ω 0 i , and evolve it back to the desired initial redshift using the well known behavior with the scale factor of CDM, baryon, photons, etc. This is not possible in our case, because we do not know a priori the exact time dependence of ω χ and ω s . For consistency with the other species, we use as input parameter the energy density that the χ particles would have had they evolved with a −3 like CDM. We dub this parameter ω d and set the initial density of χ as Clearly, ω 0 χ ̸ = ω d and the procedures described above enforce flatness. 4 The initial value of the derivatives ′ ini is also required to solve the KG equation, but this is determined as a function of the other input parameters as discussed in Sec. 3

Background behavior
The background fields exhibits a remarkably rich cosmological history, consisting of several successive epochs characterized by distinct dynamics. Its evolution directly impacts the one of the DM density, as can be seen by solving Eq. (3.1): Sinces is a monotonically decreasing function of time, 6 χ redshifts faster than CDM. We now present an analytical understanding of the solutions across cosmic time. The analytical approximations we provide are in very good agreement with numerical results, which are presented in Sec. 3.3 below.
A) Radiation dominated era:ρ s ∝ a −2 During radiation domination (RD, where H = 1/τ at zeroth order in β) the s selfinteraction can be neglected in the KG equation, V s,s ≪ρ χ (∂ log m χ (s)/∂s). We proceed perturbatively in β, making the approximations where where a solution scaling as τ −1 has been discarded, and we have used the leading-order expression a = (Ω 0 r ) 1/2 H 0 τ , with Ω 0 r = Ω 0 γ + Ω 0 ν if neutrino masses are neglected. Quantitatively, the linear dependence on a implies thats remains approximately constant until matter-radiation equality, as seen in Figs. 1 and 3. However, several insights can be extracted from Eq. (3.8). First, it determines the natural value ofs ′ ini as a function of the other input parameters. Second, via Eq. (3.3) it establishes that for the very small potentials V s relevant to this work, w s = +1 throughout RD. However, the energy density does not scale as a −6 , as it would for a decoupled field with w = +1, but redshifts asρ s ∝ a −2 instead. Consistency with the fluid equation (3.4) then requires the following quantity to be a constant of motion,ρ which can be easily checked using the approximate analytical solution and is verified numerically to high accuracy. The fraction of the total energy density stored in s grows rapidly asρ s /ρ tot ∝ a 2 until matter-radiation equality. This solution is in fact an attractor, to which trajectories with initial velocitys ′ ini larger or smaller than the one fixed by Eq. (3.8) rapidly converge. Such trajectories correspond to a nonzero coefficient for the τ −1 solution that was discarded in Eq. (3.8). If the coefficient is negative,s ′ ini is larger than the natural value andρ s undergoes an initial period of kination with a −6 scaling before converging to the attractor. If the coefficient is positive, there is an initial phase of kination that ends whenρ s reaches zero, after which the trajectory very quickly converges to the attractor.
where we defined with the last equality holding at zeroth order in β. A solution scaling as τ −3 has been discarded in order to match to the evolution in RD and the leading order expression of the scale factor a = ( Ω d + Ω 0 b )H 2 0 τ 2 /4 has been used. Since w s = +1, one immediately obtainsρ s ∝ a −3 and a simple estimate for the (constant) fraction of energy density in s (see also Ref. [21]), The same as in the left panels, but now assuming s is a Coupled Dark Energy field.
In addition, we derive perturbative solutions up to O(β) for Hubble and the DM density fraction, focusing on the 5F scenario wheres eq ≈s ini ≈ 0 in Eq. (3.10). The second Friedmann equation is written as and it is solved by where non-log-enhanced terms have been dropped, since they are numerically negligible. Thus, We also exploit Eq. (3.6) to write a perturbative expression for the DM density fraction, Figure 3. Same as Fig. 1, but for β = 0.02. The appearance of the logarithms in Eqs. (3.15) and (3.16) implies that the deviation from ΛCDM of the background evolution is parametrically enhanced by a factor ≲ log (1+z eq ) ≈ 8 with respect to the naively expected size of O(β). The log-enhanced effects encoded by these formulae match well the numerical results, as shown in Fig. 5. If s is a pure 5th force field, the above solutions persist until the CC comes to dominate the total energy density of the Universe. In this final phase of CC domination the dynamics is modified, but Eq. (3.12) still provides a good estimate of the s energy density fraction today, as can be seen in the left panels of Figs. 2 and 4.

C) Late matter dominated era:ρ s ≃ constant
If s is a coupled quintessence field the above discussion applies, but at some time during MD the equation of state eventually changes from w s = +1 to w s ≃ −1, crossing zero when the equality (s ′ ) 2 = 2G s a 2 V s is satisfied. For the assumed quadratic form of V s , this corresponds to Sinces ′ ∝ β, increasing the 5th force coupling strength delays this crossing (i.e. it decreases | log a cross |) if all other parameters are kept fixed. On the other hand, in the limit m φ ≪ H 0 one would find the scaling 1 ≪s ini ∝ √ βH 0 /m φ , leading to a cross ∝ (m φ /H 0 ) 2/3 . Hence, a smaller m φ /H 0 corresponds to an earlier crossing. Equation (3.17) is an excellent estimate of the time when the transition takes place, as demonstrated by the second row of panels in Fig. 5.
After the transition to w s ≃ −1, initially the potential term is still negligible in the KG equation and the approximate solution (3.10) still applies, givingρ s ≃ m 2 φs 2 /(2G s ) ≃ const, modulo a logarithmic evolution. Eventually, the potential energy dominates in the KG equation too and s effectively becomes a decoupled quintessence field. This new regime, however, does not modify the scaling of the s background energy density because the fluid equation (3.4) now has the decoupled solutionρ s ∝ a −3(1+ws) ≃ const , since w s ≃ −1.
Finally, once Ω s becomes of O(1) the dynamics follows well-known solutions from (decoupled) quintessence models, see e.g. Ref. [57] for a review.

Numerical results
Numerical solutions for background quantities are obtained through implementation in the Einstein-Boltzmann solver CLASS [55,56]. In all cases the input parameters include {β, m φ /H 0 , Ω d } beyond the ΛCDM ones. Two distinct shooting procedures are applied depending on the scenario considered: for pure 5th force,s ini is set to 10 −4 and ω Λ is determined by imposing the closure condition today; for Coupled Dark Energy, ω Λ is set to zero ands ini is determined from the closure requirement. The background solutions have also been cross-checked using a Mathematica code where the evolution equations were transformed into a system of first-order ODE, as done in Ref. [58] by generalizing previous methods [59].
Figures from 1 to 5 show the evolution of background quantities, focusing on four benchmark cosmologies (beside ΛCDM): for the case of 5F we set m φ /H 0 = 0.1, taking β = 0.005 or 0.02; for CDE we take m φ /H 0 = 1 and again β = 0.005 or 0.02. Anticipating the results of Sec. 5, the smaller choice β = 0.005 corresponds approximately to the 95% c.l. upper bound we have set in the 5F scenario using Planck and BAO data. On the other hand, the stronger coupling β = 0.02 is well within our bounds in the CDE scenario and is instructive for the 5F case as well, where it qualitatively represents the larger values β ∼ 0.01 that may resolve the Hubble tension. For concreteness, in this section we set all standard ΛCDM parameters to their Planck best fit values [4]. In addition, we take The second row of Fig. 5 confirms that in the CDE scenario the transition from w s = +1 to w s ≈ −1 takes place during late MD, as predicted by Eq. (3.17). Furthermore, a departure from w s = +1 is observed in the 5F scenario as well, during the CC-dominated phase. This is a consequence of the choice m φ /H 0 = 0.1; for smaller mass of the scalar, w s ≈ +1 would persist until today. Finally, the last row of Fig. 5 shows the comoving distance χ(z) = z 0 dz ′ /H(z ′ ), normalized to ΛCDM. For 5F we observe the expected high-z enhancement, while for CDE the comoving distance is close to its standard value.
The numerical implementation of the background equations also reveals the presence of an upper bound on β in the 5F scenario, simply coming from the existence of a solution to Einstein's equations until z = 0. For cosmological parameters close enough to the Planck ΛCDM best fit values, increasing β to O(0.1) reduces so much the energy density in the χ fluid that the Universe never arrives at a CC dominated phase. The scalar field s becomes dominant and quickly overcloses the Universe before any accelerated expansion can begin. This result highlights the importance of consistently including the effects of new degrees of freedom both at the level of the background and of the perturbations, even though the latter are naively considered to be most affected by the new dynamics. We will see another example of the importance of a consistent background evolution in Sec. 4.2.1.
An analogous bound does not exist in the CDE case. A larger β is always compensated by a larger initial value ofs, yielding in all cases a smooth transition from a matter dominated to a CDE dominated Universe.

Cosmological Perturbations
Similarly to the cosmological background, the equations of motion at linear order in perturbation theory can be derived from the general equations for the DM fluid discussed in Sec. 2.1 and the ones for the scalar fifth force and the metric in Sec. 2.2. Expanding to first order the continuity equation (2.7) and Euler equation (2.8) for the DM fluid we obtain where ∇ i → ik i in Fourier space and we have made the definitions δ x ≡ δρ x /ρ x and θ x ≡ iv i x k i . The KG equation for the scalar fifth force in Eq. (2.12) expanded at first order is 7 As we have seen in the previous section, at the background level the presence of the scalar fifth force impacts the Hubble flow in a highly non-trivial way, which also depends on whether the scalar is accounting or not for the DE in the Universe today. Now, from the equations for the DM perturbations one actually recognizes the imprint of the new force as a modification of the strength of the gravitational potential. In Eq. (4.2) we see how the velocity field is driven by the spatial derivatives of the gravitational potential Ψ which is augmented by the presence of the 5th force potential (∂ log m χ (s)/∂s) δs. This effect dominates over time derivatives of the potential for modes whose wavelengths are much shorter than the horizon size, generating relative velocity and density perturbation between DM and the baryons that will grow over time. This is in striking contrast with the ΛCDM cosmology for which, given adiabatic initial conditions, relative velocity perturbations decay with time and relative density perturbations stay, at best, constant after recombination (ignoring reionization). Interestingly, the time dependence of the s background affects the behavior of the fluctuations, in three different ways: i) in the continuity equation (4.1) a term proportional to δs is present and controlled bys ′ ; ii) similarly, in the Euler equation (4.2) a term proportional tos ′ enters the friction term; iii) moreover, the Hubble rate H depends nontrivially ons, and this dependence cannot be dropped for a meaningful prediction of the CMB or matter power spectrum. For modes with wavelength much larger than the horizon, the terms proportional to k 2 in Eq. (4.2) can be dropped, yet the explicit dependence of the equations ons ′ remains together with the implicit one coming from H. This implies that we expect differences with respect to a ΛCDM Universe even on super-horizon scales.
In the KG equation (4.3) the first line contains the standard terms for an uncoupled scalar, whilst the second line contains the new terms generated by the interaction of s with the DM. Deep inside the horizon, the scalar field satisfies a Poisson equation of the form to be compared with the corresponding equation for the Newtonian potential, as derived from the perturbed Einstein equations (4.6) and (4.7) given below. Thus, at the fluid level, the strength of the new interaction is dictated by G s (∂ log m χ (s)/∂s), rather than by G s alone. For example, for our representative model one has G s (∂ log m χ (s)/∂s) = G s /(1+2s). Sinces does not evolve much over time, we could therefore fine-tune the initial condition ofs to some very large value, effectively decoupling the DM from the scalar field. This is indeed what happens if one wants a very light and coupled scalar, m φ ≪ H 0 , to play the role of DE at late times. In this regime the energy density in the scalar field behaves like ∼ m 2 φ M 2 Pls 2 ini /β, and to reproduce the energy density in DEs 2 ini ≫ β is needed, implying that the strength of the fifth force effectively goes to zero. For this reason we have chosen the largest possible value of the mediator mass, m φ = H 0 , as our CDE benchmark.
Finally, we write Einstein's equations in the form where σ x are the shear perturbations. To solve these equations we need to express the perturbed fluid variables for the scalar field in terms of the perturbed field variables [60], 8) while no shear perturbation is generated, since the scalar field is minimally coupled.

Initial conditions
To set up the initial conditions (IC) for the perturbations in Newtonian gauge we follow closely Ref. [43]. We assume radiation domination and solve for all the relevant variables by expanding in kτ ≪ 1. The IC for the metric perturbations and for all other species, except the DM and the scalar field, are unchanged and reported in Appendix B. In this section we describe how to consistently set up the IC for χ and s at the lowest order in kτ . Our implementation in CLASS includes the IC in synchronous gauge as well, which are discussed in Appendix B.
The simplest kind of IC are the so called adiabatic IC in the density perturbations, which can be obtained by requiring that the gauge invariant relative entropy perturbations between all species vanish: This assumption for the IC can be realized for example in models of single field inflation. Imposing Eq. (4.9) to the photon and DM fluids at lowest order gives which very deep in radiation domination reduces to δ χ = 3δ γ /4, sinces ′ τ ≪ 1 as shown by Eq. (3.8).
The same procedure could be repeated for the fluid perturbations of s with respect to the photon fluid, obtaining (recall that δ s ≡ δρ s /ρ s ) At early times we can safely assume w s = 1 and the ratio appearing in the denominator of the second term, which was already found in Eq. (3.9), is equal to −2/3. Thus δ s = δ γ /2 to lowest order. 8 We can then use the definition of the perturbed fluid variables in Eq. (4.8), together with the physical condition that the velocity divergence of s vanish when τ → 0 , to obtain the IC for the field, δs = −δ γs ′ τ /4, and for the velocity, θ s = − δ γ k(kτ )/4 = θ γ . Finally, we plug the solution for δs into the equation for θ χ and obtain θ χ = θ γ to lowest order, once again because the long range interactions are strongly suppressed bys ′ τ ≪ 1.
We have seen in Sec. 3.2 that the background dynamics of the fifth force field at early times is fully determined by its interaction with the DM. We can then ask whether the same is true at the level of the perturbations, and therefore if the IC for δ χ and θ χ fix the ones for the scalar field. We now show that this is the case. Indeed, for adiabatic IC, imposing Eq. (4.11) is redundant and the adiabaticity of s follows from the one of χ. Deep in radiation domination and outside the horizon, the KG equation reads where in the adiabatic case δ χ + 2Ψ = −δ γ /4 is constant in time. The general solution to this equation is where (δs) ad = − δ γs ′ τ /4 is the adiabatic piece and α na is a constant parametrizing the deviation from adiabaticity. The adiabatic term is found by inserting the solution into Eq. (4.8) and again imposing that the velocity divergence of s goes to zero when τ → 0. It is then trivial to check that this implies δ s = δ γ /2 and therefore the adiabaticity of s. This is true only at leading order in the kτ expansion, with the subleading corrections easily calculated following the same steps outlined above, see Eq. (B.11).
It is also interesting to study non-adiabatic scalar field fluctuations. In this case the α na term in Eq. (4.13) generates an additional velocity contribution for χ, whose IC now reads, neglecting terms suppressed bys ′ τ ≪ 1, with (θ χ ) ad = θ γ . We see that the non-adiabatic piece is proportional to ∂ log m χ (s)/∂s and is therefore suppressed ifs is large, as it happens in CDE scenarios with m φ ≪ H 0 . In addition, the initial s density perturbation is modified to where (δ s ) ad = δ γ /2 . The IC for the DM density perturbation δ χ is unchanged. We note that α na is in general a function of k, which depends on inflationary dynamics. This is analogous to the initial curvature perturbation R(k) = 2 C(k) [43], which controls the spatial behavior of the adiabatic IC.
In the presence of a dark long range 5th force we cannot therefore a priori assume that the velocity of DM and of all other species is the same even on super-horizon scales. A non-zero α na will alter dramatically the large scale evolution of cosmological perturbations, which in ΛCDM is based on the constancy of the comoving curvature perturbation outside the horizon, and as we will see in Sec. 5 it is severely constrained by current data.

The evolution of density fluctuations
As we did for the cosmological background, we have implemented to first order in perturbation theory all the relevant equations and their IC in CLASS. As a robustness test of our code we have checked that the Newtonian and synchronous gauges give identical results for physical observables, such as the CMB power spectrum. Figure 6 shows the numerical evolution of the density perturbation of χ for four different modes, going from the largest observables scales (k = 10 −3 Mpc −1 , upper left panel), to small scales (k = 1 Mpc −1 , lower right panel). To produce these plots we use the same benchmark cosmological parameters as in Sec. 3, setting in addition α na = 0 and the same IC for the perturbations in the baseline ΛCDM and in our models. We find that at early times, when modes are outside of the horizon, the evolution of DM density fluctuations is almost identical to ΛCDM, regardless of the wavelength. This indicates that the dynamics outside the horizon is still mostly set by the IC, at least in the adiabatic case. Moreover, since at the background level the evolution during radiation domination of DM and baryons is unchanged, we expect little to no difference, with respect to a standard scenario, on which scales enter the horizon before matter-radiation equality. In addition, modes that enter the horizon deep in the radiation era will evolve like in the ΛCDM case until equality, as we can see from the lower panels in Fig. 6. As soon as matter domination begins, perturbations can grow. Thanks to the presence of a dark fifth force, the fluctuations in DM grow more rapidly than the corresponding ones in ΛCDM. In particular we notice that, at least in the 5F scenario wheres ≪ 1 always, the excess power is much larger than the naive O(β) expectation, similarly to what happened at the level of the background. In the CDE case the effect of long range interaction is less dramatic, due to the extra suppression of the interaction by ∂ log m χ (s)/∂s, as further detailed in Sec. 4.2.1.
Finally, we note that at very late times, in either the CC or CDE dominated phase, the growth rate of DM perturbations is reduced compared to ΛCDM. At first, this might seem to clash with the intuition that a new attractive force should only increase the clustering of DM, but it can be understood by considering the effect of the new interaction on the cosmological background. In the 5F case it follows from the larger value of the CC required by the flatness constraint, which exponentially suppresses the growth of structure. In the CDE case it is a consequence of the faster expansion rate at very late times, see Fig. 5. In the CDE models, the z = 0 amplitude of density fluctuations can even be suppressed with respect to a standard cosmological scenario.

Sub-horizon solutions
While an analytical understanding of all the different stages of evolution of the DM perturbations shown in Fig. 6 is not possible, some regimes can be investigated perturbatively for small β. This is the case of the sub-horizon evolution in the matter dominated era. In particular, we would like to quantitatively understand why the effect of the fifth force is much larger than the simple O(β) counting.
During matter domination and in the sub-horizon regime k/H ≫ 1, we can reduce the evolution of DM and baryons to a system of coupled second-order differential equations: whereρ m ≡ρ χ +ρ b , f χ ≡ρ χ /ρ m , and Ω m ≡ρ m /ρ tot . 9 Notice that all these quantities are time-dependent, with Ω m ≈ 1 deep in MD. The equation for δ χ is modified in two 9 In the special case of exponential form of the field dependent mass, mχ(s) = mχe − √ 2βs , Eqs. (4.16) recover the expressions given in Ref. [36].
ways relative to ΛCDM: i) the friction term receives O(β) corrections due to background modifications; ii) the last term contains a new contribution to the potential, which can be traced back to Eq. (4.4). Exact analytical solutions to Eqs. (4.16) are difficult to obtain, due to the complicated time-dependence of the background. However, the system can be solved perturbatively in β. It is convenient to take as new variables where the first is the total matter perturbation and the second the relative perturbation. Expanding as δ i = δ where H is expanded up to at most O(β) and reads as derived from Eq. (3.14) for the 5F scenario. These equations are valid for modes that enter the horizon well before matter-radiation equality. To solve for the relative perturbation it is sufficient to retain the leading order expression of Hubble, yielding the solution which grows linearly with the scale factor. The relative velocity between Dark Matter and baryons also grows accordingly. On the other hand, to solve for the total matter density up to first order we need to include the first correction to H. Plugging this into Eq. (4.19) and expanding as δ m = δ m , we arrive at the differential equation with solution The strong impact of background corrections is further highlighted by Fig. 7, demonstrating the good agreement of our analytical solution Eq. (4.24) with the numerical result. For completeness we also report the individual solutions for χ and the baryons during MD, Thus far, in the discussion of sub-horizon solutions we have focused on the 5F scenario.
For the CDE case one can derive analogous analytical solutions, but the increase of power is far less pronounced.

The CMB and the total matter power spectrum
The evolution of density perturbations over time and at fixed wavelength, discussed in the previous section and shown in Fig. 6, is useful to understand the physical effects generated by the fifth force, but it is not directly observable. In particular, in this work we are interested in constraining new long range forces with CMB and BAO data.  dashed lines to β = 0.02. The color coding is also the same as before, with the prediction for the 5F scenario in light blue, and the one for CDE in light orange. For reference the corresponding ΛCDM prediction is shown in gray. The shape of the CMB primary power spectra is mainly determined by the physical scales relevant for the tight coupling between baryons and photons, and by their projection onto observed angles on the sky. As we have seen in Sections 3 and 4.2, before matter-radiation equality a Universe containing a fifth force with small β is for all practical purposes indistinguishable from a ΛCDM one. Since last scattering happens close enough to the onset of the new interaction, we expect small differences in the physical scales at recombination, which are governed by the values of ω b , ω χ , and ω γ . For example, the sound horizon changes only by O(β). On the other hand the projection on the sky depends on the angular diameter distance between z = 0 and last scattering, and it is severely affected by the new force, see Fig. 5. To first approximation we thus expect, compared to the ΛCDM case, a shift in the location of the peaks and troughs of the CMB power spectrum. This is indeed what we see in the bottom panels of Fig. 8: the residuals around the standard model oscillate around unity. At sufficiently large multipoles we also start seeing the difference in the angular projection of diffusion damping, and therefore an increase (decrease) of power in the 5F (CDE) scenario.
The difference at large angular scales in the CMB power spectra is again due to projection, more specifically to the Integrated Sachs-Wolfe effect (ISW). While the Bardeen potentials in our model are never constant in time, even during matter domination, it is still the case that the largest difference with respect to a ΛCDM Universe happens at low redshift, which maps the late ISW anisotropies to large angular separation.
Finally, the lower right panel in Fig. 8 shows the CMB lensing power spectrum. The amplitude of the lensing power spectrum depends on the total amount of matter in the Universe, which is reduced compared to a standard cosmological model by the long range force, as shown in Fig. 5. This will act towards decreasing the overall amplitude of the lensing power spectrum. On the other hand, the CMB lensing kernel peaks at relatively high redshift, at which the clustering of DM is enhanced by the new force. This effect partially cancels the lower matter density, and it results in a relatively small difference in the lensing power spectrum when compared to a ΛCDM one. It is also worth pointing out that different multipoles L receive contributions from different physical scales k, e.g. larger multipoles are mostly sourced by smaller scales, hence the scale dependence we see in the lensing power spectrum in Fig. 8 is compatible with the one we observed in Fig. 6. It should be kept in mind that the numerical evaluation of the lensing power spectrum requires a prescription for the fully nonlinear matter power spectrum. This is not available for the models described in this work, hence the precise value of the high-L lensing power spectrum should be taken with a grain of salt. For CDE models some results have appeared in Refs. [15,58,62], but they are not yet at the same maturity level of the prediction in the standard model. Figure 9 shows the linear matter power spectrum for different redshifts. While not directly observable, the total matter power spectrum is closely related to the galaxy power spectrum measured in redshift surveys and to the cosmic shear and galaxy-galaxy lensing measurements of imaging surveys. Roughly speaking, we expect the difference in the matter power spectrum between a ΛCDM model and a cosmology with dark long range force to be twice the effect we saw in Fig. 6 and Sec. 4.2. For this reason we will not repeat here the same discussion, but rather highlight a few key features of the models described in this work. First of all, as presented in Sec. 4.2.1, the effect of the new force is much larger than the naive O(β) counting. On the other hand, the new long range interaction primarily changes the amplitude of the power spectrum, with only a mild scale dependence, roughly of O(β). In comparison with a ΛCDM model, and at fixed H 0 , the excess power is larger at higher redshift, and then decreases towards z = 0. As discussed in Sec. 4.2, this is due to the difference in the expansion rate at very low redshift between a Universe with a fifth force in the dark sector and a standard one. Finally, in light of the discussion in the next section on the use of BAO data, it is useful to look at the amplitude of the relative density and velocity power spectra. Figure 10 shows, at z = 0 and in the ΛCDM case, the matter power spectrum (gray dotted line), the cross power spectrum between the total matter and the relative density perturbations, ⟨δ m δ r ⟩ (gray solid line), and the cross power spectrum between the total matter and the relative velocity divergence perturbations, ⟨δ m θ r ⟩ (gray dashed line). 11 As it is well known, in a ΛCDM model relative density perturbations are at most 1% of the total matter and relative velocities are negligible. In light blue and light orange colors we show the same cross power spectra for the 5F and CDE scenarios with β = 0.02. The cross power spectrum between the total matter and the relative density perturbations in cosmologies with a dark fifth force,  Figure 10. Cross power spectra at z = 0 between the total matter and relative density perturbations δ r (solid lines) and between the total matter and relative velocity divergence perturbations θ r (dashed lines). The total matter power spectrum in ΛCDM is also shown for reference (dotted line). The light blue, light orange, and gray lines show the 5F, CDE, and ΛCDM scenario, respectively. solid lines, is approximately an order of magnitude larger than the corresponding quantity in the ΛCDM case, and it can reach 10% of the total power spectrum. The BAO oscillations are also out of phase with those in the matter power spectrum. The cross power spectrum between the total matter and the relative velocity divergence perturbations, dashed lines, is O(100 -1000) times larger than in the standard cosmological model. For comparison, at z = 10 the relative density perturbations in cosmologies with dark fifth forces are still more than two times larger than the corresponding quantity in a ΛCDM scenario, which becomes a factor of 30 for the relative velocities.

Constraints and Discussion
In this section we present our constraints on the cosmological models discussed in this work. We employ a Markov Chain Monte Carlo approach, and in particular the Metropolis-Hastings algorithm, to scan the parameter space of interest until the standard criteria for the convergence of the chains are reached. Our implementation relies on the MontePython code [63,64]. Datasets used in this work include the most recent Planck temperature and polarization data of the CMB [4,65], and measurements of the Baryon Acoustic Oscillations scale in spectroscopic galaxy surveys [5][6][7][8], to which we collectively refer as BAO. We will also briefly comment on the possibility that long range forces in the dark sector could alleviate the well known tension between the CMB and local determinations of H 0 (see for example Refs. [66][67][68] for extensive reviews of the problem and of the proposed solutions to it). We will therefore show additional constraints that include a Gaussian prior on the Hubble constant, H 0 = 74.03 ± 1.42 km/s/Mpc at 68% c.l., from Ref. [69]. 12 We vary 5 ΛCDM parameters {ω b , H 0 , n s , A s , τ }, where τ is the reionization optical depth, to which we add the dark sector parameters { Ω d , β}. Regardless of the combination of datasets, we find that most ΛCDM parameters are unchanged, with the exception of the Hubble constant. We will therefore only show the constraints for Ω d , β, and H 0 . Due to the marked differences in the time evolution of the background and perturbations between a 5F cosmology and a CDE one, we present the parameter constraints separately, in Sections 5.1 and 5.2 respectively. We refer the reader to Sections 3 and 4 for a detailed comparison of the two scenarios and recall here only the key distinctions.
In the 5F scenario the cosmological energy density in the scalar field is negligible at all times. This is obtained for initial field background valuess ini ≪ 1, which in turn implies that ∂ log m χ (s)/∂s ≃ 1. The dependence on the scalar mass m φ is very weak provided it is smaller than H 0 . It is therefore the fifth force coupling β that primarily controls the size of the differences with respect to the ΛCDM model.
In the CDE scenario the scalar field is required to provide the accelerated expansion at very late times. This implies a non-trivial interplay betweens ini , m φ and β, as picking a value for two of them fixes the third one. At fixed β ≃ 10 −2 and m φ = H 0 , we find s ini ∼ O(1) and therefore ∂ log m χ (s)/∂s < 1, as seen in Sec. 3. If m φ is reduced below H 0 , thens ini increases to keep the potential energy in the scalar field the same, which in turn yields a stronger suppression of ∂ log m χ (s)/∂s ≪ 1. Thus, in a CDE Universe the effective strength of the fifth force is generically reduced compared to the 5F case, hence weaker constraints on β are expected.
In Sec. 5.3 we present the constraints on the parameter α na , defined in Sec. 4.1, which quantifies departures from the adiabatic initial conditions for DM velocity perturbations. Figure 11 shows the constraints on our model parameters in a 5F scenario with m φ /H 0 = 0.1. Different colors display different datasets, and for comparison we plot as dashed lines and contours the corresponding bounds in ΛCDM from CMB-only data. Notice that in the standard model Ω d simply reduces to the present day value of the CDM energy density, Ω CDM . Using Planck data alone (blue lines and contours in Fig. 11) the bound on the strength of the fifth force normalized to gravity is β < 0.011 at 95% credible levels (c.l.). This is quite a strong constraint considering that the cosmological model at recombination is mostly unchanged with respect to the standard scenario (see Sections 3 and 4), and thus most of the constraining power comes from the geometric projection of physical scales at last scattering onto observed angles on the sky. 12 We note that more recent measurements of H0 from the SH0ES + Pantheon teams have been presented in [70,71], but they would not change quantitatively our results. As emphasized by the SH0ES team and by others [72][73][74], a more robust way to implement the local distance ladder would be through a prior on the Supernovae absolute magnitude M b , and then to refit the Pantheon data assuming a given background cosmological model. We will return to these issues in a future publication.

Pure Fifth Force and the Hubble tension
Degeneracies between β and both H 0 and Ω d are clearly visible, and can be understood in the following way. In a flat ΛCDM model, keeping the physical density of photons, baryons and CDM fixed, a change in the value of the Hubble constant can be compensated, to ensure flatness, by a different value of the CC energy density Ω Λ , which in turn implies a different value of the present day matter density and hence a strong degeneracy in the Ω m -H 0 plane. This geometric degeneracy is not exact in flat models, and it is internally broken by CMB data, as shown by the dashed contour in Fig. 11. In particular one cannot take arbitrarily low values of Ω m in a ΛCDM model and still fit the data. In the presence of a dark fifth force, and for fixed values of the energy densities at recombination, we are now allowed to increase the Hubble constant to keep the distance to last scattering the same as in ΛCDM, see the light blue lines in Fig. 5. This explains the positive correlation between H 0 and β we observe in Fig. 11. The other degeneracy lines are now just consequences of this fact, because to ensure flatness with a larger value of H 0 , a larger Ω Λ , i.e. a lower Ω d , is needed. Percent level strengths for β correspond to very low values for Ω d , and the resulting degeneracy between the present day DM density (approximately given by Ω d ) and the Hubble parameter is much larger than in the standard cosmological model. Using only CMB data, the posterior distribution for all three parameters is quite non Gaussian, which, combined with the larger degeneracies discussed above, results in large errorbars for the parameters of interest. In particular we find 66.40 < H 0 /(km/s/Mpc) < 72.90 and 0.226 < Ω d < 0.277 at 95% c.l..
There are several ways to break geometric degeneracies in primary CMB data. The most constraining one is to include BAO data. BAO measures, in units of the sound horizon at the baryon drag epoch r d , the Hubble parameter and the angular diameter distance to the redshift of a given galaxy sample, therefore allowing to constrain the matter density Ω m independently from the CMB. BAO measurements are robust to changes in the background and are pretty insensitive to the broadband shape of the galaxy power spectrum and correlation function [76][77][78].
However, all BAO measurements assume that the way r d is imprinted in the distribution of galaxies follows the distribution of the total matter overdensity δ m . Equivalently, relative density and velocity perturbations between DM and baryons are typically neglected in BAO analyses. It was long ago realized that this approximation might lead to systematic biases in BAO analyses [79,80], since in ΛCDM the oscillatory features in δ r and ⃗ v r are out of phase with the ones in δ m . Of particular concern is the fact that relative perturbations enter the prediction of the galaxy power spectrum multiplied by unknown free parameters, which can easily be of order one or larger [81,82]. As said earlier in Sec. 4, in the ΛCDM model relative density perturbations do not grow and relative velocities decay, and recent work in Refs. [61,83,84] indicates that, for current surveys and in the standard model, biases to the distances inferred via the BAO method are smaller than the typical measurement error. However, in the presence of long range interactions acting only on one species, like the ones we consider here, relative perturbations actually grow with time, as discussed in Sec. 4 and shown in Fig. 10. As a benchmark, for β = 0.005, which is roughly the 68% c.l. bound from Planck data, the relative densities are approximately a factor of 3 larger, on BAO scales and at z = 0, than the corresponding fluctuations in the standard model (a  Figure 11. Bounds on the pure 5th force scenario with m φ /H 0 = 0.1 ands ini = 10 −4 . The blue shaded region is favored by Planck CMB data [4,65] at 1σ and 2σ. The orange shaded region also includes BAO measurements [5][6][7][8], while the red shaded region adds CMB lensing measurements from Planck [75]. The green shaded contours combine Planck CMB data with a H 0 prior from local measurements [69], whereas the purple region also includes lensing. The dashed contour shows the 2σ region favored by Planck CMB data in ΛCDM. factor of 10 at the turnaround of the power spectrum). The relative velocity divergence, θ r , is now a factor of 100 larger than in ΛCDM, again at z = 0. At k = 0.1 Mpc/h, we find δ r ∼ 0.02 δ m for β = 0.005, and the ratio increases linearly with β. For reference, the current best BAO measurements have few percent accuracy [7].
On the other hand, these new relative fluctuations are, to first order, proportional to δ m , see Eq. (4.21), and are thus not expected to produce additional phase shifts of the BAO. They add to the off-phase piece produced at recombination, and eventually dominate over the latter after redshift z ∼ 10 -15. Making these statements more quantitative would require a dedicated BAO analysis framework that includes the effects of relative perturbations in the ΛCDM model and beyond. Such pipeline is not currently available, and preparing one certainly goes beyond the scope of this work. Given the numbers in Ref. [61], we roughly expect that β ∼ 0.01 is the ballpark value at which relative perturbations effects could become significant at BAO scales. This corresponds to the 95% c.l. upper bound on β from CMB data alone in the 5F scenario.
It is therefore reasonable, albeit with the caveats mentioned above, to combine Planck data, which do not favor a percent value of β, with BAO data to further improve the bounds.
The constraints from Planck plus BAO are shown in orange in Fig. 11. As expected, BAO break the degeneracy between H 0 and Ω d , both parameters move closer to their ΛCDM counterparts and are much more Gaussian distributed. The allowed range for the coupling strength also shrinks, with now β < 0.0054 at 95% c.l.. To date, this is the strongest bound on dark fifth forces from cosmological data, corresponding to g D < 2 × 10 −32 m χ /(10 −3 eV) when expressed in terms of the coupling constant of the underlying field theory description.
Given the tail at large values of the Hubble constant in the H 0 posterior from Planck data alone, we can also ask whether a dark fifth force can provide a solution to the Hubble tension. We therefore run the CMB likelihood including a prior on H 0 from Ref. [69]. The results are shown in green in Fig. 11. In the Ω d -H 0 plane the green contours occupy the leftmost region that was allowed by the CMB-only fits, i.e., a large H 0 requires small Ω d . The best fit value for the Hubble constant is H 0 = 72.41 ± 1.40 km/s/Mpc at 68% c.l., clearly compatible with the local measurements. Interestingly, there is now evidence, at more than 3.5 σ level, for a non-vanishing 5th force strength in the dark sector, with β = 0.0102 ± 0.0028 at 68% c.l.. Clearly, such large values of β results in very different distances to the mean redshift of galaxy surveys and in a different matter power spectrum at late times compared to a ΛCDM model, see for example Fig. 9. However, as discussed above, a rigorous analysis of late time data, such as the BAO and the Full Shape of the galaxy power spectrum, in presence of nonzero long range forces requires a number of new tools that are yet to be developed. For these reasons we decided not to combine a local H 0 prior with BAO data. We intend to return to these issues in a forthcoming publication, hoping that our results will further motivate the community to develop more general analysis pipelines of Large Scale Structure data.
Another way to partially lift parameter degeneracies is to include CMB lensing information. This comes with the caveats mentioned in Sec. 4 about matter non linearities in the presence of a fifth force. However, given the current uncertainties of the Planck lensing power spectrum and the relatively tight constraints on β from primary CMB, we do not expect large biases in the inferred parameters by adding Planck lensing data. This might not be the case for ground based CMB experiments, targeting the lensing power spectrum at much smaller scales [85]. In a ΛCDM model, CMB lensing is primarily sensitive to a combination of H 0 , Ω m and σ 8 , the amplitude of matter fluctuations at 8 Mpc/h. A measurement of the angular size of the sound horizon at last scattering with primary CMB data then breaks this degeneracy, which in turn helps to break the geometric degeneracy in the primary CMB. The inclusion of the Planck CMB lensing power spectrum data on our constraints is shown in Fig. 11 by the red lines and contours, for the primary CMB plus BAO combination, and in purple for the primary CMB plus local H 0 prior case. Unfortunately, in the extended parameter space of our models, Planck lensing information does not yet provide meaningful improvements of the uncertainties over other datasets, as one can see from the minor difference between the orange and red lines, or between the green and purple ones.
Finally, since we discussed a possible explanation of the Hubble tension by a dark fifth force, a comment is warranted about the impact on σ 8 . As is well known, this is affected by a moderate but persistent discrepancy between CMB and Large Scale Structure 95% c.l. measurements, with the latter giving smaller results [86,87]. At face value (in particular, for fixed bias parameters), in the 5F scenario studied here the σ 8 tension is expected to worsen.

Coupled Dark Energy
Similar considerations about degeneracies between the different parameters apply to the bounds in the CDE scenario. In this case the effective strength of the fifth force, which contains powers of ∂ log m χ (s)/∂s, is suppressed by the fact thats ≳ O(1) in order for the new mediator to play the role of the Dark Energy at late times. We thus expect the bound on β to be weaker than in the 5F scenario. Figure 12 shows the constraints on our three reference parameters in a CDE Universe with m φ /H 0 = 1. From the discussion in Sections 3.2 and 4 it is clear that the bound on β weakens for smaller mass of the scalar field, hence Fig. 12 corresponds to the most constrained CDE scenario. The blue lines and contours show results using only CMB data. We find that the contours are closer to their ΛCDM counterparts than the ones for 5F. In particular, the H 0 posterior does not have a tail at larger values, implying we should not combine CMB data with a local prior on the Hubble parameter. The bound on β is significantly weaker than in the 5F case, with β < 0.034 at 95% c.l.. Adding BAO data does not improve the constraints as much as in the 5F case, again because the tuning of the initial conditions to reproduce the CC at late times renders the effect of the fifth force vanishing small, see Figs. 8 and 9. The upper bound on the coupling parameter becomes β < 0.029 at 95% c.l.. We do not show results including lensing information, because this does not add significant constraining power. A summary of the main parameter constraints, for both the 5F and CDE scenarios, is provided in Table 1.
It is worthwhile to comment on how our CDE scenario differs from previous studies of DM-DE interactions. Formally, our setup can be described within the parametrization developed in Refs. [32,33], 13 but it is not straightforward to compare the numerical results 13 In the notation of Refs. [32,33], our CDE model corresponds to as well as αM = αB = αT = 0, reflecting the fact that gravity per se is not modified. Furthermore, we have the identifications π = δs/ṡ and vx = −a θx/k 2 , and Φ there = Ψ and Ψ there = −Φ.  [4,65]. The orange shaded region also includes BAO measurements [5][6][7][8]. The dashed contour shows the 2σ region favored by Planck CMB data in ΛCDM.
of those works to ours, because of the assumptions made there about the time evolution of the cosmological parameters [33]. In particular, Hubble was taken to evolve as in wCDM (with constant w) and the coupling β γ was taken to be time-independent, neither of which applies in our framework. In general, we emphasize that the simplest microscopic model discussed here -just an ultralight scalar with a trilinear interaction with the DM fieldleads to a complex cosmological evolution, where the dynamics of both the background and the fluctuations are modified.
General forms for the DM-DE interaction were also introduced in Ref. [30], 14 which however limited the discussion to a qualitative illustration of the expected impact on the CMB and matter power spectra. Notice also that the results in Ref. [30] were presented assuming an exponential form of the field dependent mass m χ (s), which is not considered in this work and will be discussed elsewhere. 14 In the notation of Ref. [30], our CDE model corresponds to a Type-1 theory with α = log mχ(φ) and F = Y + V (φ). Note that their perturbed KG and continuity equations (91) and (92) are incomplete: each of them is missing a term ∼ α ϕϕ , given in our Eqs.

Non-adiabatic initial conditions
Finally, we can use cosmological data to constrain α na , which parametrizes possible deviations from adiabatic initial conditions. As discussed in Sec. 4.1, a non-zero α na sources a non zero relative velocity between the DM and the other species on super-horizon scales. We assume α na to be scale-independent and, for simplicity, we only show in Fig. 13 the bounds in the 5F scenario, combining CMB and BAO data. As expected, the constraint in CDE models is weaker because the effective non-adiabaticity turns out to be proportional to α na (∂ log m χ (s)/∂s), see Eq. (4.14). For 5F we find α na = 0.0021 ± 0.0065 at 68% c.l., which is clearly compatible with zero. For such small values, the constraints on the other model parameters are basically unchanged. We can understand this by noticing that β does not modify much the evolution of super-horizon modes, as illustrated by Fig. 6. The constraint on α na is quantitatively of the same order of the ones on other kinds of non-adiabatic initial conditions [4], e.g. isocurvature perturbations.

Interplay with Visible Fifth Forces
So far we assumed the long range fifth force to be entirely confined to the dark sector. However, it is interesting to consider the possibility that the same fifth force couples to the visible matter. In this situation, depending on the specific realization, one can study the interplay between the new bounds on dark fifth forces derived in Sec. 5 and the existing constraints on visible long range forces. Since we focus on scenarios where m φ ≲ H 0 , the mediator can be taken as effectively massless on all the scales relevant for experimental tests of gravity in the visible sector. We assume the DM mass to satisfy the naturalness bound in Eq. (2.20), which requires it to be smaller than approximately 0.01 eV. In this regime the DM is an ultralight boson, possibly produced in the early Universe via the misalignment mechanism [37][38][39].
For m φ ≲ H 0 the effect of the visible fifth force on the interaction between two bodies with masses m A and m B at distance r can be parametrized as where the Yukawa factor e −mφr can be dropped for such a light mediator, so that no deviation from the Newtonian potential would arise. 15 The coupling α A can be related to the field dependence of the macroscopic masses, α A = (4πG N ) −1/2 (∂ log m A (φ)/∂φ). In general this coupling contains a universal part, which results in an unobservable shift of the Newton constant, and a non-universal part which depends non-trivially on the composition of the material and results in effective violation of the EP on macroscopic scales. The latter effect can be measured as the difference in the acceleration of two bodies with different compositions and has been extensively tested experimentally (see e.g. Ref. [3] and references therein). The macroscopic coupling α A can be expressed in terms of the microscopic couplings of the light scalar with the light Standard Model (SM) fields, which can be written in general as where e and g 3 are the electromagnetic and QCD couplings respectively, β 3 = −b 3 g 3 3 /16π 2 with b 3 = (11 − 2N f /3) is the beta function encoding the evolution of the QCD coupling constant with energy, and γ mq are the anomalous dimensions of the quark masses. Notice that a non-canonical normalization has been adopted for the electromagnetic gauge field. The analysis of Refs. [88,89] calculated the leading effects of the microscopic coefficients d i on the macroscopic parameter α A , where d * g = d g + 0.093(dm − d g ) + 2.7 × 10 −4 d e corresponds to the composition independent part and we defined dm = (d mu m u + d m d m d )/(m u + m d ). The coupling of the scalar to the 15 The experimental tests of Newton's inverse square law range from length scales of 10 −6 m to few au and are sensitive to mediator masses down to 10 −18 eV [1]. Below this mass, deviations from the inverse square law decouple like mφLexp where Lexp is the size of the experimental apparatus. This would result in tiny deviations at the 10 −15 level or below for the range of mediator masses considered in this work. electromagnetic field strength, d e , is numerically suppressed in the matching to α A compared to d g and dm, due to the weak dependence of nuclear binding energies on electromagnetism. Q ′m ≃ −0.036 A −1/3 − 1.4 × 10 −4 Z(Z − 1)A −4/3 and Q ′ e ≃ 7.7 × 10 −4 Z(Z − 1)A −4/3 in Eq. (6.3) are material dependent "dilaton charges" (A is the atomic mass number and Z the atomic number) as approximately derived in Refs. [88,89]. Experimental tests of the EP place important constraints on new long range interactions in the visible sector. The best Earth-based limit comes from the Eöt-Wash experiment [90], which set a bound on the parameter combination |d * g (dm − d g + 0.22 d e )| < 5.1 × 10 −11 at 2σ level. An even stronger limit has been obtained by the MICROSCOPE space mission [91,92], which found at 2σ This corresponds to the constraints |d g | < 1.3×10 −6 , |dm| < 4.2×10 −6 and |d e | < 9.9×10 −5 if one coupling is switched on at a time in Eq. (6.2). Again we see that the electromagnetic coupling is subject to a sizably weaker bound. One should also keep in mind that, as pointed out in Ref. [93], the constraints on EP violation can be significantly relaxed if the microscopic couplings conspire to give macroscopic effects that resemble very much the universality of Newtonian gravity (at least for the class of elements that have been tested experimentally so far).
The direct couplings of φ to the SM in Eq. (6.2) induce temporal variations of the SM parameters. These can be tested by comparing the transition frequencies of different atomic clocks, as explored in Refs. [94][95][96]. Since the scalar masses considered here are much smaller than the inverse of the typical run time of the atomic clock experiments (∼ years), the only observable effect induced by the variation of φ is a steady drift in the frequency ratio between the two clocks. For two different optical transitions, the drift rate is controlled by the change of the fine structure constant α, while for one hyperfine microwave transition and one optical transition there is also sensitivity to the change of µ A /µ B , the ratio between the nuclear magnetic moment and the Bohr magneton, which is linearly proportional to the inverse of µ = m p /m e [97].
The drifts in the fundamental constants can be related through Eq. (6.2) to the time evolution of the scalar field, where dm is defined below Eq. (6.3) and M A was estimated in Ref. [97] for several nuclei. Next, recalling Eq. (3.3) we writeφ = − ρ s (1 + w s ), where the negative sign is the correct one in our setup. This allows us to finally express the present day values of the drifts as (see also Ref. [98]) These equations make it apparent that the atomic clock constraints can be weakened if the present day energy density of the scalar field is small (Ω 0 s ≪ 1), or its equation of state is near the CC one (w 0 s ≈ −1). The best experimental constraints on the drifts of the fine structure constant and proton-to-electron mass ratio come from measurements with ytterbium ion clocks and caesium clocks [99,100]. The resulting 2σ bounds (α/α) 0 = (1.8 ± 5.0) × 10 −19 yr −1 and (μ/µ) 0 = (−8 ± 72) × 10 −18 yr −1 constitute improvements by a factor of 90 and 2, respectively, with respect to previous measurements [101,102]. In terms of the couplings in Eq. (6.2), the 2σ bounds read |d e | < 6.1 × 10 −9 where the central values of the measurements were neglected for simplicity.
In the 5F scenario, w 0 s is far from −1 (it is ≈ +1 if m φ ≪ H 0 ) and the present day fraction of energy density in the scalar field is approximately given by Ω 0 s ≃ βf 2 χ /3 , where β is the coupling of the scalar to DM, see Eq. (3.12). In light of the constraints on β from CMB+BAO data presented in Sec. 5.1, we thus expect Ω 0 s ≲ 10 −3 . Accounting for this suppression, we find by comparison with the MICROSCOPE bounds reported below Eq. (6.4) that the sensitivity of atomic clocks to the scalar-photon coupling is stronger than EP tests by at least two orders of magnitude. On the other hand, for the scalargluon coupling MICROSCOPE remains the most sensitive probe. In the CDE scenario, Ω 0 s ∼ O(1) and w 0 s can deviate at the 10% level from the CC equation of state, as shown in Fig. 5, only leading to a mild weakening of the atomic clock constraints. In this case atomic clocks always dominate over EP tests.
In the remainder of this section we discuss in explicit scenarios the interplay of the MICROSCOPE and atomic clock constraints with the dark fifth force bounds presented in Sec. 5. This depends on how the couplings in Eq. (6.2) scale with β. In Sec. 6.1 we discuss the most generic scenario where the fifth force coupling to the SM is only induced through DM loops, while in Sec. 6.2 we consider a model where the fifth force couples directly to both DM and the SM.

Visible fifth force from Dark Matter loops
We consider a scenario where the fifth force is sequestered within the dark sector and interacts only with DM at tree level. On the other hand, we assume that DM couples to the SM directly and we want to estimate the size of the fifth force coupling to the visible sector induced by DM loops, see Fig. 14. Specifically, having in mind QCD axion models, we consider the simple case of pseudoscalar DM a coupled to photons and gluons at one loop and to the dark fifth force field φ at tree level, a f a G a µν G µν a − g D m a φa 2 . (6.9) Figure 14. Generic interplay between dark and visible fifth forces. According to Eq. (6.9), the light force mediator is coupled to DM at tree level, while the coupling to the SM only arises via DM loops.
In this simple scenario, taking E/N ∼ O(1) a DM loop would induce a coupling between the fifth force and the SM photons and gluons with the approximate size (6.10) where the couplings to the gauge field strengths were defined in Eq. (6.2). If a nonzero value of β were to be found by cosmological measurements, for instance β ≃ 0.01 as motivated by the present Hubble tension (see Sec. 5.1), the constraint on EP violation in the visible sector given in Eq. As a final remark on the parametrics of Eqs. (6.10) and (6.11), we notice that the weakness of the loop-induced visible fifth force is essentially related to the smallness of m a /f a , which is constrained to be less than about 10 −18 if the mass of the axion is fixed on the QCD line m a = 5.7 µeV (10 12 GeV/f a ) [104]. This suggests that models with larger m a /f a could realize parametrically larger visible fifth forces. As discussed in Sec. 2.3, however, heavier DM would require either a symmetry mechanism in the dark sector (such as supersymmetry) to preserve the naturalness of the fifth force mass, or to accept a large amount of fine-tuning. Under this non-minimal assumption, the implications of the existence of a dark fifth force on the direct detection prospects for WIMP DM were previously investigated in Refs. [105,106].

Visible and dark fifth forces from the radial mode
We now consider a simple model where the fifth force couples directly to both DM and the SM, see Fig. 15. The scalar DM χ is an axion, while the fifth force mediator φ is the radial Figure 15. Second scenario illustrating the interplay between dark and visible fifth forces. The light force mediator is coupled to DM at tree level, Eq. (6.13), while the coupling to the SM arises from the effective interaction in Eq. (6.15).
mode of a complex scalar field Φ with Lagrangian where the U (1) symmetry acting on the complex field (Φ → e iα Φ) is both spontaneously broken by the Mexican-hat potential and explicitly broken by the last term controlled by δm 2 . The latter preserves the Z 2 symmetry Φ → Φ * that stabilizes the DM χ. Parametrizing the complex field as Φ = (f + φ)e iχ/f / √ 2 with the vacuum expectation value given by f = m 2 /λ , we obtain Matching to Eq. (2.16) leads to the following parameter identifications, , (6.14) where from the equations above we can see that an extremely small quartic coupling and a super Planckian decay constant are required to realize the dark fifth force mediator as the radial mode of a vanilla axion DM model. For β ∼ 10 −2 and m φ ∼ H 0 we need λ ∼ 10 −123 . The trilinear coupling between the axion DM and the fifth force mediator is given by δm 2 /f = β/2 m 2 χ /M Pl . Furthermore, several interactions are generated beyond those we have considered in our cosmological analysis, including a derivative (φ/f )(∂χ) 2 coupling as well as a φ 2 χ 2 term. However, these are unlikely to substantially change the results obtained in Sec. 5.
We note that, since δm 2 ≫ λf 2 , we cannot view the last term in Eq. (6.12) as a small perturbation to the U (1) invariant Lagrangian for Φ. Therefore, the choice to include only the specific explicit breaking structure (Φ−Φ * ) 2 appears to be rather ad-hoc. In particular, there is no symmetry argument justifying the suppression of the explicit breaking operator (Φ + Φ * ) 2 , which would induce a mass term for φ of order m 2 χ . In this respect, the theory discussed here should be considered as a toy model, of somewhat limited theoretical interest but still useful to compare the parametric sensitivity of different fifth force probes. Now assuming the presence of some electrically charged chiral fermions with Yukawa couplings to Φ, the effective operator is generated with c ψ ∼ O(1) , together with the usual anomalous coupling of the axion DM to photons. The corresponding value of the effective coupling in Eq. (6.2) is d e = c ψ α √ β/[(1 + 2s 0 )π] and the resulting 2σ bound from MICROSCOPE is β < 0.0018 (1 + 2s 0 ) 2 /c 2 ψ . Interestingly, this is comparable to our new bound from dark fifth forces. However, atomic clocks give much more stringent constraints both in the 5F and CDE scenarios, requiring β < 4.5 × 10 −6 1 Finally, we briefly go back to the issue of fine-tuning discussed in Sec. 2.3. The introduction of new couplings of the ultralight scalar field to the SM raises again the concern of how the tiny m φ can be stabilized against quantum corrections. In the simple setup presented here this issue can be explicitly quantified by estimating the one-loop contribution to the quartic coupling λ induced by the non-renormalizable interaction in Eq. (6.15), ∆λ ≃ c 2 ψ α 2 Λ 4 UV /(4πf ) 4 . The ultraviolet cutoff Λ UV corresponds to the masses of the new electrically charged chiral fermions, which must be heavier than the electroweak scale for O(1) charges, thereby generating a severe fine-tuning problem for the scalar field mass. The issue of introducing testable fifth forces coupled to the SM without fine-tuning is a longstanding one, as reviewed e.g. in Ref. [95]. At present, a convincing dynamical solution does not seem to exist, though interesting attempts have been made [14,49,107,108]. Remarkably, for a fifth force that is sequestered within the dark sector and whose couplings to the SM arise only from DM loops, as discussed in Sec. 6.1, this fine-tuning issue is solved by taking the DM mass to be sufficiently light, as quantified in Eq. (2.20).

Outlook
This work starts a systematic investigation of the phenomenology of dark long range interactions in cosmology. Here we briefly summarize our findings and present a number of future directions.
In this first exploration we focused on scalar dark fifth forces, and derived the corresponding equations for background and perturbations. We focused on the cosmology of dark fifth forces that can be mapped to minimal microscopic theories with natural parameter choices. The analytical results we derived in Sections 3 and 4 show the build-up of the effects of dark fifth forces over time, and clarify why cosmology is such a powerful probe of long range dynamics in the dark sector: even a parametrically small effect can become relevant after more than 13 billion years.
A very light scalar mediator could also account for the Dark Energy, with a very different phenomenology compared to the case where the new degree of freedom purely mediates a new interaction. We dub the former the Coupled Dark Energy (CDE) scenario, and the latter the pure 5th force (5F) scenario. The fundamental parameter of the theory is β, the ratio between the strength of the new interaction and the Newton constant. In the 5F case, we find β < 0.0053 at 95% c.l. by combining CMB and BAO data. This bound weakens by approximately a factor of 8 in the CDE case with m φ = H 0 , due to the fine-tuning required to match the value of the Cosmological Constant. For smaller masses of the CDE field the bound weakens further, scaling approximately at least as m φ /H 0 . These are the strongest constraints on Equivalence Principle violation in the dark sector.
Our work opens up two sets of different but complementary questions: i) on the phenomenological side, it would be important to explore the nonlinear cosmology of dark fifth forces; ii) on the theoretical side, it would be interesting to explore non-minimal scalar theories, vector mediators and the interplay with inflationary models. We now briefly present these future directions in turn.
The rich phenomenology beyond linear theory of our models is yet to be investigated. For example, we made the compelling case for a careful study of the mildly nonlinear regime of structure formation, since we have shown that dark fifth forces could alleviate the wellknown H 0 tension between CMB and local measurements by simply invalidating the use of BAO information. This was possible because the traditional BAO analyses cannot be applied to cosmologies with large density or velocity perturbations between the different species, as it is the case if DM violates the Equivalence Principle. The development of a dedicated BAO analysis pipeline is therefore of the highest priority. More generally, the shape of the power spectrum is much more affected by the presence of new dark long range interactions than the CMB anisotropies. Analyzing galaxy clustering data however requires extending next-to-leading order perturbative calculations of Large Scale Structure beyond the current state of the art. For reference, a value of β = 0.005, close to the current 95% c.l. upper bound from CMB plus BAO data, yields approximately 5% differences in the shape of the power spectrum compared to the ΛCDM case.
Moving to smaller scales, many astrophysical systems could be affected by the presence of a new dark long range interaction, and we list here just a few examples.
At very high redshifts, it is well known that relative velocities between DM and baryons have an important effect on the formation of the first stars and galaxies, and can vastly modify the shape of the 21 cm power spectrum at cosmic dawn. In our models, the relative velocities are one order of magnitude larger than in a ΛCDM scenario, and it would be therefore interesting to study their consequences for the high redshift Universe. Other examples concern the collapse of DM halos and their internal structure. We expect more halos in cosmology with a fifth force, simply because matter can accrete faster.
Their internal structure could also be different than in ΛCDM, and this may be tested for instance with strong lensing data. The abundance and the profile of small mass halos are also related to the number of satellite galaxies, which could be studied assuming a model for the baryonic physics in the presence of a 5th force. Finally, the dynamics of stars and of stellar streams is sensitive to the DM distribution inside galaxies, and could be tested with kinematics data of JWST (the former) and of Gaia (the latter). Our modified version of CLASS provides the baseline for all the above investigations.
The other set of questions are related to the theoretical foundations of the framework. While we discussed initial conditions from the perspective of initializing the linear theory equations, a period of cosmological inflation can provide a dynamical origin for the initial conditions of the Universe. It is therefore interesting to ask how to embed a new very light and interacting degree of freedom into inflation. Actually, a general prediction of inflation models with light scalar spectator fields is the generation of non-adiabaticity in the initial conditions and of Primordial non-Gaussianities, with important observational consequences.
One could also consider more general interactions with scalar fields than the ones studied in this work. One example are dilaton-like couplings, which are interesting because they would evade the suppression of the effective strength of the fifth force in CDE scenarios that we discussed here.
Finally, it is tempting to look into light Abelian vector mediators, whose very small mass can be technically natural. Abelian gauge theories are screened at large distances if there is no net dark charge in the Universe, therefore one would expect a less rich phenomenology than in the scalar case. However, at the level of the perturbations, light interacting vector bosons will generate vector perturbations, which are severely constrained by data. Clearly, if the total dark charge in the Universe is not zero we can anticipate major differences with ΛCDM at all scales.
We hope to hit the road in all these various directions in future work.
DM. The starting point is the microscopic Lagrangian, which can be written as We define the χ fluid as "containing the interaction", so that the equations expressing the non-conservation of the energy-momentum tensor for each species read Notice that the transfer four-vector appearing on the right-hand sides is proportional to the scalar field four-velocity, ∂ ν s ∝ u ν s . Writing the DM stress tensor as [109] T µν and RF denoting rest-frame quantities, the energy density and pressure are given by In the rest frame we also have Assuming χ to be pressureless in its rest frame, P RF χ = 0 , we thus find ρ RF χ = 2(V χ + V int ) and finally ρ χ − 3P χ = ρ RF χ = ∂V int /∂s ∂ log m χ (s)/∂s , (A.9) where the second equality can be easily verified for the class of potentials V χ , V int considered in this work, which are both quadratic in χ. In the pressureless limit, the above equation reduces to Eq. (2.14). On the background, the time component of the first equation in (A.4) takes the form which making use of the above results we rewrite aṡ ρ χ + 3H(ρ χ + P χ ) = (ρ χ − 3P χ ) ∂ log m χ (s) ∂sṡ . (A.11) This reduces to Eq. (3.1) in the pressureless limit.
Equation (A.11) can also be derived from the particle viewpoint, by integrating the Vlasov equation (2.6) in d 3 p E/(2π) 3 , performing integration by parts, and recalling the definitions ρ χ = d 3 pEf χ /(2π) 3 and P χ = d 3 p p 2 f χ /[(2π) 3 3E] . We note that there is yet another way to obtain the above equation from the particle perspective [22], by regarding the DM as a collection of point particles with number density n χ related to energy density and pressure by ρ χ = γm χ (s)n χ and P χ = γv 2 m χ (s)n χ /3. Therefore ρ χ − 3P χ = √ 1 − v 2 m χ (s)n χ , which manifests the expected decoupling in the ultra-relativistic limit v → 1. The connection with the field perspective is made through the relations n χ = γm χ (s)χ 2 for real scalar DM, n χ = 2γm χ (s)χ * χ for complex scalar DM, and n χ = γχχ for Dirac fermion DM.
The field approach can also be extended to first order, by making use of the relation between the perturbed fluid and field variables. In Newtonian gauge, which should be contrasted with the first of Eqs. (4.8). In addition the fluid velocity is given by (ρ χ + P χ )v i χ = −χ∇ i δχ/a. By employing these relations and the field equations of motion at the background and first-order levels, in the pressureless limit the continuity and Euler equations (4.1) and (4.2) are straightforwardly obtained.

(B.5)
As is well known, the above equations are obtained from those in NG by setting The Einstein's equations take the standard form [43]. It may be useful to recall the relation between our definitions of the NG potentials and those in Ref. [43]: Ψ here = ψ MB and Φ here = − ϕ MB .

B.1 Initial conditions in synchronous gauge
The adiabatic initial conditions in SG for all fluids except χ, s, as well as for the gravitational potentials, are found using standard methods [43]. By definition we have θ c = 0 and, retaining only the leading terms in the kτ ≪ 1 expansion, with C a dimensionless normalization constant and R ν ≡ρ ν /(ρ γ +ρ ν ). The initial condition for δ χ is found to satisfy adiabaticity, namely δ χ ≃ 3 4 δ γ up to corrections strongly suppressed bys ′ τ ≪ 1, as obtained from the RD solution in Eq. Thus for α na = 0 we find δ s = 3δ γ /2 in SG, which differs from δ s = δ γ /2 found at leading order in NG. There is, however, no inconsistency, as we now show by going to the next order. The SG potentials read h = C(kτ ) 2 + a h (kτ ) 4 , η = 2C − C 6 5 + 4R ν 15 + 4R ν (kτ ) 2 + a η (kτ ) 4 , (B.10) where a h , a η are constants whose explicit expressions are not needed for our purposes. The relative entropy perturbation between the scalar field and the photons in the two gauges is then, upon application of the background relation Eq. (3.9), neglecting O(k 4 τ 4 ) terms. Thus S sγ vanishes only at leading order, whereas it is nonzero (and gauge invariant as it must be) at next-to-leading order. A similar observation was made in Ref. [110].
Turning to the SG velocity potentials, for s we find Finally, the initial condition for the χ velocity potential is obtained from the Euler equation (B.4), ∂ log m χ (s) ∂s k(kτ ) . (B.14) Thus, for α na ̸ = 0 we find θ χ ∼ O(k 2 τ ), potentially much larger than for the photons and baryons which have θ ∼ O(k 4 τ 3 ).

B.2 Initial conditions in Newtonian gauge, summary
The adiabatic initial conditions in NG for all fluids except χ, s, and for the gravitational potentials, are at leading order in the kτ ≪ 1 expansion [43]  where the adiabatic sound speed squared is defined as This should not be confused with the sound speed squared c 2 sφ ≡ (δP s /δρ s ) RF , which always equals 1 for a scalar field. Instead, the adiabatic sound speed satisfies the gaugeinvariant relation The equation for the s velocity perturbation is found to be The second term in the parenthesis of the friction term contains the doom factor defined in Eq. (3.9), and we now see the reason for this terminology. If the doom factor is large and positive, there exists, even in the absence of sources, a solution to the Euler equation which goes as τ n , where n is proportional to the doom factor. In this scenario, perturbations become unstable at very early times [25,26]. For the case of Yukawa interactions discussed in this work, which rests on a microscopic Lagrangian description, the doom factor is negative and therefore does not lead to any runaway instabilities of the perturbations (in fact, the θ s term in Eq. (C.4) vanishes during early RD). We believe the same conclusion applies to any physical model of long range forces in the dark sector. We also rewrite the equations for χ in the same variables: Eq. (4.1) as and Eq. (4.2) as It is tedious but straightforward to check that these equations agree with Refs. [25,111], once the appropriate identifications are made. The SG equations are simply obtained by performing the replacements (B.6) in the above expressions.