Introduction
The recent DESI data analysis [1 , 2 ] provided tantalizing hints for dynamical dark energy, i.e. dark energy with a variable equation of state (EoS). The equation of state alone, however, does not capture the entire phenomenology of dark energy, and models based on scalar fields, somewhat analogous to inflationary scalar fields, have been proposed and examined [3 , 4 , 5 , 6 , 7 , 8 , 9 , 10 , 11 , 12 , 13 , 14 ] . In order to mimic the DESI best fit evolving EoS w ( z ) w(z) w ( z ) , however, these models need to reproduce the so-called “phantom crossing”, i.e. a transition from w < − 1 w < -1 w < − 1 to w > − 1 w > -1 w > − 1 at some recent epoch. This behavior cannot be obtained in ordinary scalar field models without at the same time creating instabilities [15 ] . Interestingly, non-minimally coupled scalar field models that introduce an interaction of the scalar field with matter seem to provide a good fit to current cosmological data while naturally producing an effective phantom crossing without introducing ghost or gradient instabilities [16 , 17 ] . Such models, however, automatically generate a fifth force acting on baryonic matter, that is tightly constrained by laboratory experiments and Solar System tests [18 , 19 ] (except if new ingredients are added to screen the fifth force on small scales [20 , 21 , 22 ] ). An elegant solution to evade those constraints while still solving the “phantom crossing” problem, is to couple only dark matter to the scalar field, as shown in [23 , 24 , 25 , 26 , 27 ] .
Perhaps the simplest and one of the most widely studied models of this type is a dark matter-dark energy interaction, namely the coupled dark energy model. In its basic incarnation, this model depends on a single constant coupling Q d m Q_{\mathrm{dm}} Q dm and on a potential, for instance a power law or an exponential. At the perturbation level, this model introduces, on sub-horizon scales, a Yukawa correction to the effective gravitational potential [28 , 29 ] , which for a point source of mass M M M has the following form in real space
Ψ ^ ( r ) = − G M r ( 1 + 2 Q d m 2 e − m r ) . (1) \hat{\Psi}(r)=-\frac{GM}{r}\left(1+2Q_{\mathrm{dm}}^{2}e^{-mr}\right). \tag*{(1)} Ψ ^ ( r ) = − r GM ( 1 + 2 Q dm 2 e − m r ) . (1)
Here, M M M and the scalar field mass m m m generally depend on the scalar field background value ϕ \phi ϕ , which, in turn, obeys a Klein-Gordon equation. The value of Q d m Q_{\mathrm{dm}} Q dm has been constrained in the literature several times, and a peak away from zero in its marginalized posterior was already found in some analyses before DESI [30 , 31 , 32 , 33 ] . Analyses including DESI data have also found hints of a nonzero coupling, although their significance depends on the datasets and assumptions adopted [24 , 25 , 34 , 35 , 36 , 37 ] . If this evidence is confirmed by additional data, e.g. Euclid, the dark-dark interaction could become a new standard cosmological model.
The coupled dark energy model contains a startling prediction: that the weak equivalence principle is broken at cosmological scales. That is, baryons and dark matter fall differently in a gravitational potential: while baryons remain uncoupled and feel standard gravity, following the geodesics of the metric, dark matter feels an extra interaction mediated by the scalar field and therefore deviates from the geodesics of the metric g μ ν g_{\mu\nu} g μν . Among the various effects induced by this extra interaction, the dark matter particle mass depends on the field ϕ \phi ϕ and therefore varies with time, cf. Refs [24 , 38 ] (and in principle also with space, but this would occur only beyond linear theory).
Recently, a direct way to test the EP in large scale structure has been identified [39 , 40 , 41 , 42 ] . The main idea is to employ the gravitational redshift effect to measure the gravitational potential and compare this with the effective potential Ψ ~ \tilde{\Psi} Ψ ~ mentioned above. This requires carefully incorporating horizon-scale relativistic corrections in the galaxy number count fluctuations [43 , 44 , 45 , 46 ] , which can be extracted from the two-point correlation function or the power spectrum by cross-correlating two (or more) sub-populations of galaxies with different biases [47 , 48 , 49 , 50 , 51 , 52 ] . This technique, multi-tracer with relativistic corrections, has been exploited in Ref. [41 ] to provide the first forecasts of EP-violation in present and future surveys, in a way that is independent of (i.e., marginalized over) the power spectrum shape, the background cosmology, and the perturbation growth rate.
This paper is devoted to forecasting the performance of this cosmological EP test in the context of coupled dark energy (CDE) and its generalization to doubly coupled dark energy (DCDE). In the latter model, the scalar field is conformally coupled to dark matter and nonminimally coupled to gravity. We consider the observable E P E_P E P identified in Ref. [41 ] , which provides a null test of the EP and can be measured in a model-independent way from galaxy power spectra. We express this quantity in the CDE and DCDE models, obtaining a general form that helps understand the roles of the various sectors of the theory’s Lagrangian. The model reduces to CDE when the nonminimal gravitational coupling is switched off. We then perform Fisher forecasts on these couplings for DESI, Euclid and SKA2 through a projection from E P E_P E P and additional background quantities. Our constraints on these couplings are model-independent in the sense that the underlying analysis leaves the shape of the power spectrum free to vary in several k k k -bins, and the galaxy bias also free to vary in redshift bins. This preserves the robustness of the analysis and allows us to estimate how well present and future surveys can constrain these couplings and E P E_P E P with galaxy multi-tracing and relativistic corrections, and hence assess the amount of EP-violation in our Universe.
Evolution equations in a doubly coupled dark energy model
In this section, we derive the evolution equations for the DCDE model, where dark matter and the metric are conformally coupled to a scalar field ϕ \phi ϕ . As mentioned above, a violation of the weak equivalence principle is already present when the scalar field couples only to dark matter, and it does not require the second coupling to the metric. However, the advantage of studying this more general doubly coupled model is that it provides insight about the impact of each coupling on cosmological observables and on their interplay regarding our test of the equivalence principle. In practice, the double coupling is introduced via two conformal factors f 1 f_1 f 1 and f 2 f_2 f 2 . To construct the simplest model of this kind with constant couplings, we choose appropriate functions f 1 , f 2 f_1, f_2 f 1 , f 2 . More general field-dependent couplings are of course possible, cf. Refs. [28 , 53 ] . The corresponding action of a scalar-tensor theory, where a dynamical scalar field ϕ \phi ϕ is nonminimally coupled to both gravity (cf. Ref. [54 ] ) and the dark matter sector (cf. Ref. [55 ] ), is given by:
S = ∫ d 4 x − g [ 1 2 κ f 1 ( ϕ ) R + L ϕ ( g μ ν ) + L b ( g μ ν , ψ b ) + L d m ( f 2 ( ϕ ) g μ ν , ψ D M ) ] , (2) S = \int d^4x\sqrt{-g}\left[\frac{1}{2\kappa}f_1(\phi)R + \mathcal{L}_{\phi}(g_{\mu\nu}) + \mathcal{L}_{b}(g_{\mu\nu},\psi_b) + \mathcal{L}_{\mathrm{dm}}(f_2(\phi)g_{\mu\nu},\psi_{\mathrm{DM}})\right], \tag*{(2)} S = ∫ d 4 x − g [ 2 κ 1 f 1 ( ϕ ) R + L ϕ ( g μν ) + L b ( g μν , ψ b ) + L dm ( f 2 ( ϕ ) g μν , ψ DM ) ] , (2)
where ψ b \psi_b ψ b and ψ d m \psi_{\mathrm{dm}} ψ dm are the fields of baryon and dark matter particles considered to be fermions, κ = 8 π G \kappa= 8\pi G κ = 8 π G is the inverse of the squared, reduced Planck mass, and
L ϕ ( g μ ν ) = − 1 2 ( 1 − 6 Q 2 ) f 1 ( ϕ ) ∇ μ ϕ ∇ μ ϕ − V ( ϕ ) , (3) \mathcal{L}_{\phi}(g_{\mu\nu}) = -\frac{1}{2}(1-6Q^2)f_1(\phi)\nabla_{\mu}\phi\nabla^{\mu}\phi- V(\phi), \tag*{(3)} L ϕ ( g μν ) = − 2 1 ( 1 − 6 Q 2 ) f 1 ( ϕ ) ∇ μ ϕ ∇ μ ϕ − V ( ϕ ) , (3)
where Q Q Q is a constant. The form of Eq. (3) is chosen so that Q Q Q will appear as the strength of the fifth-force induced by the coupling f 1 f_1 f 1 , as will be shown later. Specifically, f 1 f_1 f 1 induces a nonminimal coupling to the Ricci scalar, thereby modifying the gravitational sector, while f 2 f_2 f 2 introduces a direct coupling to dark matter through a conformal rescaling of the metric. While Q Q Q represents a universal modification of gravity, the equivalence principle is explicitly violated only when f 2 ≠ c o n s t f_2 \ne\mathrm{const} f 2 = const .
Varying the action (2) with respect to the metric g μ ν g_{\mu\nu} g μν and the scalar field ϕ \phi ϕ , one obtains the modified Einstein field equations and the modified Klein-Gordon equation:
f 1 G μ ν + g μ ν □ f 1 − ∇ μ ∇ ν f 1 = κ ( T μ ν d m + T μ ν b + T μ ν ϕ ) , (4) f_1G_{\mu\nu} + g_{\mu\nu}\Box f_1 - \nabla_{\mu}\nabla_{\nu}f_1 = \kappa\left(T^{\mathrm{dm}}_{\mu\nu} + T^{b}_{\mu\nu} + T^{\phi}_{\mu\nu}\right), \tag*{(4)} f 1 G μν + g μν □ f 1 − ∇ μ ∇ ν f 1 = κ ( T μν dm + T μν b + T μν ϕ ) , (4)
Figure 1. Schematic illustration of interactions between the Ricci scalar R R R , dark matter/baryon Lagrangian L d m / L b L_{\rm dm}/L_{\rm b} L dm / L b , and the scalar field ϕ \phi ϕ . The dashed arrow indicates minimal coupling via gravity, while the solid arrows represent direct couplings through the functions Q ( ϕ ) Q(\phi) Q ( ϕ ) and Q d m ( ϕ ) Q_{\rm dm}(\phi) Q dm ( ϕ ) .
□ ϕ + 1 2 f 1 , ϕ f 1 ∇ μ ϕ ∇ μ ϕ = 1 ( 1 − 6 Q 2 ) f 1 ( V , ϕ − 1 2 κ f 1 , ϕ R − 1 2 f 2 , ϕ f 2 T d m ) (5) \Box\phi+\frac{1}{2}\frac{f_{1,\phi}}{f_1}\nabla_\mu\phi\nabla^\mu\phi=\frac{1}{(1-6Q^2)f_1}\left(V_{,\phi}-\frac{1}{2\kappa}f_{1,\phi}R-\frac{1}{2}\frac{f_{2,\phi}}{f_2}T^{\rm dm}\right) \tag*{(5)} □ ϕ + 2 1 f 1 f 1 , ϕ ∇ μ ϕ ∇ μ ϕ = ( 1 − 6 Q 2 ) f 1 1 ( V , ϕ − 2 κ 1 f 1 , ϕ R − 2 1 f 2 f 2 , ϕ T dm ) (5)
where G μ ν = R μ ν − 1 2 R g μ ν G_{\mu\nu}=R_{\mu\nu}-\frac{1}{2}Rg_{\mu\nu} G μν = R μν − 2 1 R g μν ; □ ϕ ≡ ∇ μ ∇ μ ϕ \Box\phi\equiv\nabla_\mu\nabla^\mu\phi □ ϕ ≡ ∇ μ ∇ μ ϕ is the d’Alembertian operator acting on ϕ \phi ϕ ; the superscripts { b , d m } \{\mathrm{b},\mathrm{dm}\} { b , dm } denote the baryons and the dark matter components, respectively; T d m T^{\rm dm} T dm is the trace of the dark matter energy-momentum tensor; and f i , ϕ = d f i / d ϕ f_{i,\phi}=df_i/d\phi f i , ϕ = d f i / d ϕ ; V , ϕ = d V / d ϕ V_{,\phi}=dV/d\phi V , ϕ = d V / d ϕ .
We consider a flat, linearly perturbed Friedmann-Lemaître-Robertson-Walker (FLRW) metric in the conformal Newtonian gauge, where the two scalar degrees of freedom are characterized by the gravitational potentials Ψ \Psi Ψ and Φ \Phi Φ ,
d s 2 = a 2 [ − ( 1 + 2 Ψ ) d η 2 + ( 1 − 2 Φ ) d x i d x i ] (6) ds^2=a^2\left[-(1+2\Psi)d\eta^2+(1-2\Phi)dx^i dx_i\right] \tag*{(6)} d s 2 = a 2 [ − ( 1 + 2Ψ ) d η 2 + ( 1 − 2Φ ) d x i d x i ] (6)
where a a a is the scale factor and η \eta η is the conformal time. Substituting the metric Eq. (6) and assuming pressureless dark matter and baryons, one can obtain the Friedmann equation,
3 H 2 f 1 = κ [ ρ d m + ρ b + 1 2 ( 1 − 6 Q 2 ) f 1 H 2 ϕ ′ 2 + V ( ϕ ) ] − 3 H 2 f 1 , ϕ ϕ ′ (7) 3H^2f_1=\kappa\left[\rho_{\rm dm}+\rho_{\rm b}+\frac{1}{2}(1-6Q^2)f_1H^2\phi'^2+V(\phi)\right]-3H^2f_{1,\phi}\phi' \tag*{(7)} 3 H 2 f 1 = κ [ ρ dm + ρ b + 2 1 ( 1 − 6 Q 2 ) f 1 H 2 ϕ ′2 + V ( ϕ ) ] − 3 H 2 f 1 , ϕ ϕ ′ (7)
where a prime denotes differentiation with respect to the e e e -folding time N = ln a N=\ln a N = ln a . Note that the background energy density associated with the scalar field is given by ρ ϕ = 1 2 ( 1 − 6 Q 2 ) f 1 ( ϕ ) H 2 ϕ ′ 2 + V ( ϕ ) \rho_\phi=\frac{1}{2}(1-6Q^2)f_1(\phi)H^2\phi'^2+V(\phi) ρ ϕ = 2 1 ( 1 − 6 Q 2 ) f 1 ( ϕ ) H 2 ϕ ′2 + V ( ϕ ) , and the “effective” dark energy density reads
ρ D E , e f f = 1 2 ( 1 − 6 Q 2 ) f 1 ( ϕ ) H 2 ϕ ′ 2 + V ( ϕ ) − 3 H 2 f 1 , ϕ ϕ ′ κ (8) \rho_{\rm DE,eff}=\frac{1}{2}(1-6Q^2)f_1(\phi)H^2\phi'^2+V(\phi)-\frac{3H^2f_{1,\phi}\phi'}{\kappa} \tag*{(8)} ρ DE , eff = 2 1 ( 1 − 6 Q 2 ) f 1 ( ϕ ) H 2 ϕ ′2 + V ( ϕ ) − κ 3 H 2 f 1 , ϕ ϕ ′ (8)
such that
3 H 2 f 1 = κ ( ρ d m + ρ b + ρ D E , e f f ) (9) 3H^2f_1=\kappa\left(\rho_{\rm dm}+\rho_{\rm b}+\rho_{\rm DE,eff}\right) \tag*{(9)} 3 H 2 f 1 = κ ( ρ dm + ρ b + ρ DE , eff ) (9)
and one can define
Ω ~ d m = κ ρ d m 3 H 2 f 1 ; Ω ~ b = κ ρ b 3 H 2 f 1 ; Ω ~ D E , e f f = κ ρ D E , e f f 3 H 2 f 1 , (10) \widetilde{\Omega}_{\rm dm}=\frac{\kappa\rho_{\rm dm}}{3H^2f_1};\quad\widetilde{\Omega}_{\rm b}=\frac{\kappa\rho_{\rm b}}{3H^2f_1};\quad\widetilde{\Omega}_{\rm DE,eff}=\frac{\kappa\rho_{\rm DE,eff}}{3H^2f_1}, \tag*{(10)} Ω dm = 3 H 2 f 1 κ ρ dm ; Ω b = 3 H 2 f 1 κ ρ b ; Ω DE , eff = 3 H 2 f 1 κ ρ DE , eff , (10)
to ensure that Ω ~ d m + Ω ~ b + Ω ~ D E , e f f = 1 \widetilde{\Omega}_{\rm dm}+\widetilde{\Omega}_{\rm b}+\widetilde{\Omega}_{\rm DE,eff}=1 Ω dm + Ω b + Ω DE , eff = 1 . See also Ref. [56 ] in which these density fractions are defined in a similar way. The background Klein-Gordon equation is given by
ϕ ′ ′ + ( 3 + H ′ H ) ϕ ′ + 1 2 f 1 , ϕ f 1 ϕ ′ 2 = − 1 ( 1 − 6 Q 2 ) f 1 H 2 [ V , ϕ − 3 f 1 , ϕ κ H 2 ( 2 + H ′ H ) + f 2 , ϕ 2 f 2 ρ d m ] (11) \phi''+\left(3+\frac{H'}{H}\right)\phi'+\frac{1}{2}\frac{f_{1,\phi}}{f_1}\phi'^2=-\frac{1}{(1-6Q^2)f_1H^2}\left[V_{,\phi}-\frac{3f_{1,\phi}}{\kappa}H^2\left(2+\frac{H'}{H}\right)+\frac{f_{2,\phi}}{2f_2}\rho_{\rm dm}\right] \tag*{(11)} ϕ ′′ + ( 3 + H H ′ ) ϕ ′ + 2 1 f 1 f 1 , ϕ ϕ ′2 = − ( 1 − 6 Q 2 ) f 1 H 2 1 [ V , ϕ − κ 3 f 1 , ϕ H 2 ( 2 + H H ′ ) + 2 f 2 f 2 , ϕ ρ dm ] (11)
Defining Q d m ( ϕ ) ≡ f 2 , ϕ 2 f 2 Q_{\rm dm}(\phi)\equiv\frac{f_{2,\phi}}{2f_2} Q dm ( ϕ ) ≡ 2 f 2 f 2 , ϕ , the conservation equations for baryons, dark matter, and the effective dark energy component can also be derived from the action (2) :
∇ ν T μ ν b = 0 (12) \nabla^\nu T^{\rm b}_{\mu\nu}=0 \tag*{(12)} ∇ ν T μν b = 0 (12)
∇ ν T μ ν d m = Q d m T d m ∇ μ ϕ (13) \nabla^\nu T^{\rm dm}_{\mu\nu}=Q_{\rm dm}T^{\rm dm}\nabla_\mu\phi \tag*{(13)} ∇ ν T μν dm = Q dm T dm ∇ μ ϕ (13)
∇ ν T μ ν D E , e f f = − Q d m T d m ∇ μ ϕ + 1 κ G μ ν ∇ ν f 1 (14) \nabla^\nu T^{\rm DE,eff}_{\mu\nu}=-Q_{\rm dm}T^{\rm dm}\nabla_\mu\phi+\frac{1}{\kappa}G_{\mu\nu}\nabla^\nu f_1 \tag*{(14)} ∇ ν T μν DE , eff = − Q dm T dm ∇ μ ϕ + κ 1 G μν ∇ ν f 1 (14)
The time components of the above conservation equations at the background level then lead to
ρ b ′ + 3 ρ b = 0 , (15) \rho'_b+3\rho_b=0, \tag*{(15)} ρ b ′ + 3 ρ b = 0 , (15)
ρ d m ′ + 3 ρ d m = Q d m ρ d m ϕ ′ , (16) \rho'_{\rm dm}+3\rho_{\rm dm}=Q_{\rm dm}\rho_{\rm dm}\phi', \tag*{(16)} ρ dm ′ + 3 ρ dm = Q dm ρ dm ϕ ′ , (16)
ρ D E , e f f ′ + 3 ( ρ D E , e f f + p D E , e f f ) = − Q d m ρ d m ϕ ′ + 3 H 2 κ f 1 , ϕ ϕ ′ . (17) \rho'_{\rm DE,eff}+3(\rho_{\rm DE,eff}+p_{\rm DE,eff})=-Q_{\rm dm}\rho_{\rm dm}\phi'+\frac{3\mathcal{H}^2}{\kappa}f_{1,\phi}\phi'. \tag*{(17)} ρ DE , eff ′ + 3 ( ρ DE , eff + p DE , eff ) = − Q dm ρ dm ϕ ′ + κ 3 H 2 f 1 , ϕ ϕ ′ . (17)
Solving Eq. (16) , we have
ρ d m ( N ) = m d m ( N ) n d m ( N ) = ρ d m , 0 e − 3 N e ∫ ϕ 0 ϕ Q d m ( ϕ ^ ) d ϕ ^ , (18) \rho_{\rm dm}(N)=m_{\rm dm}(N)n_{\rm dm}(N)=\rho_{{\rm dm},0}e^{-3N}e^{\int_{\phi_0}^{\phi}Q_{\rm dm}(\hat{\phi})d\hat{\phi}}, \tag*{(18)} ρ dm ( N ) = m dm ( N ) n dm ( N ) = ρ dm , 0 e − 3 N e ∫ ϕ 0 ϕ Q dm ( ϕ ^ ) d ϕ ^ , (18)
where a subscript ”0” denotes the quantities at the current epoch. Assuming a conserved total number of dark matter particles, the particle mass varies due to the presence of the coupling Q d m Q_{\rm dm} Q dm , and evolves as m d m ( N ) = m d m , 0 e ∫ ϕ 0 ϕ Q d m ( ϕ ^ ) d ϕ ^ m_{\rm dm}(N)=m_{{\rm dm},0}e^{\int_{\phi_0}^{\phi}Q_{\rm dm}(\hat{\phi})d\hat{\phi}} m dm ( N ) = m dm , 0 e ∫ ϕ 0 ϕ Q dm ( ϕ ^ ) d ϕ ^ .
At the perturbation level, the perturbed scalar field ϕ ( p ) \phi_{\rm(p)} ϕ ( p ) , the density ρ m , ( p ) \rho_{\rm m,(p)} ρ m , ( p ) and four-velocity u m , ( p ) μ u^\mu_{\rm m,(p)} u m , ( p ) μ of the matter contents, on the other hand, are expanded as
ϕ ( p ) = ϕ + φ ; ρ m , ( p ) = ρ m ( 1 + δ m ) ; u m , ( p ) μ = 1 a ( 1 − Ψ , v m i ) (19) \phi_{\rm(p)}=\phi+\varphi\ ;\quad\rho_{\rm m,(p)}=\rho_{\rm m}(1+\delta_{\rm m})\ ;\quad u^\mu_{\rm m,(p)}=\frac{1}{a}(1-\Psi,v^i_{\rm m}) \tag*{(19)} ϕ ( p ) = ϕ + φ ; ρ m , ( p ) = ρ m ( 1 + δ m ) ; u m , ( p ) μ = a 1 ( 1 − Ψ , v m i ) (19)
with a subscript (p) denoting the perturbed quantities, v i = d x i / d η v^i=dx^i/d\eta v i = d x i / d η the peculiar velocity defined using the comoving coordinates x i x^i x i , and m = { b , d m } m=\{\mathrm{b},\mathrm{dm}\} m = { b , dm } . Inserting Eq. (6) into the matter conservation equations (12) and (13) , the continuity and Euler equations in Fourier space3 , are given by
δ b ′ = − θ b + 3 Φ ′ , (20) \delta'_b=-\theta_b+3\Phi', \tag*{(20)} δ b ′ = − θ b + 3 Φ ′ , (20)
θ b ′ + θ b ( 1 + H ′ H ) − Ψ λ 2 = 0 , (21) \theta'_b+\theta_b\left(1+\frac{\mathcal{H}'}{\mathcal{H}}\right)-\frac{\Psi}{\lambda^2}=0, \tag*{(21)} θ b ′ + θ b ( 1 + H H ′ ) − λ 2 Ψ = 0 , (21)
δ d m ′ + θ d m − 3 Φ ′ = Q d m φ ′ + Q d m , ϕ ϕ ′ φ , (22) \delta'_{\rm dm}+\theta_{\rm dm}-3\Phi'=Q_{\rm dm}\varphi'+Q_{{\rm dm},\phi}\phi'\varphi, \tag*{(22)} δ dm ′ + θ dm − 3 Φ ′ = Q dm φ ′ + Q dm , ϕ ϕ ′ φ , (22)
θ d m ′ + θ d m ( 1 + H ′ H + Q d m ϕ ′ ) − Q d m φ + Ψ λ 2 = 0 , (23) \theta'_{\rm dm}+\theta_{\rm dm}\left(1+\frac{\mathcal{H}'}{\mathcal{H}}+Q_{\rm dm}\phi'\right)-\frac{Q_{\rm dm}\varphi+\Psi}{\lambda^2}=0, \tag*{(23)} θ dm ′ + θ dm ( 1 + H H ′ + Q dm ϕ ′ ) − λ 2 Q dm φ + Ψ = 0 , (23)
where θ m \theta_m θ m is the rescaled velocity divergence θ m = i k i v m i / H \theta_m=ik_iv^i_m/\mathcal{H} θ m = i k i v m i / H with wave vector k k k and H = a H \mathcal{H}=aH H = a H being the conformal Hubble function. We also introduce λ = H / k \lambda=\mathcal{H}/k λ = H / k to characterize the scale dependence and define Q d m , ϕ ≡ d Q d m / d ϕ Q_{{\rm dm},\phi}\equiv dQ_{\rm dm}/d\phi Q dm , ϕ ≡ d Q dm / d ϕ . On the other hand, from the perturbed Einstein equations we obtain
Φ = − λ 2 2 f 1 [ 2 V ^ Ψ + φ V ^ , ϕ − 3 φ f 1 , ϕ ϕ ϕ ′ + 3 ( δ b Ω b + δ d m Ω d m ) + Φ ′ ( f 1 + 3 ϕ ′ f 1 , ϕ ) − 3 f 1 , ϕ ( φ + φ ′ ) + 6 Ψ ( Ω b + Ω d m ) + κ ( 1 − 6 Q 2 ) ( f 1 ϕ ′ φ ′ + f 1 , ϕ ϕ ′ φ 2 2 ) ] + φ f 1 , ϕ 2 f 1 , (24) \Phi=-\frac{\lambda^2}{2f_1}\left[2\hat{V}\Psi+\varphi\hat{V}_{,\phi}-3\varphi f_{1,\phi\phi}\phi'+3(\delta_b\Omega_b+\delta_{\rm dm}\Omega_{\rm dm})+\Phi'(f_1+3\phi'f_{1,\phi})-3f_{1,\phi}(\varphi+\varphi')+6\Psi(\Omega_b+\Omega_{\rm dm})+\kappa(1-6Q^2)\left(f_1\phi'\varphi'+\frac{f_{1,\phi}\phi'\varphi^2}{2}\right)\right]+\varphi\frac{f_{1,\phi}}{2f_1}, \tag*{(24)} Φ = − 2 f 1 λ 2 [ 2 V ^ Ψ + φ V ^ , ϕ − 3 φ f 1 , ϕϕ ϕ ′ + 3 ( δ b Ω b + δ dm Ω dm ) + Φ ′ ( f 1 + 3 ϕ ′ f 1 , ϕ ) − 3 f 1 , ϕ ( φ + φ ′ ) + 6Ψ ( Ω b + Ω dm ) + κ ( 1 − 6 Q 2 ) ( f 1 ϕ ′ φ ′ + 2 f 1 , ϕ ϕ ′ φ 2 ) ] + φ 2 f 1 f 1 , ϕ , (24)
Ψ = Φ − φ f 1 , ϕ f 1 , (25) \Psi=\Phi-\varphi\frac{f_{1,\phi}}{f_1}, \tag*{(25)} Ψ = Φ − φ f 1 f 1 , ϕ , (25)
with V ^ = V a 2 / H 2 \hat{V}=Va^2/\mathcal{H}^2 V ^ = V a 2 / H 2 . From now on, we adopt the reduced Planck units where κ = c = 1 \kappa=c=1 κ = c = 1 . Also, to avoid confusion, notice that Ω m ( N ) = κ ρ m / 3 H 2 = f 1 Ω ~ d m \Omega_m(N)=\kappa\rho_m/3\mathcal{H}^2=f_1\widetilde{\Omega}_{\rm dm} Ω m ( N ) = κ ρ m /3 H 2 = f 1 Ω dm is a time-dependent quantity, not the present values. So far, we have allowed Q d m Q_{\rm dm} Q dm to vary with time. Hereafter we specialize to constant couplings Q = − f 1 , ϕ 2 f 1 Q=-\frac{f_{1,\phi}}{2f_1} Q = − 2 f 1 f 1 , ϕ and Q d m = f 2 , ϕ 2 f 2 Q_{\rm dm}=\frac{f_{2,\phi}}{2f_2} Q dm = 2 f 2 f 2 , ϕ , because in this case one obtains a Yukawa correction in Eq. (1) with a constant coupling. Then the general solutions are f 1 ( ϕ ) = f 1 , 0 e − 2 Q ( ϕ − ϕ 0 ) f_1(\phi)=f_{1,0}e^{-2Q(\phi-\phi_0)} f 1 ( ϕ ) = f 1 , 0 e − 2 Q ( ϕ − ϕ 0 ) , f 2 ( ϕ ) = f 2 , 0 e 2 Q d m ( ϕ − ϕ 0 ) f_2(\phi)=f_{2,0}e^{2Q_{\rm dm}(\phi-\phi_0)} f 2 ( ϕ ) = f 2 , 0 e 2 Q dm ( ϕ − ϕ 0 ) , with ϕ 0 \phi_0 ϕ 0 the present-day field value and f i , 0 ≡ f i ( ϕ 0 ) f_{i,0}\equiv f_i(\phi_0) f i , 0 ≡ f i ( ϕ 0 ) . We adopt the normalization f 1 , 0 = f 2 , 0 = 1 f_{1,0}=f_{2,0}=1 f 1 , 0 = f 2 , 0 = 1 . In addition, we can further set ϕ 0 = 0 \phi_0=0 ϕ 0 = 0 without loss of generality. With these conventions, the coupling functions become f 1 = e − 2 Q ϕ f_1=e^{-2Q\phi} f 1 = e − 2 Qϕ and f 2 = e 2 Q d m ϕ f_2=e^{2Q_{\rm dm}\phi} f 2 = e 2 Q dm ϕ .
On sub-horizon scales (λ ≪ 1 \lambda\ll1 λ ≪ 1 ), we adopt the quasi-static assumption (QSA), as described in [57 ] , by implying that spatial derivatives dominate over the time derivatives of slowly evolving perturbations, i.e. X ∼ X ′ ∼ X ′ ′ X\sim X'\sim X'' X ∼ X ′ ∼ X ′′ for
f ~ ( k , t ) = ∫ d 3 x f ( x , t ) e − i k ⋅ x ; f ( x , t ) = ∫ d 3 k ( 2 π ) 3 f ~ ( k , t ) e i k ⋅ x . \tilde{f}(\mathbf{k},t)=\int d^3x f(\mathbf{x},t)e^{-i\mathbf{k}\cdot\mathbf{x}};\quad f(\mathbf{x},t)=\int\frac{d^3k}{(2\pi)^3}\tilde{f}(\mathbf{k},t)e^{i\mathbf{k}\cdot\mathbf{x}}. f ~ ( k , t ) = ∫ d 3 x f ( x , t ) e − i k ⋅ x ; f ( x , t ) = ∫ ( 2 π ) 3 d 3 k f ~ ( k , t ) e i k ⋅ x .
all perturbation quantities X ∈ { δ , θ , Ψ , Φ } X \in\{\delta,\theta,\Psi,\Phi\} X ∈ { δ , θ , Ψ , Φ } . In this case X ′ X' X ′ , X ′ ′ ∼ X ≪ λ − 2 X X'' \sim X \ll\lambda^{-2}X X ′′ ∼ X ≪ λ − 2 X .4 Consequently, the perturbed Klein-Gordon equation reads,
( 1 − 6 Q 2 + a 2 m ϕ 2 k 2 ) φ = − 3 λ 2 Q d m Ω ~ d m δ d m − Q ( 2 Ψ − 4 Φ ) (26) \left(1-6Q^2+\frac{a^2m_\phi^2}{k^2}\right)\varphi=-3\lambda^2Q_{\rm dm}\widetilde{\Omega}_{\rm dm}\delta_{\rm dm}-Q(2\Psi-4\Phi) \tag*{(26)} ( 1 − 6 Q 2 + k 2 a 2 m ϕ 2 ) φ = − 3 λ 2 Q dm Ω dm δ dm − Q ( 2Ψ − 4Φ ) (26)
where m ϕ 2 = ( V , ϕ ϕ + 2 Q V , ϕ + 2 Q Q d m ρ d m ) / f 1 ( ϕ ) m_\phi^2=(V_{,\phi\phi}+2QV_{,\phi}+2QQ_{\rm dm}\rho_{\rm dm})/f_1(\phi) m ϕ 2 = ( V , ϕϕ + 2 Q V , ϕ + 2 Q Q dm ρ dm ) / f 1 ( ϕ ) . Moreover, under the QSA Eq. (24) reduces to
Φ = − 3 λ 2 2 ( δ b Ω ~ b + δ d m Ω ~ d m ) − Q φ (27) \Phi=-\frac{3\lambda^2}{2}(\delta_b\widetilde{\Omega}_b+\delta_{\rm dm}\widetilde{\Omega}_{\rm dm})-Q\varphi \tag*{(27)} Φ = − 2 3 λ 2 ( δ b Ω b + δ dm Ω dm ) − Qφ (27)
and Eq. (26) can be brought into the convenient form,
φ = − 3 λ 2 Y ( k ) [ ( Q + Q d m ) Ω ~ d m δ d m + Q Ω ~ b δ b ] ; Y ( k ) = k 2 k 2 + a 2 m ϕ 2 . (28) \varphi=-3\lambda^2Y(k)\left[(Q+Q_{\rm dm})\widetilde{\Omega}_{\rm dm}\delta_{\rm dm}+Q\widetilde{\Omega}_b\delta_b\right]\ ;\quad Y(k)=\frac{k^2}{k^2+a^2m_\phi^2}. \tag*{(28)} φ = − 3 λ 2 Y ( k ) [ ( Q + Q dm ) Ω dm δ dm + Q Ω b δ b ] ; Y ( k ) = k 2 + a 2 m ϕ 2 k 2 . (28)
Thus, Eq. (25) can be re-written for the scalar field Ψ \Psi Ψ :
Ψ = − 3 λ 2 2 [ ( 2 Q Y ( k ) ( Q + Q d m ) + 1 ) Ω ~ d m δ d m + ( 2 Q 2 Y ( k ) + 1 ) Ω ~ b δ b ] . (29) \Psi=-\frac{3\lambda^2}{2}\left[(2QY(k)(Q+Q_{\rm dm})+1)\widetilde{\Omega}_{\rm dm}\delta_{\rm dm}+(2Q^2Y(k)+1)\widetilde{\Omega}_b\delta_b\right]. \tag*{(29)} Ψ = − 2 3 λ 2 [ ( 2 Q Y ( k ) ( Q + Q dm ) + 1 ) Ω dm δ dm + ( 2 Q 2 Y ( k ) + 1 ) Ω b δ b ] . (29)
These are the Poisson equations for the perturbation quantities φ \varphi φ and Ψ \Psi Ψ .
To determine the growth of matter perturbations at deep subhorizon scales, the continuity equations (20) and (22) can be simplified to δ b ′ = − θ b \delta'_b=-\theta_b δ b ′ = − θ b and δ d m ′ = − θ d m \delta'_{\rm dm}=-\theta_{\rm dm} δ dm ′ = − θ dm , respectively, because from the Poisson equations one sees that Ψ , φ ∼ λ 2 δ m \Psi,\varphi\sim\lambda^2\delta_m Ψ , φ ∼ λ 2 δ m . Taking the derivative of (20) and (22) and replacing θ b , θ d m \theta_b,\theta_{\rm dm} θ b , θ dm and their derivatives, the growth equations for baryons and dark matter takes the following form,
δ b ′ ′ + δ b ′ ( 1 + H ′ H ) + Ψ λ 2 = 0 (30) \delta_b''+\delta_b'\left(1+\frac{\mathcal{H}'}{\mathcal{H}}\right)+\frac{\Psi}{\lambda^2}=0 \tag*{(30)} δ b ′′ + δ b ′ ( 1 + H H ′ ) + λ 2 Ψ = 0 (30)
δ d m ′ ′ + δ d m ′ ( 1 + H ′ H + Q d m ϕ ′ ) + Ψ + Q d m φ λ 2 = 0. (31) \delta_{\rm dm}''+\delta_{\rm dm}'\left(1+\frac{\mathcal{H}'}{\mathcal{H}}+Q_{\rm dm}\phi'\right)+\frac{\Psi+Q_{\rm dm}\varphi}{\lambda^2}=0. \tag*{(31)} δ dm ′′ + δ dm ′ ( 1 + H H ′ + Q dm ϕ ′ ) + λ 2 Ψ + Q dm φ = 0. (31)
Plugging in the gravitational potentials we have
δ b ′ ′ + ( 1 + H ′ H ) δ b ′ − 3 2 [ ( 1 + 2 Q Y ( Q + Q d m ) ) Ω ~ d m δ d m + ( 1 + 2 Q 2 Y ) Ω ~ b δ b ] = 0 (32) \delta_b''+\left(1+\frac{\mathcal{H}'}{\mathcal{H}}\right)\delta_b'-\frac{3}{2}\left[(1+2QY(Q+Q_{\rm dm}))\widetilde{\Omega}_{\rm dm}\delta_{\rm dm}+(1+2Q^2Y)\widetilde{\Omega}_b\delta_b\right]=0 \tag*{(32)} δ b ′′ + ( 1 + H H ′ ) δ b ′ − 2 3 [ ( 1 + 2 Q Y ( Q + Q dm )) Ω dm δ dm + ( 1 + 2 Q 2 Y ) Ω b δ b ] = 0 (32)
δ d m ′ ′ + ( 1 + H ′ H + Q d m ϕ ′ ) δ d m ′ − 3 2 { [ 1 + 2 Y ( Q + Q d m ) 2 ] Ω ~ d m δ d m + [ 1 + 2 Q Y ( Q + Q d m ) ] Ω ~ b δ b } = 0. (33) \delta_{\rm dm}''+\left(1+\frac{\mathcal{H}'}{\mathcal{H}}+Q_{\rm dm}\phi'\right)\delta_{\rm dm}'-\frac{3}{2}\left\{[1+2Y(Q+Q_{\rm dm})^2]\widetilde{\Omega}_{\rm dm}\delta_{\rm dm}+[1+2QY(Q+Q_{\rm dm})]\widetilde{\Omega}_b\delta_b\right\}=0. \tag*{(33)} δ dm ′′ + ( 1 + H H ′ + Q dm ϕ ′ ) δ dm ′ − 2 3 { [ 1 + 2 Y ( Q + Q dm ) 2 ] Ω dm δ dm + [ 1 + 2 Q Y ( Q + Q dm )] Ω b δ b } = 0. (33)
As expected, for Q = 0 Q=0 Q = 0 , these equations reduce to the growth equations in the CDE scenario, see e.g., Refs. [32 , 38 ] .
The effective potential acting on dark matter test particles mentioned in the introduction is simply Ψ ^ d m ≡ Ψ + Q d m φ \widehat{\Psi}_{\rm dm}\equiv\Psi+Q_{\rm dm}\varphi Ψ dm ≡ Ψ + Q dm φ , while for baryon test particles it is Ψ ^ b ≡ Ψ \widehat{\Psi}_b\equiv\Psi Ψ b ≡ Ψ . As anticipated, in real space the gravitational potential is modified by the term Y ( k ) Y(k) Y ( k ) and the Newtonian potential around a point-like dark matter particle acquires a Yukawa correction 2 ( Q + Q d m ) 2 e − m ( φ ) r 2(Q+Q_{\rm dm})^2e^{-m(\varphi)r} 2 ( Q + Q dm ) 2 e − m ( φ ) r , while the correction for a baryon particle is 2 Q 2 e − m ( φ ) r 2Q^2e^{-m(\varphi)r} 2 Q 2 e − m ( φ ) r . We present this explicitly in App. A. One can then think of dark matter particles as carrying a “charge” for the field interaction Q + Q d m Q+Q_{\rm dm} Q + Q dm , and of baryons as carrying a charge Q Q Q .
Testing the Weak Equivalence Principle
In this section, we derive a combination of quantities in the model described in Sec. II that can be employed to test the EP with galaxy clustering, as shown in Refs. [41 , 42 ] .
For a generic model, the Poisson equation for the gravitational potential Ψ \Psi Ψ takes the form,
Ψ = − 3 2 ∑ m μ G , m ( k , z ) λ 2 Ω m δ m (34) \Psi=-\frac{3}{2}\sum_m\mu_{G,m}(k,z)\lambda^2\Omega_m\delta_m \tag*{(34)} Ψ = − 2 3 m ∑ μ G , m ( k , z ) λ 2 Ω m δ m (34)
where the summation runs over the matter species m = { d m , b } \mathrm{m} = \{\mathrm{dm}, \mathrm{b}\} m = { dm , b } . Here, μ G , m \mu_{G,\mathrm{m}} μ G , m is the effective gravitational “constant” that could possibly depend on both time and scale, with μ G , m ≠ 1 \mu_{G,\mathrm{m}} \ne1 μ G , m = 1 quantifying departures from general relativity. The previously derived Poisson equation (29) can be used to identify that,
μ G , d m ( k , z ) = 1 f 1 [ Y ( k ) ( 2 Q 2 + 2 Q d m Q ) + 1 ] , (35) \mu_{G,\mathrm{dm}}(k,z) = \frac{1}{f_1}\left[Y(k)(2Q^2+2Q_{\mathrm{dm}}Q)+1\right], \tag*{(35)} μ G , dm ( k , z ) = f 1 1 [ Y ( k ) ( 2 Q 2 + 2 Q dm Q ) + 1 ] , (35)
μ G , b ( k , z ) = 1 f 1 [ Y ( k ) 2 Q 2 + 1 ] . (36) \mu_{G,\mathrm{b}}(k,z) = \frac{1}{f_1}\left[Y(k)2Q^2+1\right]. \tag*{(36)} μ G , b ( k , z ) = f 1 1 [ Y ( k ) 2 Q 2 + 1 ] . (36)
In Ref. [41 ] , it was shown that a parameter, denoted as E P E_P E P , can be identified and measured in galaxy clustering surveys when relativistic corrections (in particular, terms associated with the gravitational redshift) are taken into account (see details in Sec. IV). This parameter E P E_P E P measures the deviation from the EP between baryons and dark matter, and can be estimated in a model-independent way, i.e. without choosing the initial conditions of the Universe or a particular cosmological model, either at the background or at the perturbation level. As we now show, for this DCDE model, E P E_P E P can be explicitly constructed from the Poisson and the Euler equations for dark matter and baryons.
Proposed in [39 ] , the most general cosmological Euler equation for dark matter in the modified gravity scenario can be expressed by introducing two phenomenological free functions (see, e.g., Refs [41 , 42 , 58 ] ), that become nonzero when deviating from the standard Λ \Lambda Λ CDM scenario,
θ d m ′ + θ d m [ 1 + H ′ H + Θ d m ( k , z ) ] − [ 1 + Γ d m ( k , z ) ] λ − 2 Ψ = 0. (37) \theta'_{\mathrm{dm}}+\theta_{\mathrm{dm}}\left[1+\frac{\mathcal{H}'}{\mathcal{H}}+\Theta_{\mathrm{dm}}(k,z)\right]-\left[1+\Gamma_{\mathrm{dm}}(k,z)\right]\lambda^{-2}\Psi=0. \tag*{(37)} θ dm ′ + θ dm [ 1 + H H ′ + Θ dm ( k , z ) ] − [ 1 + Γ dm ( k , z ) ] λ − 2 Ψ = 0. (37)
Here, the functions Θ d m , Γ d m \Theta_{\mathrm{dm}},\Gamma_{\mathrm{dm}} Θ dm , Γ dm correspond to an additional friction term and to the fifth force acting on dark matter, respectively. By comparing with Eq. (23) one can find that [28 ]
Θ d m = Q d m ϕ ′ ; Γ d m = Q d m ϕ Ψ . (38) \Theta_{\mathrm{dm}}=Q_{\mathrm{dm}}\phi';\qquad\Gamma_{\mathrm{dm}}=\frac{Q_{\mathrm{dm}}\phi}{\Psi}. \tag*{(38)} Θ dm = Q dm ϕ ′ ; Γ dm = Ψ Q dm ϕ . (38)
As expected, Θ \Theta Θ and Γ \Gamma Γ both vanish when Q d m → 0 Q_{\mathrm{dm}}\to0 Q dm → 0 , i.e. when the EP is satisfied. Substituting ϕ \phi ϕ and Ψ \Psi Ψ with Eqs. (28) and (29) , Γ d m \Gamma_{\mathrm{dm}} Γ dm becomes
Γ d m = 2 Q d m Y ( k ) ( Q + Q d m ) Ω d m + Q Ω b b [ 2 Q Y ( k ) ( Q + Q d m ) + 1 ] Ω d m + [ 2 Q 2 Y ( k ) + 1 ] Ω b b . (39) \Gamma_{\mathrm{dm}}=2Q_{\mathrm{dm}}Y(k)\frac{(Q+Q_{\mathrm{dm}})\Omega_{\mathrm{dm}}+Q\Omega_{\mathrm{b}}b}{\left[2QY(k)(Q+Q_{\mathrm{dm}})+1\right]\Omega_{\mathrm{dm}}+\left[2Q^2Y(k)+1\right]\Omega_{\mathrm{b}}b}. \tag*{(39)} Γ dm = 2 Q dm Y ( k ) [ 2 Q Y ( k ) ( Q + Q dm ) + 1 ] Ω dm + [ 2 Q 2 Y ( k ) + 1 ] Ω b b ( Q + Q dm ) Ω dm + Q Ω b b . (39)
with the time-dependent bias function b ( z ) b(z) b ( z ) defined by δ b = b ( z ) δ d m \delta_{\mathrm{b}}=b(z)\delta_{\mathrm{dm}} δ b = b ( z ) δ dm . On linear scales, it can be obtained by integrating Eqs. (30) and (31) . Notice that this bias is exclusively due to the EP violation and does not include fluidodynamical effects that are expected to be important at non-linear scales. Referring only to the dark matter component, E P E_P E P is defined in [41 ] in the following way,
E P = 1 + Θ d m − 3 Ω d m μ G , d m Γ d m 2 f g , (40) E_P=1+\Theta_{\mathrm{dm}}-\frac{3\Omega_{\mathrm{dm}}\mu_{G,\mathrm{dm}}\Gamma_{\mathrm{dm}}}{2f_g}, \tag*{(40)} E P = 1 + Θ dm − 2 f g 3 Ω dm μ G , dm Γ dm , (40)
where f g ≡ δ d m ′ / δ d m f_g\equiv\delta'_{\mathrm{dm}}/\delta_{\mathrm{dm}} f g ≡ δ dm ′ / δ dm is the dark matter growth rate. It is important to remark that neither Θ d m \Theta_{\mathrm{dm}} Θ dm nor Γ d m \Gamma_{\mathrm{dm}} Γ dm are directly observable, unless one first assume a specific cosmology, while E P E_P E P is a directly observable quantity, as will be shown in Sect. IV. Whenever E P ≠ 1 E_P\ne1 E P = 1 , the EP is violated, whether originating from non-standard gravitational interactions or non-minimal couplings in the dark sector. Combining Eqs. (35) , (38) and (39) , E P E_P E P for dark matter can be expressed in the most general form,
E P = 1 + Q d m ϕ ′ − 3 Ω d m Q d m Y ( k ) 2 f g f 1 [ 2 ( Q + Q d m ) Ω d m + 2 Q Ω b b ] [ 2 Q Y ( k ) ( Q + Q d m ) + 1 ] [ 2 Q Y ( k ) ( Q + Q d m ) + 1 ] Ω d m + [ 2 Q 2 Y ( k ) + 1 ] Ω b b , (41) E_P=1+Q_{\mathrm{dm}}\phi'-\frac{3\Omega_{\mathrm{dm}}Q_{\mathrm{dm}}Y(k)}{2f_gf_1}\frac{\left[2(Q+Q_{\mathrm{dm}})\Omega_{\mathrm{dm}}+2Q\Omega_{\mathrm{b}}b\right]\left[2QY(k)(Q+Q_{\mathrm{dm}})+1\right]}{\left[2QY(k)(Q+Q_{\mathrm{dm}})+1\right]\Omega_{\mathrm{dm}}+\left[2Q^2Y(k)+1\right]\Omega_{\mathrm{b}}b}, \tag*{(41)} E P = 1 + Q dm ϕ ′ − 2 f g f 1 3 Ω dm Q dm Y ( k ) [ 2 Q Y ( k ) ( Q + Q dm ) + 1 ] Ω dm + [ 2 Q 2 Y ( k ) + 1 ] Ω b b [ 2 ( Q + Q dm ) Ω dm + 2 Q Ω b b ] [ 2 Q Y ( k ) ( Q + Q dm ) + 1 ] , (41)
where we recall that f 1 = e − 2 Q ϕ f_1=e^{-2Q\phi} f 1 = e − 2 Qϕ . This is the main result of this section. Assuming baryons to be sub-dominant, E P E_P E P reduces to
E P = 1 + Q d m ϕ ′ − 3 Q d m Y ( k ) Ω d m f 1 Q + Q d m f g = 1 + Q d m ϕ ′ − 3 Q d m Ω ~ d m Y ( k ) Q + Q d m f g . (42) \begin{aligned}
E_P=1+Q_{\mathrm{dm}}\phi'-3Q_{\mathrm{dm}}Y(k)\frac{\Omega_{\mathrm{dm}}}{f_1}\frac{Q+Q_{\mathrm{dm}}}{f_g} \\
&=1+Q_{\mathrm{dm}}\phi'-3Q_{\mathrm{dm}}\widetilde{\Omega}_{\mathrm{dm}}Y(k)\frac{Q+Q_{\mathrm{dm}}}{f_g}. \tag*{(42)}
\end{aligned} E P = 1 + Q dm ϕ ′ − 3 Q dm Y ( k ) f 1 Ω dm f g Q + Q dm = 1 + Q dm ϕ ′ − 3 Q dm Ω dm Y ( k ) f g Q + Q dm . (42)
Recall that Ω ~ d m = Ω d m / f 1 \widetilde{\Omega}_{\rm dm}=\Omega_{\rm dm}/f_1 Ω dm = Ω dm / f 1 . As expected, E P → 1 E_P\to1 E P → 1 when Q d m → 0 Q_{\rm dm}\to0 Q dm → 0 . We introduce the dimensionless variables
x ( N ) = ϕ ′ 6 , y ( N ) = 1 H V 3 f 1 . (43) x(N)=\frac{\phi'}{\sqrt{6}},\qquad y(N)=\frac{1}{H}\sqrt{\frac{V}{3f_1}}. \tag*{(43)} x ( N ) = 6 ϕ ′ , y ( N ) = H 1 3 f 1 V . (43)
Neglecting baryons, the effective dark energy and dark matter density parameters can then be expressed as
Ω ~ D E , e f f = ρ D E , e f f 3 H 2 f 1 = ( 1 − 6 Q 2 ) x 2 + y 2 + 2 6 Q x , (44) \widetilde{\Omega}_{{\rm DE,eff}}=\frac{\rho_{{\rm DE,eff}}}{3H^2f_1}=(1-6Q^2)x^2+y^2+2\sqrt{6}Qx, \tag*{(44)} Ω DE , eff = 3 H 2 f 1 ρ DE , eff = ( 1 − 6 Q 2 ) x 2 + y 2 + 2 6 Q x , (44)
Ω ~ d m = 1 − Ω ~ D E , e f f = 1 − ( 1 − 6 Q 2 ) x 2 − y 2 − 2 6 Q x , (45) \widetilde{\Omega}_{\rm dm}=1-\widetilde{\Omega}_{{\rm DE,eff}}=1-(1-6Q^2)x^2-y^2-2\sqrt{6}Qx, \tag*{(45)} Ω dm = 1 − Ω DE , eff = 1 − ( 1 − 6 Q 2 ) x 2 − y 2 − 2 6 Q x , (45)
and E P E_P E P is given by
E P = 1 + 6 Q d m x − 3 Q d m [ 1 − ( 1 − 6 Q 2 ) x 2 − y 2 − 2 6 Q x ] Y ( k ) Q + Q d m f g . (46) E_P=1+\sqrt{6}Q_{\rm dm}x-3Q_{\rm dm}\left[1-(1-6Q^2)x^2-y^2-2\sqrt{6}Qx\right]Y(k)\frac{Q+Q_{\rm dm}}{f_g}. \tag*{(46)} E P = 1 + 6 Q dm x − 3 Q dm [ 1 − ( 1 − 6 Q 2 ) x 2 − y 2 − 2 6 Q x ] Y ( k ) f g Q + Q dm . (46)
We now want to investigate how strong the couplings Q d m Q_{\rm dm} Q dm and Q Q Q can be constrained by present and future surveys. To this end, two scenarios are considered: the CDE case in which only Q d m Q_{\rm dm} Q dm is nonzero, and the DCDE case in which both Q Q Q and Q d m Q_{\rm dm} Q dm are nonzero.
The first scenario is interesting since it corresponds to the simplest non-trivial coupled model, in which f 1 = 1 f_1=1 f 1 = 1 , consequently Q = 0 Q=0 Q = 0 and Ω d m = Ω ~ d m \Omega_{\rm dm}=\widetilde{\Omega}_{\rm dm} Ω dm = Ω dm . In addition, we further assume a 2 m ϕ 2 ≪ k 2 a^2m_\phi^2\ll k^2 a 2 m ϕ 2 ≪ k 2 , so that Y → 1 Y\to1 Y → 1 . This is expected, as the mass of the dark energy scalar field needs to be extremely small in order to drive the late-time accelerated expansion of the Universe. The main motivation to examine closely this particular case, besides its simplicity, is that it has been shown to provide a good fit to the recent DESI data (cf. Refs [24 , 25 , 34 , 35 , 36 , 37 ] ). Under these assumptions, E P E_P E P can be further simplified and becomes k k k -independent:
E P = 1 + Q d m ϕ ′ − 3 Ω d m Q d m 2 f g = 1 + 6 Q d m x − 3 ( 1 − x 2 − y 2 ) Q d m 2 f g . (47) \begin{aligned}
E_P&=1+Q_{\rm dm}\phi'-3\frac{\Omega_{\rm dm}Q_{\rm dm}^2}{f_g}\\
&=1+\sqrt{6}Q_{\rm dm}x-3\frac{(1-x^2-y^2)Q_{\rm dm}^2}{f_g}. \tag*{(47)}
\end{aligned} E P = 1 + Q dm ϕ ′ − 3 f g Ω dm Q dm 2 = 1 + 6 Q dm x − 3 f g ( 1 − x 2 − y 2 ) Q dm 2 . (47)
Since we expect the behavior of the model to be close to that of Λ \Lambda Λ CDM to match observations, we can adopt the Λ \Lambda Λ CDM-like approximation f g ≈ Ω d m γ f_g\approx\Omega_{\rm dm}^{\gamma} f g ≈ Ω dm γ for the growth rate, with Ω d m \Omega_{\rm dm} Ω dm evaluated using the background evolution of the CDE model and γ = 0.545 \gamma=0.545 γ = 0.545 .
By further introducing the dimensionless coupling{^5} β d m = − 6 2 Q d m \beta_{\rm dm}=-\frac{\sqrt{6}}{2}Q_{\rm dm} β dm = − 2 6 Q dm , we finally obtain E P E_P E P for the coupled dark energy scenario:
E P ( N ) = 1 − 2 β d m x − 2 β d m 2 ( 1 − x 2 − y 2 ) 1 − γ . (48) E_P(N)=1-2\beta_{\rm dm}x-2\beta_{\rm dm}^2(1-x^2-y^2)^{1-\gamma}. \tag*{(48)} E P ( N ) = 1 − 2 β dm x − 2 β dm 2 ( 1 − x 2 − y 2 ) 1 − γ . (48)
The evolution of x x x and y y y depends on the shape of the potential. For generality, we employ the analytical approximate solutions for x x x and y y y valid around the present time derived in Ref. [59 ] for a linear potential V ≈ V 0 − V 0 μ ( ϕ − ϕ 0 ) V\approx V_0-V_0\mu(\phi-\phi_0) V ≈ V 0 − V 0 μ ( ϕ − ϕ 0 ) (where μ \mu μ is constant, and we set ϕ 0 = 0 \phi_0=0 ϕ 0 = 0 as before), which corresponds to the first-order Taylor expansion of a generic potential and satisfies V ′ ( N ) = − V 0 μ ϕ ′ = − V 0 μ 6 x V'(N)=-V_0\mu\phi'=-V_0\mu\sqrt{6}x V ′ ( N ) = − V 0 μ ϕ ′ = − V 0 μ 6 x . This yields
x ( N ) = 1 6 e 3 N E ~ ( N ) { 4 6 β d m [ e 3 N i E ~ ( N i ) + Ω m 0 Ω Λ 0 ln ( e 3 N i / 2 ( E ~ ( N i ) − Ω Λ 0 ) e 3 N / 2 ( E ~ ( N ) − Ω Λ 0 ) ) ] + μ [ Ω m 0 Ω Λ 0 ln ( e 3 N / 2 ( E ~ ( N ) − Ω Λ 0 ) e 3 N i / 2 ( E ~ ( N i ) − Ω Λ 0 ) ) + e 3 N E ~ ( N ) − e 3 N i E ~ ( N i ) ] } . (49) x(N)=\frac{1}{\sqrt{6}e^{3N}\widetilde{E}(N)}\left\{\frac{4}{\sqrt{6}}\beta_{\rm dm}\left[e^{3N_i}\widetilde{E}(N_i)+\frac{\Omega_{m0}}{\sqrt{\Omega_{\Lambda0}}}\ln\left(\frac{e^{3N_i/2}(\widetilde{E}(N_i)-\sqrt{\Omega_{\Lambda0}})}{e^{3N/2}(\widetilde{E}(N)-\sqrt{\Omega_{\Lambda0}})}\right)\right]+\mu\left[\frac{\Omega_{m0}}{\sqrt{\Omega_{\Lambda0}}}\ln\left(\frac{e^{3N/2}(\widetilde{E}(N)-\sqrt{\Omega_{\Lambda0}})}{e^{3N_i/2}(\widetilde{E}(N_i)-\sqrt{\Omega_{\Lambda0}})}\right)+e^{3N}\widetilde{E}(N)-e^{3N_i}\widetilde{E}(N_i)\right]\right\}. \tag*{(49)} x ( N ) = 6 e 3 N E ( N ) 1 { 6 4 β dm [ e 3 N i E ( N i ) + Ω Λ0 Ω m 0 ln ( e 3 N /2 ( E ( N ) − Ω Λ0 ) e 3 N i /2 ( E ( N i ) − Ω Λ0 ) ) ] + μ [ Ω Λ0 Ω m 0 ln ( e 3 N i /2 ( E ( N i ) − Ω Λ0 ) e 3 N /2 ( E ( N ) − Ω Λ0 ) ) + e 3 N E ( N ) − e 3 N i E ( N i ) ] } . (49)
V ( N ) = V 0 { 1 − 4 6 β d m N μ + μ ( 8 6 β d m − 2 μ 3 Ω Λ 0 ) [ Ω Λ 0 Ω m 0 e 3 N i E ~ ( N i ) ( E ~ ( N ) − 1 ) − ln ( e 3 N i / 2 ( E ~ ( N i ) − Ω Λ 0 ) 1 − Ω Λ 0 ) + E ~ ( N ) ln ( e 3 N i / 2 ( E ~ ( N i ) − Ω Λ 0 ) e 3 N / 2 ( E ~ ( N ) − Ω Λ 0 ) ) ] } . (50) V(N)=V_0\left\{1-\frac{4}{\sqrt{6}}\beta_{\rm dm}N\mu+\mu\left(\frac{\frac{8}{\sqrt{6}}\beta_{\rm dm}-2\mu}{3\sqrt{\Omega_{\Lambda0}}}\right)\left[\frac{\sqrt{\Omega_{\Lambda0}}}{\Omega_{m0}}e^{3N_i}\widetilde{E}(N_i)\left(\widetilde{E}(N)-1\right)-\ln\left(\frac{e^{3N_i/2}(\widetilde{E}(N_i)-\sqrt{\Omega_{\Lambda0}})}{1-\sqrt{\Omega_{\Lambda0}}}\right)+\widetilde{E}(N)\ln\left(\frac{e^{3N_i/2}(\widetilde{E}(N_i)-\sqrt{\Omega_{\Lambda0}})}{e^{3N/2}(\widetilde{E}(N)-\sqrt{\Omega_{\Lambda0}})}\right)\right]\right\}. \tag*{(50)} V ( N ) = V 0 { 1 − 6 4 β dm N μ + μ ( 3 Ω Λ0 6 8 β dm − 2 μ ) [ Ω m 0 Ω Λ0 e 3 N i E ( N i ) ( E ( N ) − 1 ) − ln ( 1 − Ω Λ0 e 3 N i /2 ( E ( N i ) − Ω Λ0 ) ) + E ( N ) ln ( e 3 N /2 ( E ( N ) − Ω Λ0 ) e 3 N i /2 ( E ( N i ) − Ω Λ0 ) ) ] } . (50)
where N i N_i N i is the initial time chosen at N i = − 3 N_i=-3 N i = − 3 , Ω Λ 0 = 1 − Ω m 0 \Omega_{\Lambda0}=1-\Omega_{m0} Ω Λ0 = 1 − Ω m 0 , V 0 = 3 H 0 2 Ω Λ 0 V_0=3H_0^2\Omega_{\Lambda0} V 0 = 3 H 0 2 Ω Λ0 , and E ~ ( N ) = Ω m 0 e − 3 N + 1 − Ω m 0 \widetilde{E}(N)=\sqrt{\Omega_{m0}e^{-3N}+1-\Omega_{m0}} E ( N ) = Ω m 0 e − 3 N + 1 − Ω m 0 . The initial condition of x x x is solved from the scaling solution during the ϕ \phi ϕ MDE phase, with x i ( N = − 3 ) = 2 3 β d m x_i(N=-3)=\frac{2}{3}\beta_{\rm dm} x i ( N = − 3 ) = 3 2 β dm .
Next, for the doubly coupled dark energy, where both Q Q Q and Q d m Q_{\rm dm} Q dm are active, we further define the dimensionless coupling β = − 6 2 Q \beta=-\frac{\sqrt{6}}{2}Q β = − 2 6 Q associated with the universal coupling. Using the same approximation of Y → 1 Y\to1 Y → 1 and f g ≈ Ω ~ E P γ f_g\approx\widetilde{\Omega}_{\rm EP}^{\gamma} f g ≈ Ω EP γ , E P E_{\rm P} E P is given by
E P = 1 − 2 β d m x − 2 β d m ( β + β d m ) [ 1 − ( 1 − 4 β 2 ) x 2 − y 2 + 4 β x ] 1 − γ . (51) E_{\rm P}=1-2\beta_{\rm dm}x-2\beta_{\rm dm}(\beta+\beta_{\rm dm})\left[1-(1-4\beta^2)x^2-y^2+4\beta x\right]^{1-\gamma}. \tag*{(51)} E P = 1 − 2 β dm x − 2 β dm ( β + β dm ) [ 1 − ( 1 − 4 β 2 ) x 2 − y 2 + 4 β x ] 1 − γ . (51)
The quantities x x x and V V V in the doubly coupled model can be analytically approximated by:
x ( N ) = e 3 N i E ~ ( N i ) e 3 N E ~ ( N ) 2 ( β + β d m ) 3 − 4 β ( β + β d m ) + 1 3 ( 1 − 4 β 2 ) [ 4 β + 3 2 μ ] [ 1 − e 3 N i E ~ ( N i ) e 3 N E ~ ( N ) ] + Ω m 0 3 ( 1 − 4 β 2 ) 1 − Ω m 0 e 3 N E ~ ( N ) [ 2 β d m − 2 β − 3 2 μ ] × ln [ e 3 N i / 2 ( E ~ ( N i ) − 1 − Ω m 0 ) e 3 N / 2 ( E ~ ( N ) − 1 − Ω m 0 ) ] . (52) \begin{aligned}
x(N)={}&\frac{e^{3N_i}\widetilde{E}(N_i)}{e^{3N}\widetilde{E}(N)}
\frac{2(\beta+\beta_{\rm dm})}{3-4\beta(\beta+\beta_{\rm dm})}
+\frac{1}{3(1-4\beta^2)}
\left[4\beta+\sqrt{\frac{3}{2}}\mu\right]
\left[1-\frac{e^{3N_i}\widetilde{E}(N_i)}{e^{3N}\widetilde{E}(N)}\right]\\
&+\frac{\Omega_{m0}}{3(1-4\beta^2)\sqrt{1-\Omega_{m0}e^{3N}\widetilde{E}(N)}}
\left[2\beta_{\rm dm}-2\beta-\sqrt{\frac{3}{2}}\mu\right]\\
&\qquad\times\ln\left[
\frac{e^{3N_i/2}\left(\widetilde{E}(N_i)-\sqrt{1-\Omega_{m0}}\right)}
{e^{3N/2}\left(\widetilde{E}(N)-\sqrt{1-\Omega_{m0}}\right)}
\right]. \tag*{(52)}
\end{aligned} x ( N ) = e 3 N E ( N ) e 3 N i E ( N i ) 3 − 4 β ( β + β dm ) 2 ( β + β dm ) + 3 ( 1 − 4 β 2 ) 1 [ 4 β + 2 3 μ ] [ 1 − e 3 N E ( N ) e 3 N i E ( N i ) ] + 3 ( 1 − 4 β 2 ) 1 − Ω m 0 e 3 N E ( N ) Ω m 0 [ 2 β dm − 2 β − 2 3 μ ] × ln e 3 N /2 ( E ( N ) − 1 − Ω m 0 ) e 3 N i /2 ( E ( N i ) − 1 − Ω m 0 ) . (52)
V ( N ) = V 0 { 1 − 2 6 μ ( β + β d m ) 3 ( 1 − 4 β 2 ) N + 2 6 μ e 3 N i E ~ ( N i ) 3 Ω m 0 [ x ( N i ) − 4 β + 3 2 μ 3 ( 1 − 4 β 2 ) ] [ E ~ ( N ) − 1 ] + 2 6 μ [ 2 β d m − 2 β − 3 2 μ ] 9 ( 1 − 4 β 2 ) Ω Λ 0 [ E ~ ( N ) ln ( e 3 N i / 2 [ E ~ ( N i ) − Ω Λ 0 ] e 3 N / 2 [ E ~ ( N ) − Ω Λ 0 ] ) − ln ( e 3 N i / 2 [ E ~ ( N i ) − Ω Λ 0 ] 1 − Ω Λ 0 ) ] } . (53) \begin{aligned}
V(N)=V_0\Bigg\{1
&-\frac{2\sqrt{6}\mu(\beta+\beta_{\rm dm})}{3(1-4\beta^2)}N
+\frac{2\sqrt{6}\mu e^{3N_i}\widetilde{E}(N_i)}{3\Omega_{m0}}
\left[x(N_i)-\frac{4\beta+\sqrt{\frac{3}{2}}\mu}{3(1-4\beta^2)}\right]
\left[\widetilde{E}(N)-1\right]\\
&+\frac{2\sqrt{6}\mu\left[2\beta_{\rm dm}-2\beta-\sqrt{\frac{3}{2}}\mu\right]}
{9(1-4\beta^2)\sqrt{\Omega_{\Lambda0}}}
\Bigg[
\widetilde{E}(N)\ln\left(
\frac{e^{3N_i/2}\left[\widetilde{E}(N_i)-\sqrt{\Omega_{\Lambda0}}\right]}
{e^{3N/2}\left[\widetilde{E}(N)-\sqrt{\Omega_{\Lambda0}}\right]}
\right)\\
&\qquad\qquad\qquad\qquad
-\ln\left(
\frac{e^{3N_i/2}\left[\widetilde{E}(N_i)-\sqrt{\Omega_{\Lambda0}}\right]}
{1-\sqrt{\Omega_{\Lambda0}}}
\right)
\Bigg]\Bigg\}. \tag*{(53)}
\end{aligned} V ( N ) = V 0 { 1 − 3 ( 1 − 4 β 2 ) 2 6 μ ( β + β dm ) N + 3 Ω m 0 2 6 μ e 3 N i E ( N i ) x ( N i ) − 3 ( 1 − 4 β 2 ) 4 β + 2 3 μ [ E ( N ) − 1 ] + 9 ( 1 − 4 β 2 ) Ω Λ0 2 6 μ [ 2 β dm − 2 β − 2 3 μ ] [ E ( N ) ln e 3 N /2 [ E ( N ) − Ω Λ0 ] e 3 N i /2 [ E ( N i ) − Ω Λ0 ] − ln 1 − Ω Λ0 e 3 N i /2 [ E ( N i ) − Ω Λ0 ] ] } . (53)
The derivation of Eqs. (52) and (53) is provided in App. § 2.2 . Moreover, the phase-space system for the doubly coupled model including baryons is summarized in App. § 2.3 . The expressions Eqs. (48) –(53) will be implemented in our forecasts. Next, we describe our forecasting methodology in Sec. IV and provide the forecast results for the constant couplings β \beta β and β d m \beta_{\rm dm} β dm in Sec. VI.
Fisher forecast for the couplings
In Ref. [41 ] , it has been shown that the auto- and cross-matter linear power spectra for two galaxy populations (with the subscripts B, F, for ”bright” and ”faint”, respectively) up to order λ = H / k \lambda=\mathcal{H}/k λ = H / k can be written as
P F F = ( 1 + β F μ 2 ) 2 β B β F S g 2 B P . (54) P_{\rm FF}=(1+\beta_{\rm F}\mu^2)^2\frac{\beta_{\rm B}}{\beta_{\rm F}}S_g^2BP. \tag*{(54)} P FF = ( 1 + β F μ 2 ) 2 β F β B S g 2 B P . (54)
P B B = ( 1 + β B μ 2 ) 2 β F β B S g 2 B P . (55) P_{\rm BB}=(1+\beta_{\rm B}\mu^2)^2\frac{\beta_{\rm F}}{\beta_{\rm B}}S_g^2BP. \tag*{(55)} P BB = ( 1 + β B μ 2 ) 2 β B β F S g 2 B P . (55)
P F B = [ ( 1 + β F μ 2 ) ( 1 + β B μ 2 ) + i λ μ ( τ 1 + μ 2 τ 2 ) ] S g 2 B P . (56) P_{\rm FB}=\left[(1+\beta_{\rm F}\mu^2)(1+\beta_{\rm B}\mu^2)+i\lambda\mu(\tau_1+\mu^2\tau_2)\right]S_g^2BP. \tag*{(56)} P FB = [ ( 1 + β F μ 2 ) ( 1 + β B μ 2 ) + iλ μ ( τ 1 + μ 2 τ 2 ) ] S g 2 B P . (56)
in which μ \mu μ is the angle between the galaxy pair and the line of sight to the observer, β B , F = f g / b B , F \beta_{\rm B,F}=f_g/b_{\rm B,F} β B , F = f g / b B , F (b B , F b_{\rm B,F} b B , F being the time-dependent linear galaxy bias for each population), P P P is the power spectrum at z = 0 z=0 z = 0 , and B ≡ b B ( z ) b F ( z ) G 2 ( z ) B\equiv b_{\rm B}(z)b_{\rm F}(z)G^2(z) B ≡ b B ( z ) b F ( z ) G 2 ( z ) , with G ≡ δ m ( z ) / δ m 0 G\equiv\delta_m(z)/\delta_{m0} G ≡ δ m ( z ) / δ m 0 denoting the linear growth factor. The spectra are further multiplied by factors accounting for the Finger-of-God effect
S g 2 ( k , μ , z ) = e − 1 2 k 2 μ 2 σ g 2 . (57) S_g^2(k,\mu,z)=e^{-\frac{1}{2}k^2\mu^2\sigma_g^2}. \tag*{(57)} S g 2 ( k , μ , z ) = e − 2 1 k 2 μ 2 σ g 2 . (57)
in which the damping strengths σ g \sigma_g σ g are left as free parameters in each redshift bin. Moreover, the parameters τ 1 \tau_1 τ 1 and τ 2 \tau_2 τ 2 appearing in the cross-spectrum are given by
τ 1 = β B α B − β F α F − E P ( β B − β F ) , (58) \tau_1=\beta_{\rm B}\alpha_{\rm B}-\beta_{\rm F}\alpha_{\rm F}-E_{\rm P}(\beta_{\rm B}-\beta_{\rm F}), \tag*{(58)} τ 1 = β B α B − β F α F − E P ( β B − β F ) , (58)
τ 2 = β B β F ( α B − α F ) . (59) \tau_2=\beta_{\rm B}\beta_{\rm F}(\alpha_{\rm B}-\alpha_{\rm F}). \tag*{(59)} τ 2 = β B β F ( α B − α F ) . (59)
where we have defined for each galaxy population
α ≡ 1 − 5 s + 5 s − 2 H r − H ′ H + f e v o l . (60) \alpha\equiv1 - 5s + \frac{5s - 2}{\mathcal{H}r} - \frac{\mathcal{H}'}{\mathcal{H}} + f^{\mathrm{evol}}. \tag*{(60)} α ≡ 1 − 5 s + H r 5 s − 2 − H H ′ + f evol . (60)
Here, s s s is the magnification bias, accounting for the fact that surveys are flux limited, and f e v o l f^{\mathrm{evol}} f evol is the evolution bias, encoding the evolution of galaxies. Here we simply set f e v o l = 0 f^{\mathrm{evol}} = 0 f evol = 0 , as in our previous analysis [41 , 42 ] . Moreover, we assume in the following that α B , F \alpha_{B,F} α B , F can be accurately constrained, as both the magnification and the evolution bias can be measured from the average number of galaxies in a given sample [60 , 61 ] . By grouping the power spectra Eqs. (54) –(56) into polynomials of μ a λ b \mu^a\lambda^b μ a λ b , one can identify measurable model-independent quantities. From τ 1 \tau_1 τ 1 in Eq. (58) , we see then explicitly that E P E_P E P can only be measured if the two galaxy populations have different biases. This is the fundamental reason why we need two distinct source populations.
The spectra in Eqs. (54) –(56) encode also the Alcock-Paczynski (AP) effect through k k k and μ \mu μ , relating to the reference cosmology through
k = α A P k r ; μ = μ r h α A P ; α A P = 1 d μ r 2 ( h 2 d 2 − 1 ) + 1 , (61) k = \alpha_{\mathrm{AP}} k_r; \qquad\mu= \frac{\mu_r h}{\alpha_{\mathrm{AP}}}; \qquad\alpha_{\mathrm{AP}} = \frac{1}{d}\sqrt{\mu_r^2(h^2d^2-1)+1}, \tag*{(61)} k = α AP k r ; μ = α AP μ r h ; α AP = d 1 μ r 2 ( h 2 d 2 − 1 ) + 1 , (61)
with k r k_r k r and μ r \mu_r μ r evaluated in the reference cosmology. These expressions depend on h ≡ E / E r h \equiv E/E_r h ≡ E / E r and on d ≡ L A / L A r d \equiv L_A/L_{Ar} d ≡ L A / L A r , in which E ( z ) ≡ H ( z ) / H 0 E(z) \equiv H(z)/H_0 E ( z ) ≡ H ( z ) / H 0 is the dimensionless Hubble function and L A ( z ) ≡ H 0 D A ( z ) L_A(z) \equiv H_0D_A(z) L A ( z ) ≡ H 0 D A ( z ) is the dimensionless comoving angular diameter distance. Since d d d is degenerate with the power spectrum at the linear level, we will consider the combination h d ≡ h d h_d \equiv hd h d ≡ h d as a parameter in our analysis. In Ref. [41 ] , a forecast analysis has been performed for a set of parameters { h d , d , B , σ g , β B , β F , E P } \{h_d,d,B,\sigma_g,\beta_B,\beta_F,E_P\} { h d , d , B , σ g , β B , β F , E P } at each redshift. In addition, the power spectrum shape was parametrized in several wavebands in the first redshift bin z 1 z_1 z 1 and evolved with the free function B ( z ) B(z) B ( z ) . The k k k -range was chosen to be k ∈ ( 0.01 − 0.12 ) h / M p c k \in(0.01-0.12)\,h/\mathrm{Mpc} k ∈ ( 0.01 − 0.12 ) h / Mpc with intervals Δ k = 0.01 h / M p c \Delta k = 0.01\,h/\mathrm{Mpc} Δ k = 0.01 h / Mpc , for a total of twelve P ( k ) P(k) P ( k ) waveband parameters. Following Ref. [41 ] , we assume uninformative flat priors for all parameters except d d d , with a 1% Gaussian prior on the latter, based on current constraints (see, e.g., the Pantheon+ supernova catalog [62 ] and App. § 2.5 ), while fixing B ( z 1 ) = 1 B(z_1) = 1 B ( z 1 ) = 1 to break the degeneracy with P ( k , z 1 ) P(k,z_1) P ( k , z 1 ) (in other words, we only measure ratios of B B B at different redshifts).
Here, we follow the methodology of Refs. [41 , 42 ] , with only slightly different survey specifications (see § 1.4 ), to obtain constraints on { h d , d , B , σ g , β B , β F , E P } \{h_d,d,B,\sigma_g,\beta_B,\beta_F,E_P\} { h d , d , B , σ g , β B , β F , E P } within each redshift bin. We then project these redshift-dependent constraints onto our constant model parameters. For the CDE scenario, we project onto { β d m , Ω m 0 } \{\beta_{\mathrm{dm}},\Omega_{m0}\} { β dm , Ω m 0 } , while for the DCDE scenario, the parameters are instead { β , β d m , Ω m 0 } \{\beta,\beta_{\mathrm{dm}},\Omega_{m0}\} { β , β dm , Ω m 0 } . The couplings β \beta β and β d m \beta_{\mathrm{dm}} β dm enter both E P ( z ) E_P(z) E P ( z ) (see Eqs. (48) and (51) ) and the expansion rate E ( z ) E(z) E ( z ) (see Eqs. (120) and (78) ), with the latter modifying h d h_d h d and d d d through the AP effect. We therefore project the marginalised constraints on { E P , h d , d } \{E_P,h_d,d\} { E P , h d , d } in each redshift bin onto the CDE and DCDE model parameters, respectively6 . We summarize the details of the Fisher projection in App. § 2.4 .
Survey specifications and fiducials
Here we consider three upcoming galaxy surveys and forecast how well they can constrain the coupling strengths: the DESI Bright Galaxy Sample (BGS) [63 ] , a Euclid-like stage-IV Hα \alpha α spectroscopic survey [64 ] , and the SKA Phase 2 HI spectroscopic galaxy survey [65 ] . In addition, we compute the fiducial matter power spectrum at the lowest redshift bin, background quantities, and linear growth rate assuming a standard flat Λ \Lambda Λ CDM cosmology based on the Planck 2018 results [66 ] . To forecast the couplings, we also need to specify the initial time N i N_i N i and the fiducial values of β \beta β , β d m \beta_{\mathrm{dm}} β dm , and Ω m 0 \Omega_{m0} Ω m 0 in the expressions of x x x and y y y . Following Ref. [59 ] , we adopt N i = − 3 N_i=-3 N i = − 3 , deep in the matter-dominated era, where x i x_i x i can be determined from the ϕ \phi ϕ MDE scaling solution (see App. § 2.2 ). Finally, we take the fiducial values Ω m 0 = 0.315 \Omega_{m0}=0.315 Ω m 0 = 0.315 , β = 0.1 \beta=0.1 β = 0.1 , β d m = 0.2 \beta_{\mathrm{dm}}=0.2 β dm = 0.2 , V 0 = 3 ( 1 − Ω m 0 ) ≈ 2.055 V_0=3(1-\Omega_{m0})\approx2.055 V 0 = 3 ( 1 − Ω m 0 ) ≈ 2.055 and μ = 0.145 \mu=0.145 μ = 0.145 . The fiducial values for the couplings are larger than current constraints obtained by comparison with CMB and other cosmological data (see e.g. [24 , 26 , 32 ] ). These analyses, however, assume a constant coupling throughout the entire cosmic evolution, while here we only consider late-time data.
For each survey, we split the full galaxy sample into two populations (bright and faint) within each redshift bin. This is done by choosing a flux threshold that divides the sample equally between the two populations. As in [41 ] , we choose the fiducial biases to be b g , B = b g + 0.5 b_{g,B}=b_g+0.5 b g , B = b g + 0.5 and b g , F = b g − 0.5 b_{g,F}=b_g-0.5 b g , F = b g − 0.5 , where b g b_g b g is the mean galaxy bias of the full sample. This choice is consistent with the measurements reported in [61 ] .
Table II. Survey specifications for a Euclid-like spectroscopic survey, covering the redshift range 0.9 ≤ z ≤ 1.8 0.9 \le z \le1.8 0.9 ≤ z ≤ 1.8 in four redshift bins with boundaries at z = 0.9 , 1.1 , 1.3 , 1.5 , 1.8 z = 0.9, 1.1, 1.3, 1.5, 1.8 z = 0.9 , 1.1 , 1.3 , 1.5 , 1.8 . The quantities below are derived from Refs. [60 , 67 ] .
The remaining survey specifications for DESI BGS, Euclid, and SKA2 are identical to those in Ref. [42 ] . We list them explicitly in Table I –Table III . These include the redshift bins, survey volumes, galaxy number densities, and average galaxy biases, as well as the treatment of galaxy redshift uncertainties (associated with σ g \sigma_g σ g ) and the prescriptions for magnification and evolution biases.
Table I. Survey specifications for DESI BGS covering the redshift range 0 ≤ z ≤ 0.5 0 \le z \le0.5 0 ≤ z ≤ 0.5 in bins of width Δ z = 0.1 \Delta z = 0.1 Δ z = 0.1 . Here we list the central redshift of each bin. In this table and hereafter, we express the galaxy number density n g n_g n g in units of 10 − 3 ( h / M p c ) 3 10^{-3}(h/\mathrm{Mpc})^3 1 0 − 3 ( h / Mpc ) 3 , assuming h = 0.67 h = 0.67 h = 0.67 , and the survey volumes V V V in ( G p c / h ) 3 (\mathrm{Gpc}/h)^3 ( Gpc / h ) 3 . The fiducial bias are assumed as b g , B , b g , F = b g ± 0.5 b_{g,B},b_{g,F} = b_g \pm0.5 b g , B , b g , F = b g ± 0.5 . The galaxy number density n g n_g n g and the magnification bias s B , s F s_B,s_F s B , s F , are obtained from [60 ] . The observed absolute magnitude limit is chosen to be m c = 19.5 m_c = 19.5 m c = 19.5 . The evolution bias, f e v o l f_{\mathrm{evol}} f evol , is set to zero.
Table III. SKA2 survey specifications covering the redshift range 0.2 ≤ z ≤ 0.8 0.2 \le z \le0.8 0.2 ≤ z ≤ 0.8 in bins of width Δ z = 0.1 \Delta z = 0.1 Δ z = 0.1 . The quantities below are derived from Refs. [60 , 65 ] .
The forecasting results presented in the next section are based on the aforementioned survey specifications.
Results and discussion
We present our forecasting results in Tables IV–VII for DESI, Euclid, SKA2, and, just to illustrate the effect of the number density, a survey mimicking SKA2 but with ten times as many galaxies. All the results shown in these tables are absolute errors obtained after marginalising over all other parameters. No priors are taken into account except those on d d d , as illustrated in Sec. IV. In each table, we present both model-independent constraints on { h d , d , E P } \{h_d,d,E_P\} { h d , d , E P } within each redshift bin and model-independent constraints on the parameters of the CDE and DCDE models obtained by projection. We emphasize that here we only employ the forms of E P E_P E P and the expansion rate E E E (which enters h d h_d h d and d d d ) corresponding to the models onto which the constraints are projected, while the underlying analysis remains
Table V. Same as Table. Table IV, but for the Euclid spectroscopic survey.
general as far as the cosmological model is concerned. By doing so, we expect our results to be more robust than those obtained from analyses carried out assuming specific models, at the price of possibly weaker constraints. Moreover, these constraints are fully obtained from the late-time Universe, in contrast to tighter constraints obtained including CMB measurements and assuming constant couplings throughout the evolution.
Among the three surveys, SKA2 could provide the tightest constraints on h d h_d h d , d d d and E P E_{\mathrm{P}} E P , consistent with the results reported in Refs. [41 , 42 ] , thanks to its large number of detected galaxies. Projecting these constraints onto the CDE model, we obtain relative 1 σ 1\sigma 1 σ errors of 22% and 3.5% on β d m \beta_{\mathrm{dm}} β dm and Ω m 0 \Omega_{\mathrm{m0}} Ω m0 , respectively. Let us first discuss the constraints on β d m \beta_{\mathrm{dm}} β dm , which directly reflects EP violation. From Eq. (48) , the uncertainty of β d m \beta_{\mathrm{dm}} β dm inferred from E P E_{\mathrm{P}} E P at a single redshift bin, while keeping Ω m 0 \Omega_{\mathrm{m0}} Ω m0 fixed, can be approximated by
σ β d m ≈ ∣ 1 ∂ E P / ∂ β d m ∣ f i d σ E P . (62) \sigma_{\beta_{\mathrm{dm}}} \approx\left|\frac{1}{\partial E_{\mathrm{P}}/\partial\beta_{\mathrm{dm}}}\right|_{\mathrm{fid}} \sigma_{E_{\mathrm{P}}}. \tag*{(62)} σ β dm ≈ ∂ E P / ∂ β dm 1 fid σ E P . (62)
For SKA2, ∣ ∂ E P / ∂ β d m ∣ f i d ∼ 0.9 − 1.12 |\partial E_{\mathrm{P}}/\partial\beta_{\mathrm{dm}}|_{\mathrm{fid}} \sim0.9 - 1.12 ∣ ∂ E P / ∂ β dm ∣ fid ∼ 0.9 − 1.12 across the six redshift bins, evaluated at the fiducial. Assuming independent bins, the 1 σ 1\sigma 1 σ uncertainty of β d m \beta_{\mathrm{dm}} β dm from E P E_{\mathrm{P}} E P alone is approximately σ β d m ≈ 0.044 \sigma_{\beta_{\mathrm{dm}}} \approx0.044 σ β dm ≈ 0.044 , corresponding to a relative 1 σ 1\sigma 1 σ uncertainty of 22.1%. This estimate is comparable to our forecast constraint of σ β d m = 0.044 \sigma_{\beta_{\mathrm{dm}}} = 0.044 σ β dm = 0.044 , as listed in Table VI , despite neglecting the dependence of h d h_d h d and d d d on β d m \beta_{\mathrm{dm}} β dm and correlation between bins. We also verified that, when including this, accounting for correlations between bins, and fixing Ω m 0 \Omega_{\mathrm{m0}} Ω m0 , one can obtain σ β d m = 0.034 \sigma_{\beta_{\mathrm{dm}}} = 0.034 σ β dm = 0.034 . These results indicate that, for the upcoming Stage IV surveys considered here, E P E_{\mathrm{P}} E P measured from SKA2 can provide useful constraints on β d m \beta_{\mathrm{dm}} β dm within the CDE model. The constraints on β d m \beta_{\mathrm{dm}} β dm can be further improved at the cost of reducing
Table VI. Same as Table. Table IV, but for the SKA2.
model independence. For instance, we explored the case in which the power spectrum is fixed at fiducial, while leaving the remaining parameters free to vary. In this case, the relative 1 σ 1\sigma 1 σ error on β d m \beta_{\rm dm} β dm for SKA2 is reduced to 19%.
On the other hand, β d m \beta_{\rm dm} β dm not only affects perturbations by modifying the Euler equation and the growth rate, but also modifies the background expansion rate E E E . This introduces a degeneracy between the two parameters and increases the marginalised uncertainty on Ω m 0 \Omega_{\rm m0} Ω m0 when projecting the constraints on h d h_d h d and d d d onto the model parameters. Indeed, the correlation coefficients between β d m \beta_{\rm dm} β dm and Ω m 0 \Omega_{\rm m0} Ω m0 are − 0.93 -0.93 − 0.93 , − 0.99 -0.99 − 0.99 , and − 0.65 -0.65 − 0.65 for DESI, Euclid, and SKA2, respectively. We also repeat the analysis but with a Λ \Lambda Λ CDM expansion, and the correlation coefficients become − 0.0094 -0.0094 − 0.0094 , 0.0019 0.0019 0.0019 , 0.056 0.056 0.056 for the three surveys. For SKA2, which is the most promising survey in this paper, the marginalised 1 σ 1\sigma 1 σ error on β d m \beta_{\rm dm} β dm remains approximately the same at 0.044, but the absolute error on Ω m 0 \Omega_{\rm m0} Ω m0 decreases to 0.0032.
We also present in Figure 2 the two-dimensional confidence contours between β d m \beta_{\rm dm} β dm and Ω m 0 \Omega_{\rm m0} Ω m0 for all three surveys, together with their one-dimensional marginalised posterior distributions. In all three cases, the two parameters are negatively correlated. Moreover, the eigenvector ( c 1 , c 2 ) (c_1,c_2) ( c 1 , c 2 ) associated with the smallest eigenvalue of the covariance matrix gives the linear combination C = c 1 β d m + c 2 Ω m 0 C=c_1\beta_{\rm dm}+c_2\Omega_{\rm m0} C = c 1 β dm + c 2 Ω m0 with the minimum variance. This eigenvector points along the short axis of the confidence ellipse. We find ( c 1 , c 2 ) ≃ ( 0.27 , 0.96 ) , ( 0.21 , 0.98 ) , ( 0.17 , 0.99 ) (c_1,c_2)\simeq(0.27,0.96),(0.21,0.98),(0.17,0.99) ( c 1 , c 2 ) ≃ ( 0.27 , 0.96 ) , ( 0.21 , 0.98 ) , ( 0.17 , 0.99 ) for DESI, Euclid, and SKA2, respectively. The corresponding relative 1 σ 1\sigma 1 σ uncertainties on C C C are 29%, 4.6%, 2.5%. Thus, despite the large uncertainties on the individual parameters, this particular combination can be much better constrained, especially for Euclid, as can be seen directly from Figure 2 .
Figure 2. Forecast constraints on β d m \beta_{\rm dm} β dm and Ω m 0 \Omega_{\rm m0} Ω m0 for the CDE model. Here we show both one-dimensional posterior distributions of these two parameters, together with their 2D contour plots, within the 1 σ 1\sigma 1 σ and 2 σ 2\sigma 2 σ confidence region. The fiducial values are indicated as a star.
For the DCDE model, one can further constrain the universal coupling β \beta β , as it enters E P E_P E P differently from β d m \beta_{\rm dm} β dm . The constraints are weaker, as expected, since we now have an additional parameter to vary. In the most tightly constrained case of SKA2, we obtain relative errors of 430%, 130%, 24% for β \beta β , β d m \beta_{\rm dm} β dm and Ω m 0 \Omega_{\rm m0} Ω m0 , respectively. Note that even though β \beta β also alters the expansion and the growth rate, only β d m \beta_{\rm dm} β dm is associated with the violation of the EP. To improve the constraints on β \beta β , one can add more cosmological probes, such as the CMB measurements from the early Universe.
We also explored the effect of increasing the galaxy number density. For instance, increasing the SKA2 galaxy number density by a factor of ten could substantially improve the projected constraints, as listed in Table VII. For CDE, the 1 σ 1\sigma 1 σ uncertainties decrease to σ β d m = 0.0081 \sigma_{\beta_{\rm dm}}=0.0081 σ β dm = 0.0081 and σ Ω m 0 = 0.003 \sigma_{\Omega_{\rm m0}}=0.003 σ Ω m0 = 0.003 , corresponding to relative errors of 4.1% and 0.95%, respectively. For DCDE, the 1 σ 1\sigma 1 σ uncertainties decrease to σ β = 0.21 \sigma_\beta= 0.21 σ β = 0.21 , σ β d m = 0.12 \sigma_{\beta_{\rm dm}} = 0.12 σ β dm = 0.12 , and σ Ω m 0 = 0.038 \sigma_{\Omega_{\rm m0}} = 0.038 σ Ω m0 = 0.038 , improving by roughly a factor of two. Moreover, constraints on E P E_{\rm P} E P can be significantly improved, reaching 3% at z = 0.45 z = 0.45 z = 0.45 .
Lastly, this methodology can be extended to other EP-violating models or to alternative shapes of the scalar-field potential. Once the corresponding E P E_{\rm P} E P and expansion rate E E E are identified, one can directly employ the framework presented here to investigate constraints on those model parameters.
Conclusion
In this paper, we studied model-independent cosmological tests of the weak equivalence principle in the context of a doubly coupled dark energy model, in which a scalar field couples simultaneously to gravity (associated with β \beta β ) and dark matter (associated with β d m \beta_{\rm dm} β dm ). When β \beta β vanishes, this model reduces to coupled quintessence. We derived the background and linear perturbation equations in this DCDE scenario and showed that the effective potential is modified by a Yukawa correction. Additionally, both couplings modify the background expansion, the growth equation, the Euler equation, and therefore enter E P E_{\rm P} E P . As expected, E P E_{\rm P} E P reduces to unity when β d m \beta_{\rm dm} β dm is switched off, even with a non-vanishing universal coupling. Employing a dynamical system approach, we then obtained analytical expressions for E P E_{\rm P} E P in terms of the two couplings, assuming a linear potential.
We then investigated whether these couplings, assumed here to be constant during the relevant late-time epoch, could be detected by next-generation cosmological surveys, namely DESI, Euclid and SKA2. To preserve robustness, we projected the model-independent constraints on { E P , h d , d } \{E_{\rm P}, h_{\rm d}, d\} { E P , h d , d } , obtained similarly to our previous analyses [41 , 42 ] using two different galaxy populations, onto the CDE and DCDE model parameters. This could give us, for the first time, constraints on β d m \beta_{\rm dm} β dm derived from direct cosmological tests of EP violations.
To summarise our results, our forecasts reveal the potential of EP measurements from SKA2 to constrain β d m \beta_{\rm dm} β dm within the CDE model. However, obtaining tight constraints remains challenging for DESI and Euclid, as well as for the parameters of the DCDE model. This is due both to the limited sensitivity of E P E_{\rm P} E P to the couplings and to degeneracies between the couplings and other parameters. Increasing the galaxy number density could significantly help improve the constraints.
Our results suggest several directions for further investigation. First, one could explore how the constraints depend on more complex shapes of the scalar-field potential. This requires to solve for x x x and y y y numerically. Second, constraints on E P E_{\rm P} E P could be improved by selecting galaxy populations with a larger bias difference, as shown in Ref. [41 ] . One could also investigate whether extending the analysis to more than two galaxy populations provides further improvements. Finally, adding more cosmological probes could help reduce the degeneracy between β \beta β and β d m \beta_{\rm dm} β dm and improve constraints on the latter in the DCDE model, and therefore worth investigating this model through a full MCMC analysis. We leave these investigations for future work.
Acknowledgments
The authors thank C. Bonvin and S. Castello for useful comments. LA acknowledges support by DFG under Germany’s Excellence Strategy EXC 2181/1 - 390900948 (the Heidelberg STRUCTURES Excellence Cluster) and under Project 554679582 ”GeoGrav: Cosmological Geometry and Gravity with nonlinear physics” and by Tamkeen under the NYUAD Research Institute grant CG014.
Yukawa corrections to the effective gravitational potentials
From Eqs. (28) , (29) and (31) , the effective potential for dark matter particles Ψ ^ d m ≡ Ψ + Q d m φ \hat{\Psi}_{\rm dm} \equiv\Psi+ Q_{\rm dm}\varphi Ψ ^ dm ≡ Ψ + Q dm φ is given by
Ψ ^ d m = Ψ + Q d m φ = − 3 λ 2 2 [ ( 2 Q Y ( k ) ( Q + Q d m ) + 1 ) δ ~ d m δ d m + ( 2 Q 2 Y ( k ) + 1 ) δ ~ b δ b ] − 3 λ 2 2 [ 2 Q d m Y ( k ) ( Q + Q d m ) δ ~ d m δ d m + 2 Y ( k ) Q Q d m δ ~ b δ b ] = − 3 λ 2 2 { [ 1 + 2 Y ( Q + Q d m ) 2 ] δ ~ d m δ d m + [ 1 + 2 Y Q ( Q + Q d m ) ] δ ~ b δ b } . (63) \begin{aligned}
\hat{\Psi}_{\rm dm} &= \Psi+ Q_{\rm dm}\varphi\\
&= -\frac{3\lambda^2}{2}\left[(2QY(k)(Q+Q_{\rm dm})+1)\tilde{\delta}_{\rm dm}\delta_{\rm dm}+(2Q^2Y(k)+1)\tilde{\delta}_{b}\delta_b\right] \\
&\quad-\frac{3\lambda^2}{2}\left[2Q_{\rm dm}Y(k)(Q+Q_{\rm dm})\tilde{\delta}_{\rm dm}\delta_{\rm dm}+2Y(k)QQ_{\rm dm}\tilde{\delta}_{b}\delta_b\right] \\
&= -\frac{3\lambda^2}{2}\left\{\left[1+2Y(Q+Q_{\rm dm})^2\right]\tilde{\delta}_{\rm dm}\delta_{\rm dm}+\left[1+2YQ(Q+Q_{\rm dm})\right]\tilde{\delta}_{b}\delta_b\right\}. \tag*{(63)}
\end{aligned} Ψ ^ dm = Ψ + Q dm φ = − 2 3 λ 2 [ ( 2 Q Y ( k ) ( Q + Q dm ) + 1 ) δ ~ dm δ dm + ( 2 Q 2 Y ( k ) + 1 ) δ ~ b δ b ] − 2 3 λ 2 [ 2 Q dm Y ( k ) ( Q + Q dm ) δ ~ dm δ dm + 2 Y ( k ) Q Q dm δ ~ b δ b ] = − 2 3 λ 2 { [ 1 + 2 Y ( Q + Q dm ) 2 ] δ ~ dm δ dm + [ 1 + 2 Y Q ( Q + Q dm ) ] δ ~ b δ b } . (63)
Following Ref. [28 ] , the density contrast for a dark matter particle in the DCDE model, located in empty space at the origin, is given by
Ω ~ d m δ d m = ρ ˉ d m ρ c r i t ρ d m − ρ ˉ d m ρ ˉ d m = δ ρ d m ρ c r i t = δ ρ d m 3 f 1 H 2 . (64) \widetilde{\Omega}_{\rm dm}\delta_{\rm dm}=\frac{\bar{\rho}_{\rm dm}}{\rho_{\rm crit}}\frac{\rho_{\rm dm}-\bar{\rho}_{\rm dm}}{\bar{\rho}_{\rm dm}}=\frac{\delta\rho_{\rm dm}}{\rho_{\rm crit}}=\frac{\delta\rho_{\rm dm}}{3f_1H^2}. \tag*{(64)} Ω dm δ dm = ρ crit ρ ˉ dm ρ ˉ dm ρ dm − ρ ˉ dm = ρ crit δ ρ dm = 3 f 1 H 2 δ ρ dm . (64)
On the other hand we know that
m d m ( a ) = m d m , 0 e Q d m [ ϕ ( a ) − ϕ 0 ] = a 3 ∫ d 3 x δ ρ d m ( x ) . (65) m_{\rm dm}(a)=m_{{\rm dm},0}e^{Q_{\rm dm}[\phi(a)-\phi_0]}=a^3\int d^3x\delta\rho_{\rm dm}(\boldsymbol{x}). \tag*{(65)} m dm ( a ) = m dm , 0 e Q dm [ ϕ ( a ) − ϕ 0 ] = a 3 ∫ d 3 x δ ρ dm ( x ) . (65)
Solving this gives,
δ ρ d m ( x ) = m d m a 3 δ D ( 3 ) ( x ) , (66) \delta\rho_{\rm dm}(\boldsymbol{x})=\frac{m_{\rm dm}}{a^3}\delta_D^{(3)}(\boldsymbol{x}), \tag*{(66)} δ ρ dm ( x ) = a 3 m dm δ D ( 3 ) ( x ) , (66)
with δ D ( 3 ) ( x ) \delta_D^{(3)}(\boldsymbol{x}) δ D ( 3 ) ( x ) the Dirac delta function satisfying ∫ d 3 x δ D ( 3 ) ( x ) = 1 \int d^3x\delta_D^{(3)}(\boldsymbol{x})=1 ∫ d 3 x δ D ( 3 ) ( x ) = 1 . Transforming to Fourier space we obtain
δ ρ d m ( k ) = m d m a 3 . (67) \delta\rho_{\rm dm}(\boldsymbol{k})=\frac{m_{\rm dm}}{a^3}. \tag*{(67)} δ ρ dm ( k ) = a 3 m dm . (67)
Plugging back into Eq. (64)
3 2 λ 2 Ω ~ d m δ d m = 3 2 λ 2 δ ρ d m 3 f 1 H 2 = m d m 2 f 1 k 2 a , (68) \frac{3}{2}\lambda^2\widetilde{\Omega}_{\rm dm}\delta_{\rm dm}=\frac{3}{2}\lambda^2\frac{\delta\rho_{\rm dm}}{3f_1H^2}=\frac{m_{\rm dm}}{2f_1k^2a}, \tag*{(68)} 2 3 λ 2 Ω dm δ dm = 2 3 λ 2 3 f 1 H 2 δ ρ dm = 2 f 1 k 2 a m dm , (68)
and Ψ ^ d m ← d m \widehat{\Psi}_{{\rm dm}\leftarrow{\rm dm}} Ψ dm ← dm (i.e. the effective potential sourced from this single dark matter particle) becomes
Ψ ^ d m ← d m ( k ) = − m d m 2 a f 1 [ 1 k 2 + 2 ( Q + Q d m ) 2 k 2 + a 2 m ϕ 2 ] , (69) \widehat{\Psi}_{{\rm dm}\leftarrow{\rm dm}}(\boldsymbol{k})=-\frac{m_{\rm dm}}{2af_1}\left[\frac{1}{k^2}+\frac{2(Q+Q_{\rm dm})^2}{k^2+a^2m_\phi^2}\right], \tag*{(69)} Ψ dm ← dm ( k ) = − 2 a f 1 m dm [ k 2 1 + k 2 + a 2 m ϕ 2 2 ( Q + Q dm ) 2 ] , (69)
and transforming to real space this gives (with r = a x \boldsymbol{r}=a\boldsymbol{x} r = a x denoting the physical distance)
Ψ ^ d m ← d m ( r ) = − G m d m f 1 r [ 1 + 2 ( Q + Q d m ) 2 e − m ϕ r ] . (70) \widehat{\Psi}_{{\rm dm}\leftarrow{\rm dm}}(r)=-\frac{Gm_{\rm dm}}{f_1r}\left[1+2(Q+Q_{\rm dm})^2e^{-m_\phi r}\right]. \tag*{(70)} Ψ dm ← dm ( r ) = − f 1 r G m dm [ 1 + 2 ( Q + Q dm ) 2 e − m ϕ r ] . (70)
Likewise, for a baryon particle located at the origin, its effective potential acting on a dark matter test particle is given by
Ψ ^ d m ← b ( r ) = − G m b f 1 r [ 1 + 2 Q ( Q + Q d m ) e − m ϕ r ] . (71) \widehat{\Psi}_{{\rm dm}\leftarrow{\rm b}}(r)=-\frac{Gm_{\rm b}}{f_1r}\left[1+2Q(Q+Q_{\rm dm})e^{-m_\phi r}\right]. \tag*{(71)} Ψ dm ← b ( r ) = − f 1 r G m b [ 1 + 2 Q ( Q + Q dm ) e − m ϕ r ] . (71)
On the other hand, for a baryon test particle the effective potential is simply Ψ ^ b = Ψ \widehat{\Psi}_{\rm b}=\Psi Ψ b = Ψ , one can find this from Eq. (30) . Consider a single baryon particle as the source, then the effective potential is
Ψ ^ b ← b ( k ) = − m b 2 a f 1 [ 1 k 2 + 2 Q 2 k 2 + a 2 m ϕ 2 ] , (72) \widehat{\Psi}_{{\rm b}\leftarrow{\rm b}}(\boldsymbol{k})=-\frac{m_{\rm b}}{2af_1}\left[\frac{1}{k^2}+\frac{2Q^2}{k^2+a^2m_\phi^2}\right], \tag*{(72)} Ψ b ← b ( k ) = − 2 a f 1 m b [ k 2 1 + k 2 + a 2 m ϕ 2 2 Q 2 ] , (72)
and transforming to real space this gives
Ψ ^ b ← b ( r ) = − G m b f 1 r [ 1 + 2 Q 2 e − m ϕ r ] . (73) \widehat{\Psi}_{{\rm b}\leftarrow{\rm b}}(r)=-\frac{Gm_{\rm b}}{f_1r}\left[1+2Q^2e^{-m_\phi r}\right]. \tag*{(73)} Ψ b ← b ( r ) = − f 1 r G m b [ 1 + 2 Q 2 e − m ϕ r ] . (73)
Approximate analytical solutions for x ( N ) = ϕ ′ 6 x(N)=\frac{\phi'}{\sqrt{6}} x ( N ) = 6 ϕ ′ and y ( N ) = 1 H V 3 f 1 y(N)=\frac{1}{H}\sqrt{\frac{V}{3f_1}} y ( N ) = H 1 3 f 1 V in the doubly coupled model
It has been shown in [59 ] that one can find approximate analytical solutions for x x x and y y y with a linear potential and with only the coupling β d m \beta_{\rm dm} β dm switched on. Moreover, the comparison between the analytical and numerical solutions in [59 ] (see their Fig. 3) indicates that assuming a background close to Λ \Lambda Λ CDM is a good approximation for small couplings. In this appendix, we further derive the analytical solutions for x x x and y y y in DCDE with the same linear potential, while depending on both β \beta β and β d m \beta_{\rm dm} β dm . In this appendix we also neglect the baryon contribution. For a dynamical system approach including the baryons, see the next App. C.
We first derive the expansion rate E E E in this model. Following a similar derivation as in [59 ] , taking the time derivative of Eq. (9) and plug in Eqs. (16) and (17) , we have
H ′ H = − 3 2 ( 1 + p D E , e f f 3 f 1 H 2 ) (74) \frac{H'}{H}=-\frac{3}{2}\left(1+\frac{p_{\mathrm{DE,eff}}}{3f_1H^2}\right) \tag*{(74)} H H ′ = − 2 3 ( 1 + 3 f 1 H 2 p DE , eff ) (74)
where p D E , e f f p_{\mathrm{DE,eff}} p DE , eff is effective pressure of the dark energy component. Since we set κ = 1 \kappa=1 κ = 1 , the energy-momentum tensor of the effective dark energy reads
T μ ν D E , e f f = T μ ν ϕ + ∇ μ ∇ ν f 1 − g μ ν □ f 1 (75) T_{\mu\nu}^{\mathrm{DE,eff}}=T_{\mu\nu}^{\phi}+\nabla_\mu\nabla_\nu f_1-g_{\mu\nu}\Box f_1 \tag*{(75)} T μν DE , eff = T μν ϕ + ∇ μ ∇ ν f 1 − g μν □ f 1 (75)
and therefore
p D E , e f f = 1 2 ( 1 − 6 Q 2 ) f 1 H 2 ϕ ′ 2 − V + H 2 [ f 1 ′ ′ + ( 2 + H ′ H ) f 1 ′ ] (76) p_{\mathrm{DE,eff}}=\frac{1}{2}(1-6Q^2)f_1H^2\phi'^2-V+H^2\left[f_1''+\left(2+\frac{H'}{H}\right)f_1'\right] \tag*{(76)} p DE , eff = 2 1 ( 1 − 6 Q 2 ) f 1 H 2 ϕ ′2 − V + H 2 [ f 1 ′′ + ( 2 + H H ′ ) f 1 ′ ] (76)
Expressing the pressure in terms of x x x and y y y , and using the background equations to cancel f 1 ′ ′ f_1'' f 1 ′′ and H ′ / H H'/H H ′ / H , we obtain
p D E , e f f 3 f 1 H 2 = 4 3 β 2 + ( 1 − 4 β 2 ) [ ( 1 − 4 3 β 2 ) x 2 − y 2 − 4 3 β x ] + 4 3 β β d m [ 1 − ( 1 − 4 β 2 ) x 2 − y 2 + 4 β x ] − 2 6 3 β V , ϕ V y 2 (77) \frac{p_{\mathrm{DE,eff}}}{3f_1H^2}=\frac{4}{3}\beta^2+(1-4\beta^2)\left[\left(1-\frac{4}{3}\beta^2\right)x^2-y^2-\frac{4}{3}\beta x\right]+\frac{4}{3}\beta\beta_{\rm dm}\left[1-(1-4\beta^2)x^2-y^2+4\beta x\right]-\frac{2\sqrt{6}}{3}\beta\frac{V_{,\phi}}{V}y^2 \tag*{(77)} 3 f 1 H 2 p DE , eff = 3 4 β 2 + ( 1 − 4 β 2 ) [ ( 1 − 3 4 β 2 ) x 2 − y 2 − 3 4 β x ] + 3 4 β β dm [ 1 − ( 1 − 4 β 2 ) x 2 − y 2 + 4 β x ] − 3 2 6 β V V , ϕ y 2 (77)
It is easy to find that p D E , e f f 3 f 1 H 2 = x 2 − y 2 \frac{p_{\mathrm{DE,eff}}}{3f_1H^2}=x^2-y^2 3 f 1 H 2 p DE , eff = x 2 − y 2 if β = 0 \beta=0 β = 0 . The expansion rate E E E can be solved from Eq. (74) :
E ( N ) = exp [ − 3 2 ∫ 0 N ( 1 + p D E , e f f ( N ′ ) 3 f 1 ( N ′ ) H 2 ( N ′ ) ) d N ′ ] (78) E(N)=\exp\left[-\frac{3}{2}\int_0^N\left(1+\frac{p_{\mathrm{DE,eff}}(N')}{3f_1(N')H^2(N')}\right)\mathrm{d}N'\right] \tag*{(78)} E ( N ) = exp [ − 2 3 ∫ 0 N ( 1 + 3 f 1 ( N ′ ) H 2 ( N ′ ) p DE , eff ( N ′ ) ) d N ′ ] (78)
by plugging in Eq. (77) . We emphasize that this integral is an exact solution of E E E . Next, we proceed to solve for the approximate solutions for x ( N ) x(N) x ( N ) and y ( N ) y(N) y ( N ) . Substituting ϕ ′ = 6 x \phi'=\sqrt{6}x ϕ ′ = 6 x ; − 2 6 β = Q = − f 1 , ϕ f 1 -\frac{2}{\sqrt{6}}\beta=Q=-\frac{f_{1,\phi}}{f_1} − 6 2 β = Q = − f 1 f 1 , ϕ ; − 2 6 β d m = Q d m = f 2 , ϕ f 2 -\frac{2}{\sqrt{6}}\beta_{\rm dm}=Q_{\rm dm}=\frac{f_{2,\phi}}{f_2} − 6 2 β dm = Q dm = f 2 f 2 , ϕ ; V = V 0 − V 0 μ ϕ V=V_0-V_0\mu\phi V = V 0 − V 0 μ ϕ , the Klein-Gordon equation (11) becomes
x ′ + ( 3 + H ′ H ) x + 2 β x 2 = 1 1 − 4 β 2 [ V 0 μ 6 f 1 H 2 + 2 β ( 2 + H ′ H ) + β d m Ω ~ d m ] (79) x'+\left(3+\frac{H'}{H}\right)x+2\beta x^2=\frac{1}{1-4\beta^2}\left[\frac{V_0\mu}{\sqrt{6}f_1H^2}+2\beta\left(2+\frac{H'}{H}\right)+\beta_{\rm dm}\widetilde{\Omega}_{\rm dm}\right] \tag*{(79)} x ′ + ( 3 + H H ′ ) x + 2 β x 2 = 1 − 4 β 2 1 [ 6 f 1 H 2 V 0 μ + 2 β ( 2 + H H ′ ) + β dm Ω dm ] (79)
which is still exact. Assuming the background is close to Λ \Lambda Λ CDM (which indicates small β , β d m \beta,\beta_{\rm dm} β , β dm ) and w D E , e f f = p D E , e f f / ρ D E , e f f ≈ − 1 w_{\mathrm{DE,eff}}=p_{\mathrm{DE,eff}}/\rho_{\mathrm{DE,eff}}\approx-1 w DE , eff = p DE , eff / ρ DE , eff ≈ − 1 , we have f 1 , f 2 ≈ 1 f_1,f_2\approx1 f 1 , f 2 ≈ 1 ; V 0 ≈ 3 H 0 2 ( 1 − Ω m 0 ) V_0\approx3H_0^2(1-\Omega_{\rm m0}) V 0 ≈ 3 H 0 2 ( 1 − Ω m0 ) and
Ω ~ d m ≈ Ω m 0 Ω m 0 + ( 1 − Ω m 0 ) e 3 N ; H ′ H ≈ − 3 2 Ω ~ d m ≈ − 3 Ω m 0 2 [ Ω m 0 + ( 1 − Ω m 0 ) e 3 N ] ; H 2 ≈ H 0 2 ( Ω m 0 e − 3 N + 1 − Ω m 0 ) (80) \widetilde{\Omega}_{\rm dm}\approx\frac{\Omega_{\rm m0}}{\Omega_{\rm m0}+(1-\Omega_{\rm m0})e^{3N}};\qquad\frac{H'}{H}\approx-\frac{3}{2}\widetilde{\Omega}_{\rm dm}\approx-\frac{3\Omega_{\rm m0}}{2[\Omega_{\rm m0}+(1-\Omega_{\rm m0})e^{3N}]};\qquad H^2\approx H_0^2(\Omega_{\rm m0}e^{-3N}+1-\Omega_{\rm m0}) \tag*{(80)} Ω dm ≈ Ω m0 + ( 1 − Ω m0 ) e 3 N Ω m0 ; H H ′ ≈ − 2 3 Ω dm ≈ − 2 [ Ω m0 + ( 1 − Ω m0 ) e 3 N ] 3 Ω m0 ; H 2 ≈ H 0 2 ( Ω m0 e − 3 N + 1 − Ω m0 ) (80)
In addition, within the low-redshift range considered in this paper, we can assume that the kinetic energy of the scalar field is much smaller than its potential energy, allowing us to neglect the term quadratic in x x x in Eq. (79) , which now can be reduced to
x ′ + ( 3 − 3 Ω m 0 e − 3 N 2 E ~ 2 ) x = 1 1 − 4 β 2 [ V 0 μ 6 H 0 2 E ~ 2 + 4 β + ( β d m − 3 β ) Ω m 0 e − 3 N E ~ 2 ] (81) x'+\left(3-\frac{3\Omega_{\rm m0}e^{-3N}}{2\widetilde{E}^2}\right)x=\frac{1}{1-4\beta^2}\left[\frac{V_0\mu}{\sqrt{6}H_0^2\widetilde{E}^2}+4\beta+(\beta_{\rm dm}-3\beta)\frac{\Omega_{\rm m0}e^{-3N}}{\widetilde{E}^2}\right] \tag*{(81)} x ′ + ( 3 − 2 E 2 3 Ω m0 e − 3 N ) x = 1 − 4 β 2 1 [ 6 H 0 2 E 2 V 0 μ + 4 β + ( β dm − 3 β ) E 2 Ω m0 e − 3 N ] (81)
recall that E ~ \widetilde{E} E is previously defined as E ~ ( N ) = Ω m 0 e − 3 N + 1 − Ω m 0 \widetilde{E}(N)=\sqrt{\Omega_{\rm m0}e^{-3N}+1-\Omega_{\rm m0}} E ( N ) = Ω m0 e − 3 N + 1 − Ω m0 . Eq. (81) has analytical solutions
x ( N ) = e 3 N i E ~ ( N i ) e 3 N E ~ ( N ) x i ( N i ) + 1 3 ( 1 − 4 β 2 ) [ 4 β + 3 2 μ ] [ 1 − e 3 N i E ~ ( N i ) e 3 N E ~ ( N ) ] + Ω m 0 3 ( 1 − 4 β 2 ) 1 − Ω m 0 e 3 N E ~ ( N ) [ 2 β d m − 2 β − 3 2 μ ] × ln [ e 3 N i / 2 ( E ~ ( N i ) − 1 − Ω m 0 ) e 3 N / 2 ( E ~ ( N ) − 1 − Ω m 0 ) ] . (82) \begin{aligned}
x(N)={}&\frac{e^{3N_i}\widetilde{E}(N_i)}{e^{3N}\widetilde{E}(N)}x_i(N_i)+\frac{1}{3(1-4\beta^2)}\left[4\beta+\sqrt{\frac{3}{2}}\mu\right]\left[1-\frac{e^{3N_i}\widetilde{E}(N_i)}{e^{3N}\widetilde{E}(N)}\right]\\
&+\frac{\Omega_{\rm m0}}{3(1-4\beta^2)\sqrt{1-\Omega_{\rm m0}e^{3N}\widetilde{E}(N)}}\left[2\beta_{\rm dm}-2\beta-\sqrt{\frac{3}{2}}\mu\right]\\
&\times\ln\left[\frac{e^{3N_i/2}\left(\widetilde{E}(N_i)-\sqrt{1-\Omega_{\rm m0}}\right)}{e^{3N/2}\left(\widetilde{E}(N)-\sqrt{1-\Omega_{\rm m0}}\right)}\right]. \tag*{(82)}
\end{aligned} x ( N ) = e 3 N E ( N ) e 3 N i E ( N i ) x i ( N i ) + 3 ( 1 − 4 β 2 ) 1 [ 4 β + 2 3 μ ] [ 1 − e 3 N E ( N ) e 3 N i E ( N i ) ] + 3 ( 1 − 4 β 2 ) 1 − Ω m0 e 3 N E ( N ) Ω m0 [ 2 β dm − 2 β − 2 3 μ ] × ln e 3 N /2 ( E ( N ) − 1 − Ω m0 ) e 3 N i /2 ( E ( N i ) − 1 − Ω m0 ) . (82)
Using V ′ ( N ) = − V 0 μ ϕ ′ = − V 0 μ 6 x V'(N)=-V_0\mu\phi'=-V_0\mu\sqrt{6}x V ′ ( N ) = − V 0 μ ϕ ′ = − V 0 μ 6 x , one can further obtain
V ( N ) = V 0 { 1 − 2 6 μ ( β + β d m ) 3 ( 1 − 4 β 2 ) N + 2 6 μ e 3 N i E ~ ( N i ) 3 Ω m 0 [ x ( N i ) − 4 β + 3 2 μ 3 ( 1 − 4 β 2 ) ] [ E ~ ( N ) − 1 ] + 2 6 μ [ 2 β d m − 2 β − 3 2 μ ] 9 ( 1 − 4 β 2 ) Ω Λ 0 [ E ~ ( N ) ln ( e 3 N i / 2 [ E ~ ( N i ) − Ω Λ 0 ] e 3 N / 2 [ E ~ ( N ) − Ω Λ 0 ] ) − ln ( e 3 N i / 2 [ E ~ ( N i ) − Ω Λ 0 ] 1 − Ω Λ 0 ) ] } . (83) \begin{aligned}
V(N)=V_0\Biggl\{1&-\frac{2\sqrt{6}\mu(\beta+\beta_{\rm dm})}{3(1-4\beta^2)}N+\frac{2\sqrt{6}\mu e^{3N_i}\widetilde{E}(N_i)}{3\Omega_{m0}}\left[x(N_i)-\frac{4\beta+\sqrt{\frac{3}{2}}\mu}{3(1-4\beta^2)}\right]\left[\widetilde{E}(N)-1\right]\\
&+\frac{2\sqrt{6}\mu\left[2\beta_{\rm dm}-2\beta-\sqrt{\frac{3}{2}}\mu\right]}{9(1-4\beta^2)\sqrt{\Omega_{\Lambda0}}}\Biggl[\widetilde{E}(N)\ln\left(\frac{e^{3N_i/2}\left[\widetilde{E}(N_i)-\sqrt{\Omega_{\Lambda0}}\right]}{e^{3N/2}\left[\widetilde{E}(N)-\sqrt{\Omega_{\Lambda0}}\right]}\right)\\
&\qquad-\ln\left(\frac{e^{3N_i/2}\left[\widetilde{E}(N_i)-\sqrt{\Omega_{\Lambda0}}\right]}{1-\sqrt{\Omega_{\Lambda0}}}\right)\Biggr]\Biggr\}. \tag*{(83)}
\end{aligned} V ( N ) = V 0 { 1 − 3 ( 1 − 4 β 2 ) 2 6 μ ( β + β dm ) N + 3 Ω m 0 2 6 μ e 3 N i E ( N i ) x ( N i ) − 3 ( 1 − 4 β 2 ) 4 β + 2 3 μ [ E ( N ) − 1 ] + 9 ( 1 − 4 β 2 ) Ω Λ0 2 6 μ [ 2 β dm − 2 β − 2 3 μ ] [ E ( N ) ln e 3 N /2 [ E ( N ) − Ω Λ0 ] e 3 N i /2 [ E ( N i ) − Ω Λ0 ] − ln 1 − Ω Λ0 e 3 N i /2 [ E ( N i ) − Ω Λ0 ] ] } . (83)
One can verify that, in the β = 0 \beta=0 β = 0 limit, Eqs. (82) and (B10) reduce to Eqs. (49) and (50) , as expected. The final step is to find the initial condition x i ≡ x ( N i ) x_i\equiv x(N_i) x i ≡ x ( N i ) as required by Eq. (82) . We can obtain from the fixed point that satisfies the conditions x ′ = 0 x'=0 x ′ = 0 and y = 0 y=0 y = 0 . This is similar to deriving the ϕ \phi ϕ MDE phase in the CDE scenario, and Eq. (45) reads
Ω ~ d m ( N i ) = ρ d m ( N i ) 3 f 1 ( N i ) H 2 ( N i ) = 1 − ( 1 − 6 Q 2 ) x i 2 − 2 6 Q x i = 1 − ( 1 − 4 β 2 ) x i 2 + 4 β x i , (84) \widetilde{\Omega}_{\rm dm}(N_i)=\frac{\rho_{\rm dm}(N_i)}{3f_1(N_i)H^2(N_i)}=1-\left(1-6Q^2\right)x_i^2-2\sqrt{6}Qx_i=1-\left(1-4\beta^2\right)x_i^2+4\beta x_i, \tag*{(84)} Ω dm ( N i ) = 3 f 1 ( N i ) H 2 ( N i ) ρ dm ( N i ) = 1 − ( 1 − 6 Q 2 ) x i 2 − 2 6 Q x i = 1 − ( 1 − 4 β 2 ) x i 2 + 4 β x i , (84)
Note that Eq. (84) follows from the exact constraint in the DCDE case. We no longer assume a background close to Λ \Lambda Λ CDM here or in the following derivation. For a constant non-vanishing density fraction of dark matter
Ω ~ d m ′ Ω ~ d m = ρ d m ′ ρ d m − f 1 ′ f 1 − 2 H ′ H = 0 , (85) \frac{\widetilde{\Omega}'_{\rm dm}}{\widetilde{\Omega}_{\rm dm}}=\frac{\rho'_{\rm dm}}{\rho_{\rm dm}}-\frac{f'_1}{f_1}-2\frac{H'}{H}=0, \tag*{(85)} Ω dm Ω dm ′ = ρ dm ρ dm ′ − f 1 f 1 ′ − 2 H H ′ = 0 , (85)
it is easy to find that
ρ d m ′ ρ d m = Q d m ϕ ′ − 3 = − 2 β d m x i − 3 ; f 1 ′ f 1 = f 1 , ϕ ϕ ′ f 1 = − 2 Q ϕ ′ = 4 β x i , (86) \frac{\rho'_{\rm dm}}{\rho_{\rm dm}}=Q_{\rm dm}\phi'-3=-2\beta_{\rm dm}x_i-3;\qquad\frac{f'_1}{f_1}=\frac{f_{1,\phi}\phi'}{f_1}=-2Q\phi'=4\beta x_i, \tag*{(86)} ρ dm ρ dm ′ = Q dm ϕ ′ − 3 = − 2 β dm x i − 3 ; f 1 f 1 ′ = f 1 f 1 , ϕ ϕ ′ = − 2 Q ϕ ′ = 4 β x i , (86)
and therefore
H ′ H = − 3 2 − ( β d m + 2 β ) x i . (87) \frac{H'}{H}=-\frac{3}{2}-(\beta_{\rm dm}+2\beta)x_i. \tag*{(87)} H H ′ = − 2 3 − ( β dm + 2 β ) x i . (87)
Substituting Eqs. (84) and (B14) into Eq. (79) , we find
x i = 2 ( β + β d m ) 3 − 4 β ( β + β d m ) . (88) x_i=\frac{2(\beta+\beta_{\rm dm})}{3-4\beta(\beta+\beta_{\rm dm})}. \tag*{(88)} x i = 3 − 4 β ( β + β dm ) 2 ( β + β dm ) . (88)
Setting β = 0 \beta=0 β = 0 recovers the CDE solution x i = 2 β d m / 3 x_i=2\beta_{\rm dm}/3 x i = 2 β dm /3 , as expected, while for β d m = 0 \beta_{\rm dm}=0 β dm = 0 the result agrees with Eq. 32 in Ref. [68]. In addition, the constant effective dark energy fraction during the scaling phase of the DCDE model is given by
Ω ~ D E , e f f = ( 1 − 4 β 2 ) x i 2 − 4 β x i = ( 1 − 4 β 2 ) [ 2 ( β + β d m ) 3 − 4 β ( β + β d m ) ] 2 − 8 β ( β + β d m ) 3 − 4 β ( β + β d m ) . (89) \widetilde{\Omega}_{\rm DE,eff}=(1-4\beta^2)x_i^2-4\beta x_i=(1-4\beta^2)\left[\frac{2(\beta+\beta_{\rm dm})}{3-4\beta(\beta+\beta_{\rm dm})}\right]^2-\frac{8\beta(\beta+\beta_{\rm dm})}{3-4\beta(\beta+\beta_{\rm dm})}. \tag*{(89)} Ω DE , eff = ( 1 − 4 β 2 ) x i 2 − 4 β x i = ( 1 − 4 β 2 ) [ 3 − 4 β ( β + β dm ) 2 ( β + β dm ) ] 2 − 3 − 4 β ( β + β dm ) 8 β ( β + β dm ) . (89)
We also discuss here a particular case. From Eq. (87) , we see that, in this regime, we can identify the effective EoS
w e f f ≡ − 1 − 2 3 H ′ H = 2 3 ( β d m + 2 β ) x i . (90) w_{\rm eff}\equiv-1-\frac{2}{3}\frac{H'}{H}=\frac{2}{3}(\beta_{\rm dm}+2\beta)x_i. \tag*{(90)} w eff ≡ − 1 − 3 2 H H ′ = 3 2 ( β dm + 2 β ) x i . (90)
Interestingly, for β d m = − 2 β \beta_{\rm dm}=-2\beta β dm = − 2 β , one can have an expansion rate identical to a pure matter-dominated universe, i.e. w e f f = 0 w_{\rm eff}=0 w eff = 0 , even in presence of dark energy, i.e. x i ≠ 0 x_i\ne0 x i = 0 . From the equivalent form of Eq. (31) (again neglecting baryons)
δ d m ′ ′ + ( 1 + H ′ H − 2 β d m x ) δ d m ′ − 3 2 [ 1 + 4 3 Y ( k ) ( β + β d m ) 2 ] [ 1 − ( 1 − 4 β 2 ) x 2 − y 2 + 4 β x ] δ d m = 0 , (91) \delta''_{\rm dm}+\left(1+\frac{H'}{H}-2\beta_{\rm dm}x\right)\delta'_{\rm dm}-\frac{3}{2}\left[1+\frac{4}{3}Y(k)(\beta+\beta_{\rm dm})^2\right]\left[1-(1-4\beta^2)x^2-y^2+4\beta x\right]\delta_{\rm dm}=0, \tag*{(91)} δ dm ′′ + ( 1 + H H ′ − 2 β dm x ) δ dm ′ − 2 3 [ 1 + 3 4 Y ( k ) ( β + β dm ) 2 ] [ 1 − ( 1 − 4 β 2 ) x 2 − y 2 + 4 β x ] δ dm = 0 , (91)
and plugging in this particular case with y = 0 y=0 y = 0 , Y ( k ) = 1 Y(k)=1 Y ( k ) = 1 , β d m = − 2 β \beta_{\rm dm}=-2\beta β dm = − 2 β and x = x i x=x_i x = x i it becomes
δ d m ′ ′ + ( 1 2 − 8 β 2 3 + 4 β 2 ) δ d m ′ − 9 − 4 β 2 2 ( 3 + 4 β 2 ) δ d m = 0. (92) \delta''_{\rm dm}+\left(\frac{1}{2}-\frac{8\beta^2}{3+4\beta^2}\right)\delta'_{\rm dm}-\frac{9-4\beta^2}{2(3+4\beta^2)}\delta_{\rm dm}=0. \tag*{(92)} δ dm ′′ + ( 2 1 − 3 + 4 β 2 8 β 2 ) δ dm ′ − 2 ( 3 + 4 β 2 ) 9 − 4 β 2 δ dm = 0. (92)
For β 2 < 9 / 4 \beta^2 < 9/4 β 2 < 9/4 , the growing mode is given by δ d m ∝ a \delta_{\rm dm} \propto a δ dm ∝ a and thus f g = 1 f_g = 1 f g = 1 , again the same as during the matter-dominated era in a Λ \Lambda Λ CDM Universe. The E P E_P E P , however, is given by
E P = 9 ( 1 − 4 β 2 ) ( 3 + 4 β 2 ) 2 (93) E_P = \frac{9(1-4\beta^2)}{(3+4\beta^2)^2} \tag*{(93)} E P = ( 3 + 4 β 2 ) 2 9 ( 1 − 4 β 2 ) (93)
which could still deviate from unity. For this particular case, then, while the background and the linearly perturbed behavior are indistinguishable from uncoupled dark energy, the deviation from the equivalence principle can still be detected.
Phase-space system for the doubly coupled model including baryons
In this appendix we rewrite the background equations in terms of the variables x x x , y y y introduced in the main text and then identify some of the critical points. Starting with the equation for Ω ~ d m \widetilde{\Omega}_{\rm dm} Ω dm with κ = 1 \kappa= 1 κ = 1 ,
1 = ρ d m 3 f 1 H 2 + ρ b 3 f 1 H 2 + 1 2 ( 1 − 6 Q 2 ) ϕ ′ 2 3 + V ( ϕ ) 3 f 1 H 2 − f 1 , ϕ ϕ ′ f 1 = Ω ~ d m + Ω ~ b + ( 1 − 6 Q 2 ) x 2 + y 2 − 4 β x . (94) \begin{aligned}
1 &= \frac{\rho_{\rm dm}}{3f_1H^2}+\frac{\rho_b}{3f_1H^2}+\frac{1}{2}(1-6Q^2)\frac{\phi'^2}{3}+\frac{V(\phi)}{3f_1H^2}-\frac{f_{1,\phi}\phi'}{f_1} \\
&= \widetilde{\Omega}_{\rm dm}+\widetilde{\Omega}_b+(1-6Q^2)x^2+y^2-4\beta x . \tag*{(94)}
\end{aligned} 1 = 3 f 1 H 2 ρ dm + 3 f 1 H 2 ρ b + 2 1 ( 1 − 6 Q 2 ) 3 ϕ ′2 + 3 f 1 H 2 V ( ϕ ) − f 1 f 1 , ϕ ϕ ′ = Ω dm + Ω b + ( 1 − 6 Q 2 ) x 2 + y 2 − 4 β x . (94)
Introducing a third variable z = 1 H ρ b 3 f 1 z=\frac{1}{H}\sqrt{\frac{\rho_b}{3f_1}} z = H 1 3 f 1 ρ b associated with the baryon density fraction, we get
Ω ~ d m = 1 − ( 1 − 4 β 2 ) x 2 − y 2 + 4 β x − z 2 . (95) \widetilde{\Omega}_{\rm dm}=1-(1-4\beta^2)x^2-y^2+4\beta x-z^2 . \tag*{(95)} Ω dm = 1 − ( 1 − 4 β 2 ) x 2 − y 2 + 4 β x − z 2 . (95)
From
H ′ H = − 3 2 ( 1 + p D E , e f f 3 f 1 H 2 ) , (96) \frac{H'}{H}=-\frac{3}{2}\left(1+\frac{p_{\rm DE,eff}}{3f_1H^2}\right), \tag*{(96)} H H ′ = − 2 3 ( 1 + 3 f 1 H 2 p DE , eff ) , (96)
we include baryons and get
− 2 3 H ′ H = 1 + 4 3 β 2 + ( 1 − 4 β 2 ) [ ( 1 − 4 3 β 2 ) x 2 − y 2 − 4 3 β x ] + 4 3 β β d m [ 1 − ( 1 − 4 β 2 ) x 2 − y 2 + 4 β x − z 2 ] − 2 6 3 β V , ϕ V y 2 . (97) -\frac{2}{3}\frac{H'}{H}=1+\frac{4}{3}\beta^2+(1-4\beta^2)\left[\left(1-\frac{4}{3}\beta^2\right)x^2-y^2-\frac{4}{3}\beta x\right]+\frac{4}{3}\beta\beta_{\rm dm}\left[1-(1-4\beta^2)x^2-y^2+4\beta x-z^2\right]-\frac{2\sqrt{6}}{3}\beta\frac{V_{,\phi}}{V}y^2 . \tag*{(97)} − 3 2 H H ′ = 1 + 3 4 β 2 + ( 1 − 4 β 2 ) [ ( 1 − 3 4 β 2 ) x 2 − y 2 − 3 4 β x ] + 3 4 β β dm [ 1 − ( 1 − 4 β 2 ) x 2 − y 2 + 4 β x − z 2 ] − 3 2 6 β V V , ϕ y 2 . (97)
Putting now V = V 0 e − μ ϕ V=V_0e^{-\mu\phi} V = V 0 e − μ ϕ , we obtain an equation containing only x , y , z x,y,z x , y , z :
− 2 3 H ′ H = 1 + 4 3 β 2 + ( 1 − 4 β 2 ) [ ( 1 − 4 3 β 2 ) x 2 − y 2 − 4 3 β x ] + 4 3 β β d m [ 1 − ( 1 − 4 β 2 ) x 2 − y 2 + 4 β x − z 2 ] + 2 6 3 β μ y 2 . (98) -\frac{2}{3}\frac{H'}{H}=1+\frac{4}{3}\beta^2+(1-4\beta^2)\left[\left(1-\frac{4}{3}\beta^2\right)x^2-y^2-\frac{4}{3}\beta x\right]+\frac{4}{3}\beta\beta_{\rm dm}\left[1-(1-4\beta^2)x^2-y^2+4\beta x-z^2\right]+\frac{2\sqrt{6}}{3}\beta\mu y^2 . \tag*{(98)} − 3 2 H H ′ = 1 + 3 4 β 2 + ( 1 − 4 β 2 ) [ ( 1 − 3 4 β 2 ) x 2 − y 2 − 3 4 β x ] + 3 4 β β dm [ 1 − ( 1 − 4 β 2 ) x 2 − y 2 + 4 β x − z 2 ] + 3 2 6 β μ y 2 . (98)
Continuing with the Klein-Gordon equation:
ϕ ′ ′ + ( 3 + H ′ H ) ϕ ′ + 1 2 f 1 , ϕ f 1 ϕ ′ 2 = − 1 ( 1 − 4 β 2 ) f 1 H 2 [ V , ϕ − 3 f 1 , ϕ f 1 H 2 ( 2 + H ′ H ) + f 2 , ϕ 2 f 2 ρ d m ] , (99) \phi''+\left(3+\frac{H'}{H}\right)\phi'+\frac{1}{2}\frac{f_{1,\phi}}{f_1}\phi'^2=-\frac{1}{(1-4\beta^2)f_1H^2}\left[V_{,\phi}-\frac{3f_{1,\phi}}{f_1}H^2\left(2+\frac{H'}{H}\right)+\frac{f_{2,\phi}}{2f_2}\rho_{\rm dm}\right], \tag*{(99)} ϕ ′′ + ( 3 + H H ′ ) ϕ ′ + 2 1 f 1 f 1 , ϕ ϕ ′2 = − ( 1 − 4 β 2 ) f 1 H 2 1 [ V , ϕ − f 1 3 f 1 , ϕ H 2 ( 2 + H H ′ ) + 2 f 2 f 2 , ϕ ρ dm ] , (99)
which can be written as
x ′ + ( 3 + H ′ H ) x + 2 β x 2 = 3 μ y 2 ( 1 − 4 β 2 ) 6 + 2 1 − 4 β 2 β ( 2 + H ′ H ) + β d m Ω ~ d m 1 − 4 β 2 . (100) x'+\left(3+\frac{H'}{H}\right)x+2\beta x^2=\frac{3\mu y^2}{(1-4\beta^2)\sqrt{6}}+\frac{2}{1-4\beta^2}\beta\left(2+\frac{H'}{H}\right)+\beta_{\rm dm}\frac{\widetilde{\Omega}_{\rm dm}}{1-4\beta^2}. \tag*{(100)} x ′ + ( 3 + H H ′ ) x + 2 β x 2 = ( 1 − 4 β 2 ) 6 3 μ y 2 + 1 − 4 β 2 2 β ( 2 + H H ′ ) + β dm 1 − 4 β 2 Ω dm . (100)
Then, from the definition of y y y and z z z , we also have
y ′ y = − H ′ H + 1 2 V ′ V − 1 2 f 1 ′ f 1 = − H ′ H − 1 2 μ 6 x − 2 β x . (101) \begin{aligned}
\frac{y'}{y}=-\frac{H'}{H}+\frac{1}{2}\frac{V'}{V}-\frac{1}{2}\frac{f_1'}{f_1} \\
&=-\frac{H'}{H}-\frac{1}{2}\mu\sqrt{6}x-2\beta x . \tag*{(101)}
\end{aligned} y y ′ = − H H ′ + 2 1 V V ′ − 2 1 f 1 f 1 ′ = − H H ′ − 2 1 μ 6 x − 2 β x . (101)
z ′ z = − 3 2 − H ′ H − 2 β x . (102) \frac{z'}{z}=-\frac{3}{2}-\frac{H'}{H}-2\beta x . \tag*{(102)} z z ′ = − 2 3 − H H ′ − 2 β x . (102)
So the final system of autonomous equations is:
x ′ = − ( 3 + H ′ H ) x − 2 β x 2 + 6 μ y 2 2 ( 1 − 4 β 2 ) + 2 1 − 4 β 2 β ( 2 + H ′ H ) + β d m Ω ~ d m 1 − 4 β 2 , (103) x'=-\left(3+\frac{H'}{H}\right)x-2\beta x^2+\frac{\sqrt{6}\mu y^2}{2(1-4\beta^2)}+\frac{2}{1-4\beta^2}\beta\left(2+\frac{H'}{H}\right)+\beta_{\rm dm}\frac{\widetilde{\Omega}_{\rm dm}}{1-4\beta^2}, \tag*{(103)} x ′ = − ( 3 + H H ′ ) x − 2 β x 2 + 2 ( 1 − 4 β 2 ) 6 μ y 2 + 1 − 4 β 2 2 β ( 2 + H H ′ ) + β dm 1 − 4 β 2 Ω dm , (103)
y ′ = y ( − H ′ H − 1 2 μ 6 x − 2 β x ) , (104) y'=y\left(-\frac{H'}{H}-\frac{1}{2}\mu\sqrt{6}x-2\beta x\right), \tag*{(104)} y ′ = y ( − H H ′ − 2 1 μ 6 x − 2 β x ) , (104)
z ′ = z ( − 3 2 − H ′ H − 2 β x ) . (105) z'=z\left(-\frac{3}{2}-\frac{H'}{H}-2\beta x\right). \tag*{(105)} z ′ = z ( − 2 3 − H H ′ − 2 β x ) . (105)
For β = 0 \beta=0 β = 0 we get
x ′ = − ( 3 + H ′ H ) x + 6 μ y 2 2 + β d m ( 1 − x 2 − y 2 − z 2 ) , (106) x'=-\left(3+\frac{H'}{H}\right)x+\frac{\sqrt{6}\mu y^2}{2}+\beta_{\rm dm}\left(1-x^2-y^2-z^2\right), \tag*{(106)} x ′ = − ( 3 + H H ′ ) x + 2 6 μ y 2 + β dm ( 1 − x 2 − y 2 − z 2 ) , (106)
y ′ = y ( − H ′ H − 6 2 μ x ) , (107) y'=y\left(-\frac{H'}{H}-\frac{\sqrt{6}}{2}\mu x\right), \tag*{(107)} y ′ = y ( − H H ′ − 2 6 μx ) , (107)
z ′ = z ( − 3 2 − H ′ H ) . (108) z'=z\left(-\frac{3}{2}-\frac{H'}{H}\right). \tag*{(108)} z ′ = z ( − 2 3 − H H ′ ) . (108)
which coincides with the usual equations of CDE, with
H ′ H = − 3 2 ( 1 + x 2 − y 2 ) (109) \frac{H'}{H}=-\frac{3}{2}\left(1+x^2-y^2\right) \tag*{(109)} H H ′ = − 2 3 ( 1 + x 2 − y 2 ) (109)
Now we find the critical points of the system Eqs. Eqs. (103) , (C11) for y = 0 y=0 y = 0 and z = 0 z=0 z = 0 :
x ′ = − ( 3 + H ′ H ) x − 2 β x 2 + 2 β 1 − 4 β 2 ( 2 + H ′ H ) + β d m 1 − ( 1 − 4 β 2 ) x 2 + 4 β x 1 − 4 β 2 , (110) x'=-\left(3+\frac{H'}{H}\right)x-2\beta x^2+\frac{2\beta}{1-4\beta^2}\left(2+\frac{H'}{H}\right)+\beta_{\rm dm}\frac{1-(1-4\beta^2)x^2+4\beta x}{1-4\beta^2}, \tag*{(110)} x ′ = − ( 3 + H H ′ ) x − 2 β x 2 + 1 − 4 β 2 2 β ( 2 + H H ′ ) + β dm 1 − 4 β 2 1 − ( 1 − 4 β 2 ) x 2 + 4 β x , (110)
y ′ = 0 , (111) y'=0, \tag*{(111)} y ′ = 0 , (111)
z ′ = 0 , (112) z'=0, \tag*{(112)} z ′ = 0 , (112)
with
− 2 3 H ′ H = 1 + 4 3 β 2 + ( 1 − 4 β 2 ) [ ( 1 − 4 3 β 2 ) x 2 − 4 3 β x + 4 3 β β d m [ 1 − ( 1 − 4 β 2 ) x 2 + 4 β x ] . (113) -\frac{2}{3}\frac{H'}{H}=1+\frac{4}{3}\beta^2+(1-4\beta^2)[(1-\frac{4}{3}\beta^2)x^2-\frac{4}{3}\beta x+\frac{4}{3}\beta\beta_{\rm dm}[1-(1-4\beta^2)x^2+4\beta x]. \tag*{(113)} − 3 2 H H ′ = 1 + 3 4 β 2 + ( 1 − 4 β 2 ) [( 1 − 3 4 β 2 ) x 2 − 3 4 β x + 3 4 β β dm [ 1 − ( 1 − 4 β 2 ) x 2 + 4 β x ] . (113)
Then we can simplify the x ′ x' x ′ equation and obtain
x ′ = 1 2 ( ( 4 β 2 − 1 ) x 2 + 4 β x + 1 ) ( 2 β + β d m ( 4 β x + 2 ) + ( 4 β 2 − 3 ) x ) . (114) x'=\frac{1}{2}\left((4\beta^2-1)x^2+4\beta x+1\right)\left(2\beta+\beta_{\rm dm}(4\beta x+2)+(4\beta^2-3)x\right). \tag*{(114)} x ′ = 2 1 ( ( 4 β 2 − 1 ) x 2 + 4 β x + 1 ) ( 2 β + β dm ( 4 β x + 2 ) + ( 4 β 2 − 3 ) x ) . (114)
The y = 0 y=0 y = 0 critical points are then given by solving x ′ = 0 x'=0 x ′ = 0 :
x ± = 1 ± 1 − 2 β , (115) x_{\pm}=\frac{1}{\pm1-2\beta}, \tag*{(115)} x ± = ± 1 − 2 β 1 , (115)
x i = 2 ( β + β d m ) 3 − 4 β ( β + β d m ) . (116) x_i=\frac{2(\beta+\beta_{\rm dm})}{3-4\beta(\beta+\beta_{\rm dm})}. \tag*{(116)} x i = 3 − 4 β ( β + β dm ) 2 ( β + β dm ) . (116)
Since these critical points exist for y = 0 y=0 y = 0 , they are independent of the particular choice of the potential. On the two points x ± x_{\pm} x ± , Ω ~ d m = 0 \widetilde{\Omega}_{\rm dm}=0 Ω dm = 0 , so they cannot be valid early-time solutions. Instead, the last point yields a finite amount of dark matter and early dark energy, and should be taken as initial condition, as we have shown in the previous appendix. A detailed analysis of the full phase space and of the stability of the critical points is left to future work.
Details on the Fisher Matrix Projection
In this appendix, we present details on the Fisher matrix projection. More specifically, we project the model-independent constraints on the parameters { E P , h d , d } \{E_{\rm P},h_d,d\} { E P , h d , d } , obtained using a methodology similar to that in Refs. [41 , 42 ] (with each parameter left free to vary at each redshift bin), onto a set of constant parameters: { β , β d m , Ω m 0 } \{\beta,\beta_{\rm dm},\Omega_{\rm m0}\} { β , β dm , Ω m0 } . When all three parameters are varied, we obtain constraints on the doubly coupled dark energy model, while in the case β = 0 \beta=0 β = 0 , the projection corresponds to the coupled dark energy model. This can be done by transforming the original Fisher matrix, F o l d \mathcal{F}_{\rm old} F old , which contains only the information on { E P , h d , d } \{E_P,h_d,d\} { E P , h d , d } at different redshift bins, into a new 3 × 3 3\times3 3 × 3 Fisher matrix, F n e w \mathcal{F}_{\rm new} F new , from whose inverse we obtain the marginalised errors on { β , β d m , Ω m 0 } \{\beta,\beta_{\rm dm},\Omega_{\rm m0}\} { β , β dm , Ω m0 } , as follows:
F n e w = J T F o l d J , (117) \mathcal{F}_{\rm new}=J^{T}\mathcal{F}_{\rm old}J, \tag*{(117)} F new = J T F old J , (117)
with J J J the Jacobian matrix and J T J^T J T its transpose. For simplicity, hereafter we define that, in the i-th redshift bin, E P ( N i ) ≡ E P , i E_P(N_i)\equiv E_{P,i} E P ( N i ) ≡ E P , i , h d ( N i ) ≡ h d , i h_d(N_i)\equiv h_{d,i} h d ( N i ) ≡ h d , i and d ( N i ) ≡ d i d(N_i)\equiv d_i d ( N i ) ≡ d i . The E P ( N ) E_P(N) E P ( N ) expressions for our specific models are given by Eqs. (48) and (51) , while those for d ( N ) d(N) d ( N ) and h d ( N ) h_d(N) h d ( N ) are given by
d ( N ) = L A ( N ) L A r ( N ) = ∫ N 0 e − N ′ d N ′ E ( N ′ ) ∫ N 0 e − N ′ d N ′ E r ( N ′ ) , (118) d(N)=\frac{L_A(N)}{L_{A_r}(N)}=\frac{\int_N^0\frac{e^{-N'}\,dN'}{E(N')}}{\int_N^0\frac{e^{-N'}\,dN'}{E_r(N')}}, \tag*{(118)} d ( N ) = L A r ( N ) L A ( N ) = ∫ N 0 E r ( N ′ ) e − N ′ d N ′ ∫ N 0 E ( N ′ ) e − N ′ d N ′ , (118)
h d ( N ) = h ( N ) d ( N ) = E ( N ) E r ( N ) L A ( N ) L A r ( N ) = E ( N ) ∫ N 0 e − N ′ d N ′ E ( N ′ ) E r ( N ) ∫ N 0 e − N ′ d N ′ E r ( N ′ ) , (119) h_d(N)=h(N)d(N)=\frac{E(N)}{E_r(N)}\frac{L_A(N)}{L_{A_r}(N)}=\frac{E(N)\int_N^0\frac{e^{-N'}\,dN'}{E(N')}}{E_r(N)\int_N^0\frac{e^{-N'}\,dN'}{E_r(N')}}, \tag*{(119)} h d ( N ) = h ( N ) d ( N ) = E r ( N ) E ( N ) L A r ( N ) L A ( N ) = E r ( N ) ∫ N 0 E r ( N ′ ) e − N ′ d N ′ E ( N ) ∫ N 0 E ( N ′ ) e − N ′ d N ′ , (119)
where E r E_r E r is the expansion rate for the reference cosmology, i.e. the Λ \Lambda Λ CDM model. The expansion rate E E E , on the other hand, depend on β \beta β and β d m \beta_{\rm dm} β dm through x x x and y y y in Eqs. (52) and (53) , with the exact solution given by Eq. (78) . Moreover, in the CDE case, E ( N ) E(N) E ( N ) simply reduces to
E ( N ) = e − 3 / 2 ∫ 0 N ( x ( N ′ ) 2 − y ( N ′ ) 2 + 1 ) d N ′ , (120) E(N)=e^{-3/2\int_0^N(x(N')^2-y(N')^2+1)dN'}, \tag*{(120)} E ( N ) = e − 3/2 ∫ 0 N ( x ( N ′ ) 2 − y ( N ′ ) 2 + 1 ) d N ′ , (120)
with the x x x and y y y given by Eqs. (49) and (50) . This expansion is well justified, as the relative difference from E r ( N ) E_r(N) E r ( N ) remains below 0.4 % 0.4\% 0.4% over the low-redshift range considered in our forecast, as shown in Table VIII .
Table VIII. Comparison of the reference expansion rate E r E_r E r and the expanded expression for E E E at different redshifts. The last column shows the relative difference between the two expressions.
The Jacobian elements J ( i ) J^{(i)} J ( i ) within each redshift bin η i \eta_i η i are given by:
J ( i ) = ( ∂ h d , i ∂ β ∂ h d , i ∂ β d m ∂ h d , i ∂ Ω m 0 ∂ d i ∂ β ∂ d i ∂ β d m ∂ d i ∂ Ω m 0 ∂ E P , i ∂ β ∂ E P , i ∂ β d m ∂ E P , i ∂ Ω m 0 ) , (121) J^{(i)}=
\begin{pmatrix}
\frac{\partial h_{d,i}}{\partial\beta} & \frac{\partial h_{d,i}}{\partial\beta_{\rm dm}} & \frac{\partial h_{d,i}}{\partial\Omega_{\rm m0}}\\
\frac{\partial d_i}{\partial\beta} & \frac{\partial d_i}{\partial\beta_{\rm dm}} & \frac{\partial d_i}{\partial\Omega_{\rm m0}}\\
\frac{\partial E_{P,i}}{\partial\beta} & \frac{\partial E_{P,i}}{\partial\beta_{\rm dm}} & \frac{\partial E_{P,i}}{\partial\Omega_{\rm m0}}
\end{pmatrix}, \tag*{(121)} J ( i ) = ∂ β ∂ h d , i ∂ β ∂ d i ∂ β ∂ E P , i ∂ β dm ∂ h d , i ∂ β dm ∂ d i ∂ β dm ∂ E P , i ∂ Ω m0 ∂ h d , i ∂ Ω m0 ∂ d i ∂ Ω m0 ∂ E P , i , (121)
with each entry evaluated at the fiducial values, e.g.:
∂ E P ( N i ) ∂ β ∣ f i d = ∂ ∂ β [ E P ( N i ; β , β d m , Ω m 0 ) ] ∣ N = N i ; β = β f i d , β d m = β d m , f i d , Ω m 0 = Ω m 0 , f i d . (122) \left.\frac{\partial E_P(N_i)}{\partial\beta}\right|_{\rm fid}
=
\left.\frac{\partial}{\partial\beta}\left[E_P(N_i;\beta,\beta_{\rm dm},\Omega_{\rm m0})\right]\right|_{N=N_i;\,\beta=\beta_{\rm fid},\,\beta_{\rm dm}=\beta_{\rm dm,fid},\,\Omega_{\rm m0}=\Omega_{\rm m0,fid}}. \tag*{(122)} ∂ β ∂ E P ( N i ) fid = ∂ β ∂ [ E P ( N i ; β , β dm , Ω m0 ) ] N = N i ; β = β fid , β dm = β dm , fid , Ω m0 = Ω m0 , fid . (122)
In this paper we set: β f i d = 0.1 \beta_{\rm fid}=0.1 β fid = 0.1 , β d m , f i d = 0.2 \beta_{\rm dm,fid}=0.2 β dm , fid = 0.2 , and Ω m 0 , f i d = 0.315 \Omega_{\rm m0,fid}=0.315 Ω m0 , fid = 0.315 .
For the full Jacobian J J J including all the number of N N N redshift bins, the blocks are stacked as:
J = ( J ( 1 ) J ( 2 ) ⋮ J ( N ) ) . (123) J=
\begin{pmatrix}
J^{(1)}\\
J^{(2)}\\
\vdots\\
J^{(N)}
\end{pmatrix}. \tag*{(123)} J = J ( 1 ) J ( 2 ) ⋮ J ( N ) . (123)
Estimates of the relative error on the luminosity distance
We give now a rough estimate of the current uncertainty on the luminosity distance d d d that we use as prior in the main text. For N N N independent supernovae with identical Gaussian apparent-magnitude uncertainty δ m \delta m δ m and perfectly known absolute magnitude, the relative error on d d d is
σ d d = ln 10 5 δ m N . (124) \frac{\sigma_d}{d} = \frac{\ln10}{5}\frac{\delta m}{\sqrt{N}}. \tag*{(124)} d σ d = 5 ln 10 N δ m . (124)
Using the current Pantheon+ catalog [62 ] , we show the resulting uncertainty on d d d in Table IX
Table IX. Pantheon+ supernova luminosity-distance uncertainties.
These results are just a rough estimate, since we assumed a constant δ m \delta m δ m and ignored correlated systematics, peculiar velocities, and intrinsic dispersion. However, they also show that within the SKA redshift range, our prior of 1% is justified, and will be soon improved by further data.
References