Gravitational redshift constraints on the effective theory of interacting dark energy

Upcoming galaxy surveys provide the necessary sensitivity to measure gravitational redshift, a general relativistic effect that generates a dipole in galaxy clustering data when correlating two distinct populations of galaxies. Here, we study the constraining power of gravitational redshift within the framework of the effective theory of interacting dark energy. This formalism describes linear cosmological perturbations in scalar-tensor theories of gravity with a limited number of free functions, and allows each particle species to be coupled differently to the gravitational sector. In this work, we focus on Horndeski theories with a non-minimal coupling of dark matter to the scalar degree of freedom, yielding a breaking of the weak equivalence principle for this cosmic component, a scenario that is yet untested. We show that the dipole generated by gravitational redshift significantly breaks degeneracies and tightens the constraints on the parameters of the effective theory compared to the standard redshift-space distortion analysis solely based on the even multipoles in the galaxy correlation function, with an improvement of up to ∼ 50% for populations with a galaxy bias difference equal to 1. We make the Python package EF-TIGRE (Effective Field Theory of Interacting dark energy with Gravitational REdshift) developed for this work publicly available (https://github.com/Mik3M4n/EF-TIGRE).

1 Introduction In recent years, a considerable amount of both theoretical and observational work in cosmology has focused on testing gravity on scales larger than individual galaxies, in the quest to unveil the unknown mechanism driving the accelerated expansion of the Universe.A powerful observational tool to achieve this consists in employing maps of galaxies in redshift space, which are directly sensitive to the peculiar velocities of the galaxies.Since different theories of gravity generate different velocity fields, comparing these measurements with theoretical predictions provides an efficient way of testing models beyond General Relativity (GR).However, the difficulty of such studies lies in the huge number of theories that have been proposed, making it impractical to compare them against the data one by one.
A powerful approach is provided by the effective field theory (EFT) of dark energy [1][2][3][4][5][6][7][8][9], which encompasses a large class of theories in a unified formalism, describing a rich phenomenology beyond GR. 1 This formalism is based on a model-independent construction at the level of the action, providing a description of linear cosmological perturbations in scalar-tensor theories of gravity around a homogeneous and isotropic background.In particular, it covers Horndeski models [13], which constitute the most general class of Lorentz-invariant scalartensor theories with second-order equations of motion (this comprises e.g.f (R) gravity [14], quintessence [15,16], Brans-Dicke models [17], kinetic braiding [18], Galileons/generalised Galileons [19][20][21]).In this paper, we will focus on models of this class, where four functions of time fully describe the perturbations [3,6,8]: the kineticity α K , the braiding α B , the Planck-mass run rate α M and the tensor speed excess α T . 2 Observational as well as theoretical constraints have been put on these free functions, combining the Cosmic Microwave Background (CMB), large-scale structure data and gravitational-wave detections [27][28][29][30][31][32][33][34][35][36][37].However, at the moment only the parameter α T is tightly restricted to its GR value of zero by the gravitational-wave event GW170817 and its electromagnetic counterpart, GRB170817A [38].Even if this already restricts the allowed space of theories [27][28][29][30], a wide range of allowed models remains, and data from next-generation surveys are necessary to put stringent constraints on the remaining free functions [39][40][41]. 3 Moreover, the aforementioned constraints are based on the assumption of a universal coupling of the gravitational sector to the various matter species.This means that the validity of the weak equivalence principle (WEP) is preserved.However, the universality of fifth force effects, mediated by additional degrees of freedom associated to deviations from GR, is yet untested.More specifically, fifth force effects on baryons and photons are strongly constrained (see e.g.ref. [47][48][49][50][51]), while those on cold dark matter (CDM) could be significant and indeed occur in various theories [52][53][54][55][56][57][58][59][60][61], leading to a breaking of the WEP.To include such theories with distinct couplings, the effective treatment can be broadened to an extension known as effective theory of interacting dark energy.This framework has been introduced in ref. [62] in the context of Horndeski theories and further developed in refs.[26,63] to include more general classes of scalar-tensor theories.Within this formalism, ref. [39] considered a WEP-breaking scenario where photons and baryons have standard couplings to the metric, whereas CDM has a non-minimal coupling described by a single additional free function γ c .The capability of future galaxy surveys to constrain the parameters of this framework was forecasted through a Fisher matrix approach based on a combination of the galaxy power spectrum in redshift space, the tomographic weak-lensing shear power spectrum, and the correlation between the integrated Sachs-Wolfe effect and the galaxy distribution.This analysis highlighted strong degeneracies between the effects of a fifth force and the other modifications.
As shown in ref. [64] using a phenomenological approach, these degeneracies are due to the fact that the growth of structure is affected in a similar manner by modifications in the Poisson and Euler equations.Both these equations determine how matter falls into a gravitational potential and consequently how it clusters, so their respective modifications cannot be distinguished through redshift-space distortion (RSD) analyses, which only probe the growth of structure.Gravitational lensing does not help disentangling the two effects, since it does not probe the spatial distortion of the geometry and the time distortion separately [65].Luckily, ref. [64] has demonstrated that a measurement of the distortion of time through the effect of gravitational redshift on cosmological scales can break the degeneracies, since it provides a way of testing both the Euler equation (by comparing the distortion of time with the galaxy velocity), and the relation between the temporal and spatial distortion (by comparing the distortion of time with gravitational lensing [66][67][68]).
Gravitational redshift is an effect predicted by General Relativity [69] resulting in a change of wavelength when photons travel across gravitational potentials.This effect has been measured in the laboratory [70], from solar observations [71], on astrophysical scales from stacked galaxy clusters [72][73][74][75][76][77], as well as in the non-linear regime of large-scale structure [78,79].Here, we focus on cosmological scales described within linear perturbation theory and consider the two-point correlation function of galaxies as an observable.Together with other relativistic corrections, gravitational redshift generates a dipolar modulation in the correlation function [80][81][82] (or similarly an imaginary part in the power spectrum [83,84]), which only appears by cross-correlating two populations of galaxies with different clustering properties, for example bright and faint galaxies.Combining the relativistic dipole with the even multipoles generated by RSDmakes it possible to isolate the contribution due to gravitational redshift [66].While the even multipoles of the correlation function have been successfully measured with various surveys, see e.g.[85][86][87][88], the dipole is too small to be detected by current data in the linear regime [89].Detailed forecasts have however shown that it will be measurable by upcoming galaxy surveys, such as the Dark Energy Spectroscopic Instrument (DESI) [90,91] and the Square Kilometre Array (SKA) [92,93].In particular, SKA is expected to observe close to a billion galaxies up to z = 2, allowing for a detection of the dipole with signal-to-noise of 80 (see section 5.1.2).
In this work, we investigate the constraining power of gravitational redshift within the framework of the effective theory of interacting dark energy established in refs.[39,62].With this approach, we provide forecasts on constraints on fundamental parameters that enter at the level of the Lagrangian and cover a wide range of modified gravity theories, more specifically all Horndeski theories with an additional breaking of the WEP for CDM.Using a Monte Carlo Markov Chain (MCMC) analysis, we forecast the expected constraints from a survey like SKA phase 2 (SKA2) [94].We show that, without gravitational redshift, the parameters are strongly degenerated, leading to wide contours with two branches in most 2dimensional projections.Such degeneracies leading to highly non-Gaussian posteriors cannot be captured by a Fisher analysis and require an MCMC approach.We show that the inclusion of gravitational redshift significantly alleviates the degeneracies and tightens the constraints on all parameters, reducing the 1 σ error bars by up to ∼ 50 percent when the bias difference between the two populations is assumed to be 1, which is a reasonable value for current and upcoming galaxy surveys [89].This is truly remarkable, since gravitational redshift can be extracted from the same data as the growth rate of structure by measuring a dipole from the cross-correlation of galaxies in addition to the standard even multipoles.
This paper is structured as follows: In section 2, we revise the effective theory of interacting dark energy developed in refs.[39,62], while highlighting differences to the original work in our implementation of the formalism.The expressions for the observables used for our analysis, in particular the galaxy clustering dipole arising from the general relativistic time distortion, are given in section 3.In section 4, we discuss the survey specifications and our numerical implementation of the effective theory approach in the EF-TIGRE (Effective Field Theory of Interacting dark energy with Gravitational REdshift) Python code.Finally, we present the results on the constraining power of gravitational redshift in section 5, and we conclude in section 6.Moreover, in appendix A we state the relation between our model parameters and commonly employed, phenomenologically-motivated parameters, while appendix B contains supplementary plots.

Effective theory of interacting dark energy
In this work, we adopt the effective theory of interacting dark energy established in refs.[39,62].This comprehensive formalism allows to describe a wide class of modified gravity models encompassing all Horndeski theories, while allowing for different couplings of the matter sector to gravity due to interactions in the dark sector.In particular, this enables us to consider a breaking of the WEP for the unknown CDM component.We summarize our approach in the following, referring to refs.[39,62] for details and highlighting the differences with respect to the original setup.

Gravitational and matter sector
As in refs.[39,62], we assume that the gravitational sector is described by a four-dimensional metric g µν and a scalar field ϕ, with an action corresponding to the class of Horndeski theories [13].The action can be written in the unitary gauge, where the constant-time hypersurfaces coincide with the uniform scalar field ones.In ADM form [95], the metric can be written as 4ds where N is the lapse function, N i is the shift and h ij the 3-dimensional spatial metric.The gravitational action can then be expressed in the general form with K ij being the extrinsic curvature tensor and R ij the Ricci tensor associated to the constant time hypersurfaces.The contribution from the matter action can be generically written as with each of the N S particle species being minimally coupled to the gravitational sector by a generally different effective metric ǧ(I) µν .In this work, we consider contributions to the total energy density arising from baryons, CDM and radiation.We assume that the coupling of CDM to the gravitational sector is different from the one of baryons and radiation, leading to a breaking of the WEP.Thus, we cannot identify a global Jordan frame, but we choose to work in the Jordan frame of baryons and radiation, defining the corresponding metric as the gravitational metric g µν .In this frame, CDM is allowed to be non-minimally coupled to g µν .We denote baryons, CDM and radiation by the subscripts b, c and r respectively.

Background
We consider a homogeneous and isotropic background described by the Friedmann-Lemaître-Robertson-Walker (FLRW) metric where a(t) is the scale factor.In the following, a bar denotes the 0-th order quantities and we assume spatial flatness.The background evolution equations can be obtained by taking variations of the homogeneous action Sg + Sm .We use the logarithm of the scale factor ln(a) as a time variable, denoting the corresponding derivatives with a prime, and write the Friedmann equations as (see eqs. (2.12)-(2.13) of ref. [62]) (2.6) Here, ρ denotes the energy density and p the pressure, and we assume that baryons and CDM are pressureless, p b = p c = 0. H is the Hubble parameter and M denotes the effective Planck mass, which in this framework is a function of time.Its time evolution is determined by the function α M , defined by The equations governing the energy densities of baryons, CDM and radiation follow from the invariance of the matter action under arbitrary diffeomorphisms.These are given in eqs.(4.2)-(4.4) of ref. [62] and read ) ρ′ r = −4ρ r . (2.10) Hence, baryons and photons evolve in a standard way, as perfect fluids with pressure p b = 0 and p r = ρ r /3 respectively, while the evolution of CDM is modified by its non-minimal coupling to gravity, encoded in the free function γ c . 5 It is convenient to introduce the dimensionless density parameters Ω X ≡ ρX 3M 2 H 2 , and rewrite the evolution equations for ρb , ρc and ρr as ) ) In the above equations, Ω DE is defined via eq.(2.5) by The Hubble parameter can be calculated from the solutions of eq.(2.12) to (2.14) by writing eq.(2.6) as and integrating this expression while fixing as initial condition the value H 0 at present time (see section 4 for a discussion about the value of H 0 ), (2.16)

Perturbations
To compute the equations of motions for linear perturbations, the gravitational Lagrangian in eq.(2.2) must be expanded up to second order.The coefficients of the second-order expansion can be conveniently factorized in three dimensionless combinations denoted by α B , α K and α T (see section 2.3 of ref. [62] for details).The most general quadratic action leading to second-order equations of motion can consequently be written as (2.17) Here, δ 2 denotes the second-order term in a perturbative expansion.The gravity modifications entering eq.(2.17) are thus parameterized by four time-dependent functions: the tensor speed excess α T , the kineticity α K , the braiding α B , 6 and the effective Planck mass M , whose evolution is encoded in the parameter α M defined in eq.(2.7).The parameter α T is directly related to the speed of propagation of tensor modes c T ≡ 1 + α T , which was constrained to match the speed of light up to ∼ 10 −15 by the multi-messenger observation of the binary neutron star merger GW170817 in conjunction with the gamma-ray burst GRB 170817A in 2017 [38].Therefore, we set α T = 0 throughout our analysis. 7We are thus left with the three free functions α M , α K and α B , which take the values α K = α B = α M = 0 in the special case of GR.
In order to write the perturbations equations, it is convenient to adopt the Newtonian gauge for the linearly perturbed FLRW metric, (2.18) We work in the quasi-static approximation, following the approach of ref. [99], where a hierarchy of terms is introduced based on the assumption that time derivatives are smaller than spatial ones.With these assumptions, the action in eq.(2.2) leads to the following Euler equations for baryons and CDM (see eqs. (4.16) and (4.18) in ref. [62], where Φ and Ψ are interchanged),8 Here, V b,c are the velocity potentials of baryons and CDM in Fourier space, related to the velocity fields by V b,c = −ik/k V b,c , and the quantities c 2 s α and µ are defined below.The Poisson equation reads (see eq. (4.19) in ref. [62]) i.e. the Newton constant G is replaced by an effective Newton constant, Gµ, which depends on the gravity modifications and the non-minimal coupling γ c (see eq. (2.25) below).Combining eq.(2.21) with the continuity equation and the Euler equations gives rise to second-order evolution equations for the density (see eqs. (4.23) and (4.24) in ref. [62]), Here, Ω m and δ m refer to the total matter contribution, and ω c ≡ Ω c /Ω m and b c ≡ δ c /δ m denote the fraction of CDM at the level of the background and perturbations respectively.The Euler, evolution and Poisson equations depend on the function ζ given in eq.(2.6), and on the function µ defined as Moreover, the evolution equations contain the quantity c 2 s α, involving the speed of sound of scalar perturbations c s (given by eq.(2.39) of ref. [62]), The quantity α, defined in eq.(2.36) of ref. [62], contains a combination of the functions α K , α B and α D,c , where the latter arises from a disformal coupling of CDM to the metric g µν .However, we note that the combination c 2 s α in eq.(2.26) does not depend on α K and α D,c .This is a consequence of the quasi-static approximation and implies that the kineticity α K remains unconstrained when adopting this limit. 9Finally, the absence of ghost and gradient instabilities requires c 2 s α ≥ 0, and we note that a major advantage of the EFT formalism is that these theoretical constraints can be easily incorporated in the analysis, leading to a large restriction of the parameter space [37].We will include this condition in our analysis.
In the following, we will investigate the potential of upcoming galaxy surveys to constrain the functions α M , α B , γ c and w DE , highlighting the fundamental role played by gravitational redshift.

Galaxy clustering observables
The theoretical framework presented in the previous section can be constrained through observations of the distribution of galaxies across the sky.Galaxy surveys provide measurements of the number of galaxies N as a function of redshift z and angular position n in the sky.We can use these observations to construct the number count fluctuations observable, where N (z) is the average number of galaxies at redshift z.The quantity ∆ can be computed at linear order in perturbation theory [100][101][102][103][104], yielding10 where Here, a dot denotes derivatives with respect to conformal time, r is the comoving distance and H is the conformal Hubble parameter.Eqs.(3.2)-(3.4)are valid in any metric theory of gravity and only rely on the assumption that photons travel along null geodesics of the metric g µν .Eq. (3.3) contains the two dominant contributions to ∆, which we refer to as standard : the first one arises from the matter density, related to the galaxy density by the linear bias b, whereas the second one encodes the RSD due to the peculiar velocity field of galaxies [108].On top of these terms, ∆ is also affected by relativistic corrections collected in eq.(3.4): the first term encodes gravitational redshift, which changes the apparent size of a redshift bin due to the time distortion affecting galaxies situated inside of a potential well.We then have two Doppler terms depending on the velocity V and its time derivative. 11ere, we note that s is the magnification bias (see e.g.refs.[80,101]) arising from magnitude cuts imposed in the sample selection.Note that here we neglect the impact of the evolution bias, whose contribution is expected to be subdominant [111]. 12he standard approach to analyze galaxy clustering data is to apply summary statistics.In this work, we consider the galaxy two-point correlation function, 13 ξ ≡ ⟨∆(n, z)∆(n ′ , z ′ )⟩, which can be expanded in multipoles, Here, µ is the cosine angle between the line of sight and the separation between the correlated galaxies (see figure 2 in ref. [82] for an illustration), d is the comoving distance between the galaxies and P ℓ (µ) denotes the Legendre polynomial of order ℓ.
In the distant-observer regime, the standard contributions (density fluctuations and RSD) generate three even multipoles: a monopole (ℓ = 0), a quadrupole (ℓ = 2) and a hexadecapole (ℓ = 4), which were measured by various surveys, see e.g.ref. [88].These are given by ) where f = d ln(δ m )/d ln(a) is the matter growth rate, and we have defined Here, P δδ (k, z) is the matter power spectrum, which can be computed by solving eqs.(2.22)-(2.23).
The contribution to the even multipoles arising from the general relativistic effects in ∆ rel is suppressed by a factor of (H/k) 2 and is therefore negligible on sub-horizon scales.On the other hand, the cross-correlation of ∆ rel with the standard terms breaks the symmetry of the two-point function ξ [80,81,84], generating a dipole moment. 14To extract this dipole from galaxy surveys, it is necessary to cross-correlate two differently biased populations of galaxies, for example a bright (B) and a faint (F) sample.Employing the Poisson equation (2.21) and the Euler equations (2.19)-(2.20),we obtain the following expression for the dipole, Here, x c is the fraction of CDM mass in a galaxy (which in general can differ from the background CDM fraction ω c as well as from the fraction b c ≡ δ c /δ m at the level of perturbations), and The first line in eq.(3.10) arises from the first three terms of ∆ rel in eq.(3.4), which are related to the density through the modified Poisson equation (2.21) and the Euler equation for CDM (2.20).These terms vanish if γ c = 0, i.e. if CDM is minimally coupled to the scalar field, thus obeying the Euler equation.Their combination is weighted by a factor x c due to the fact that the WEP is only broken for CDM, but not for baryons.The remaining Doppler contributions in eq.(3.4) lead to the second and third line of eq.(3.10).The fourth line contains the wide-angle effect, a contamination to the dipole from RSD due to the fact that the sky is not flat [80].Even though this term is suppressed by d/r, it is of the same order as the contributions in the first three lines, since it is proportional to µ 2 (d), which is enhanced by a factor k/H with respect to ν 1 (d).We note that the dipole in eq.(3.10), including the wide-angle effect, exactly vanishes for a single galaxy population.This would not be the case if we had chosen to align the line-of-sight direction along one of the correlated galaxies instead of along the median direction between the galaxies.In this case, one would have a remaining "large-angle effect" even for a single population of galaxies, as described in refs.[89,113]). 15he galaxy velocity V appearing in eqs.(3.3) and (3.4) is a weighted average of the velocities of baryons and CDM [39], (3.12) The velocity potentials V b,c can be related to the density fields δ b,c via where f b,c (a) ≡ d ln(δ b,c )/d ln(a) is the growth rate of the respective species.Eq. (3.13) corresponds to the linear-order continuity equation, remaining valid for all particle species in the effective theory of interacting dark energy (see eqs. (4.15) and (4.17) in ref. [62]).The galaxy velocity potential then becomes where the growth rate of galaxies f is related to that of baryons and CDM through In the following sections, we forecast the constraints on the effective theory of interacting dark energy arising from measurements of the even multipoles of ξ, and we compare them with the constraints obtained when the dipole is added to the analysis.

Survey specifications and numerical implementation
The effective theory framework and the galaxy clustering observables presented in the previous sections form the foundation for our work.Here, we discuss additional preliminaries necessary for the numerical analysis with our EF-TIGRE Python package.We define the survey specifications and choice of time evolution for the parameters we constrain.Moreover, we specify the set of observables and the fiducial models used in our analysis, and outline further details of our numerical implementation.

Survey specifications
We consider a survey similar to SKA2, covering 30,000 square degrees and ranging from z = 0 to z = 1.5.We divide this redshift range into 15 bins of size ∆z = 0.1 and in each of these bins we consider the survey specifications given in Table 3 in ref. [94]. 16We split the sample of galaxies into two populations with the same number of galaxies and model the bias of each population with the fitting function from ref. [94], where P denotes either the bright (P = B) or faint (P = F) population.We vary the four parameters b 1,B , b 1,F , b 2,B and b 2,F , around fiducial values b 1,P = 0.554 and b 2,P = 0.783 following ref.[94].We fix ∆b = 1, which is consistent with the bias difference obtained from the measurement of BOSS [89].This quantity could in principle be different for SKA2, but once the data become available, various possibilities of splitting the galaxy sample can be explored in order to boost the bias difference [91,115], and thus the dipole signal from eq. (3.10).The magnification biases s B and s F will be measurable from the average number of galaxies within the SKA2 survey, and for the purpose of our forecast we fix their values according to appendix B of ref. [64].

Choices of parameterizations and time evolution of the free functions
In addition to the galaxy and magnification bias parameters, our observables depend on four free functions of time, which are the main focus of our work: the braiding α B , the running of the Planck mass α M , the parameter γ c encoding a non-minimal coupling of CDM to the metric g µν , and the dark energy equation of state parameter w DE .In order to extract constraints on these functions, it is necessary to either assume a parameterized time evolution for them, or to resort to non-parametric models, typically consisting in constraining the free functions in redshift bins [116].Due to the large dimensionality of the parameters space and the strong correlations between the free functions, in this work we choose to adopt the first strategy and assume a specific time evolution.Hence, we assume that w DE is constant in time, whereas α B , α M and γ c decay proportionally to the dark energy density when going backwards in cosmic history, in line with much previous literature (see e.g.refs.[35,39,40,64,117]).This time parameterization is motivated by the fact that the scalar field drives the accelerated expansion of the Universe, and that its impact should therefore be negligible well before acceleration started.In practice, for each parameter p ∈ {α B , α M , γ c }, we write where the superscript ΛCDM denotes quantities evaluated in ΛCDM, i.e. explicitly Note that here we parameterize the time evolution using the fiducial ΛCDM evolution for simplicity.Given that the background evolution is expected to be close to ΛCDM and that the choice of time parameterization involves a degree of arbitrariness, there is no compelling reason to introduce a more intricate evolution.With these assumptions, the degrees of freedom reduce to three free parameters α B,0 , α M,0 and γ c,0 , denoting the present-time values of the functions.The final constraints on these parameters are sensitive to our choice of time parameterization, but we will argue in section 5 that this does not affect the main message of our work, i.e. that the inclusion of gravitational redshift breaks several degeneracies among the parameters.Our choice of time parameterization is the same as in ref. [39] for α M and α B , but it differs concerning γ c .In ref. [39], the authors assumed a constant value for the combination ) and provided constraints on this quantity. 17Instead, here we assume that also γ c evolves according to eq. (4.2), thus ensuring that all free functions entering the evolution equations for the perturbations, eqs.(2.22)-(2.23),share the same time dependence.This is important, since assuming distinct time evolutions could artificially break the degeneracies between parameters, thus enhancing the constraining power of our observables.Our choice of constraining γ c instead of β γ also simplifies the numerical implementation in comparison to ref. [39], as the expression for c 2 s α in eq.(2.26), combined with ζ in eq.(2.6), is directly given in terms of the free functions.To summarize, the relevant degrees of freedom are α B,0 , α M,0 , γ c,0 and w DE , which in the following we refer to as EFT parameters.

Observables
We consider as our observables all possible auto-and cross-correlations between the bright and faint galaxy samples, which we expand in multipoles.This yields a total of three monopoles and three quadrupoles (bright-bright, faint-faint and bright-faint).Since the hexadecapole, given in eq.(3.8), is independent of population-specific biases, no additional information is obtained by considering various samples, so we only include the hexadecapole of the total galaxy population.Thus, we have in total seven even multipoles, which are evaluated at each redshift bin.We compare the constraints obtained from these even multipoles with an analysis where we also include the dipole in the cross-correlation between the bright and faint populations, which includes the effect of gravitational redshift.Note that the dipole is strictly zero for two populations with equal galaxy and magnification bias, as can be seen from eq. (3.10), and therefore no bright-bright or faint-faint dipole exists.We fix the minimum separation between galaxy pairs to d min = 20 Mpc/h, such that the impact of non-linearities is negligible [111].Moreover, after testing that going beyond this range does not improve the resulting constraints, we fix the maximum separation to d max = 160 Mpc/h.
Our choice of time parameterization for the modifications ensures that we recover GR at high redshift, well before the cosmic acceleration started.As a consequence, the constraints from the CMB on the parameters Ω m , Ω b , Ω r at early time as well as on the primordial parameters A s and n s remain valid.Hence, we fix the values of these parameters according to the Planck values [32] and start solving the background equations at recombination.
On the other hand, the dark energy equation of state parameter w DE affects the latetime evolution of the Universe and can differ from -1 in Horndeski theories.We therefore keep it as a free parameter and we cannot directly set a Planck prior on it, as the Planck constraints on this parameter were obtained assuming the validity of GR at late times.We thus introduce an indirect constraint arising from the distance to last scattering, which is inferred in a model-independent manner from the position of the peaks in the CMB angular power spectrum.More precisely, we consider the Planck measurement of the comoving angular diameter distance D M at recombination, 18 which is given by the ratio between the sound horizon at recombination, r * , and the angular scale of the first peak of the CMB, θ * .We include the quantity D M as an independent observable in addition to the multipoles, with a central value and variance determined by Planck [32].The expected value is obtained by integrating the inverse of the Hubble parameter from today up to the redshift z * of recombination, Note that this implies that, unlike for the constraints obtained in ref. [39], we need to include radiation in the background evolution, as it becomes relevant at high redshift.On the other hand, since the multipoles are only observed at late time, it is sufficient to solve the perturbation equations starting at z in = 10, where we can neglect the impact of radiation.Lastly, we fix the value of the Hubble constant today, H 0 , since it can be determined in a model-independent way from local measurements of the expansion rate.In practice, we should use the best-fit value obtained from supernova and Cepheid measurements [120], but since local measurements are currently in tension with the value inferred from Planck data, such a choice would mean that the fiducial cosmology is not ΛCDM.In particular, if H 0 differs from the Planck value, the theory of gravity should also be modified to ensure that the background evolution yields a distance to last scattering compatible with Planck.We therefore consider two separate fiducial models (see section 4.4): a ΛCDM fiducial model, based on the assumption that the H 0 tension is due to systematic effects and will disappear in the future, where we choose H 0 to be compatible with the Planck constraints; and a modified gravity fiducial model, where the values of the parameters are compatible with D M measured from Planck and with H 0 measured from supernovae and Cepheids.

Choice of fiducial models
We consider the three following cases: 1. ΛCDM fiducial with EFT parameters, Here we fix H 0 = 67.66km s −1 Mpc −1 according to the Planck measurement [32].
2. Modified gravity with large modifications (MGI).In this case, we assume that H 0 = 73.04km s −1 Mpc −1 is given by the local measurement of ref. [120].In order for the distance to last scattering to match the Planck value, other parameters affecting the background evolution must be modified, i.e. at least one among α fid M,0 , γ fid c,0 and w fid DE should be different from the ΛCDM value.At the same time, the stability condition in eq.(2.26) must be satisfied, which in general requires modifying at least two of these parameters.Once we adopt a modified gravity fiducial, the most general choice is that all parameters deviate from their ΛCDM values, unless some fine tuning is present.We thus choose fiducial values α fid M,0 , α fid B,0 and γ fid c,0 , which lie between the 1 and 2 σ constraints of the ΛCDM fiducial model.Then, the value of w fid DE is fixed such that the distance to last scattering matches the Planck measurement when assuming H 0 = 73.04km s −1 Mpc −1 .In summary, we choose for the EFT parameters Note that w DE < −1 does not lead to instabilities in Horndeski theories, since the effect of the equation of state on the sound speed can be counterbalanced by other parameters to keep c 2 s α > 0 (see eq. ( 5.1) below).
3. Modified gravity with small modifications (MGII).The values of the MGI fiducial model correspond to quite large deviations from ΛCDM.For comparison, we also consider a fiducial model with smaller deviations, to investigate the impact on the constraints with and without the gravitational redshift.We therefore include a third fiducial model defined by In this case, H 0 is fixed to the Planck value, α fid M,0 , α fid B,0 and γ fid c,0 are chosen arbitrarily, and w fid DE is determined from the distance to last scattering.
We always fix the values of A s and n s and the density parameters Ω m , Ω b , Ω r at z * according to the Planck constraints [32].We show the galaxy correlation multipoles in the three fiducial models in figure 1.The deviation of the even multipoles with respect to ΛCDM is due to a shift in the growth rate f , which is reduced in MGI and increased in MGII.The monopole is less affected by gravity modifications because it is dominated by the density contribution, which is less sensitive to deviations from GR than the velocity contribution.The dipole is additionally affected by the term due to gravitational redshift in eq.(3.10), which is zero in ΛCDM and positive in both MGI and MGII.
In addition to the equation of state and EFT parameters, the observables are affected by the fraction of CDM within a galaxy, x c , defined in eq.(3.12), which is generally unknown.The most straightforward possibility would be to assume that this fraction is equal to the fraction of CDM at the background level, as was done in ref. [39].However, this is not necessarily the case.Therefore, we include a free parameter X in our analysis, capturing the deviation of x c from the background CDM fraction ω c .In particular, we assume that the fractions of mass in CDM and baryons in each galaxy (denoted as x c , x b respectively) are related to the respective background energy fractions ω b , ω c as so that the condition x c + x b = 1 is always satisfied.The parameter X can have values X ∈ [0, 1 + ω b /ω c ], corresponding to the extreme cases x c = 1, x b = 0 for X = 1 + ω b /ω c , and x c = 0, x b = 1 for X = 0. 19 The fiducial value is assumed to be X = 1, in which case the fraction of CDM inside a galaxy corresponds to the background one.

MCMC analysis
In summary, our model is characterized by the following set of parameters: For each fiducial model in section 4.4, we generate mock values for the observables described above assuming zero noise.Then, we perform an MCMC analysis proceeding as follows.First, we evaluate the background quantities Ω m , Ω c and Ω r by solving the system of equations given in eqs.(2.12)-(2.14).We then obtain H(z) by integrating eq.(2.16) with initial condition H(z = 0) = H 0 , and compute the combination c 2 s α from eq. (2.26) with ζ given by eq.(2.6).Note that the derivative α ′ B can be expressed using the time evolution specified in eq.(4.2).If c 2 s α is negative, signaling the presence of an instability, the point is rejected.Otherwise, we insert the background solutions into the perturbations equations, solve the resulting system of equations in eqs.(2.22) and (2.23), and compute the expected values of all the observables.
The likelihood for the galaxy clustering observables is given by a multivariate Gaussian centered on the fiducial values.The respective covariance matrix includes cosmic variance due to the finite survey volume as well as shot noise arising from the finite number of observed galaxies (see appendix C of ref. [92] for detail).We include the non-vanishing cross-correlations between the even multipoles as well as between the bright and faint sample of galaxies.We also consider the pure cosmic variance contribution to the dipole from relativistic effects following ref.[111], where the robustness of the full theoretical expression for the covariance was thoroughly tested. 20For the angular diameter distance to the CMB, we assume a Gaussian likelihood centered on the Planck value D M = r * /θ * , where r * and θ * are taken from Planck 2018 [32].The variance is computed by propagating the error accordingly.We adopt large priors on all variables except for X, which is restricted to its physical range as already described in section 4.4.
We use the ensemble slice sampler zeus [121,122] with 18 walkers, ensuring that the lechangeh of the chains is at least 50 times the autocorrelation time, and further excluding a burn-in phase equal to the correlation lechangeh and thinning by the autocorrelation time.

Results
In this section, we present the constraints obtained on the EFT parameters α M,0 , α B,0 , γ c,0 and w DE .We discuss the results for the three fiducial models introduced in section 4.4 (ΛCDM and two modified gravity fiducials) and we additionally study the impact of restricting the parameter space to specific subclasses of models.We compare the analysis solely based on the even multipoles of the galaxy two-point correlation function (standard RSD analysis) with the constraints obtained when also including the dipole, which contains the contribution of gravitational redshift.

Constraints around a ΛCDM fiducial model
Choosing ΛCDM as our fiducial model, we show in figure 2 the constraints on the EFT parameters, marginalized over the biases and the parameter X.The results obtained with a standard RSD analysis are presented in blue, while those that also include gravitational redshift are in red.In both cases, we also consider the angular diameter distance in our set of observables.In Table 1, we report the marginal median values and credible intervals for the EFT parameters, both with and without the gravitational redshift ("Dipole" and "No dipole" columns).

Standard RSD analysis including even multipoles only
When considering the even multipoles only, the constraints arise from the galaxy velocity field and the angular diameter distance to last scattering, where the latter depends on the integral of the inverse Hubble parameter H −1 (z) from today to the redshift of decoupling.In the following, we discuss the respective constraints (blue contours in figure 2) for each parameter pair.
• w DE − α M,0 contour: These two parameters affect directly the evolution of H(z) and consequently the angular diameter distance to the CMB.From eq. ( 2.16) we see indeed that H(z) is directly sensitive to w DE , and that it also depends on α M through the 20 In particular, ref. [111] performed a comparison to the simulation-based jackknife covariance in the case of DESI, showing that the theoretical covariance that we also adopt here slightly underestimates the errors due to the impact of non-linearities.This only leads to a mild difference in the detectability of the dipole in the case of DESI.For SKA2, we expect a similar impact in the low-redshift sample, but a lower impact at high redshifts where non-linearities have less importance.Note that the complete expression for the covariance should also contain the impact of the survey mask.This can be evaluated once the survey is performed.Median values and 68% credible intervals for the ΛCDM fiducial model both for the standard RSD analysis (no dipole) and when adding the dipole that contains gravitational redshift.We also show the percent reduction on the credible interval obtained by including the dipole.evolution of the radiation density Ω r , see eq. (2.14).This leads to a strong degeneracy between α M,0 and w DE , as can be seen in figure 2. For example, an equation of state w DE < −1 leads to a dark energy pressure more negative than in ΛCDM, increasing the distance to last scattering.This can be compensated by a negative α M , which makes radiation decay more slowly than in ΛCDM: the increase in radiation pressure counteracts the negative pressure of dark energy, keeping the distance in agreement with CMB constraints.On the other hand, when w DE > −1, the opposite occurs, favoring a positive α M,0 .This leads to the two branches in the contours in figure 2. Note that the direction of the w DE −α M,0 degeneracy changes when w DE passes through −1.This is due to the fact that an equation of state significantly smaller than −1 can be easily compensated by increasing the pressure of radiation.On the other hand, an equation of state too close to zero would lead to no cosmic acceleration, which could not be compensated by a very low pressure from radiation.
• α M,0 − α B,0 and w DE − α B,0 contours: We obtain a first branch with w DE < −1 and α M < 0, in which case the stability of the perturbations requires α B < 0. Indeed, linearizing eq. ( 2.26) around ΛCDM and applying eq. ( 2.15), we find that (5.1) Since 1 + Ω r − 3Ω DE < 0 at low redshift, we see that α B < 0 is needed to keep c 2 s α > 0 when the first two terms are negative (note that α ′ B has the same sign as α B with our chosen time evolution).Moreover, since Ω DE and α M both decay at high redshift, while the term in front of α B increases at high redshift (it becomes less negative), a strongly negative α B,0 is necessary to maintain c 2 s α > 0 at higher redshift. 21This is clearly visible in the α M,0 − α B,0 and w DE − α B,0 contours of figure 2.
On the other hand, we have a second branch with w DE > −1 and α M > 0, which both contribute positively to the sound speed.Since the term proportional to α B in eq.(5.1) dominates over the first two at high redshift, α B,0 needs to be negative also in this branch to keep c 2 s α > 0 in the past.In this case, however, the absolute value of α B,0 does not need to be as large as in the other branch, as is visible from figure 2.
• Contours involving the parameter γ c,0 : In addition to the distance to last scattering and the stability condition, the parameters are constrained by the even multipoles, which are directly sensitive to the galaxy growth rate f .The growth rate depends on the evolution of the dark matter and baryon density given in eqs.(2.22) and (2.23), providing a constraint on the combination This explains the behavior of γ c,0 along the aforementioned two branches.Along the first branch, α M,0 and α B,0 are negative, and |α B,0 | > |α M,0 | (to ensure stability also in the past).To counterbalance this, and keep the growth rate close to ΛCDM, it is necessary to have γ c,0 > 0. Along the second branch, we have α M,0 > 0 and large, while α B,0 < 0, again requiring γ c,0 > 0 to keep the growth rate consistent with ΛCDM.This is clearly visible in the γ c,0 − α M,0 and γ c,0 − α B,0 contours of figure 2.
The degeneracy between γ c and α M can be understood from the Einstein and Euler equations: a positive γ c mainly enhances the growth of structure by changing the way dark matter falls in a gravitational potential, as can be seen from eq. (2.20).Indeed, the coupling of the scalar field to dark matter generates an additional force on dark matter particles.On the other hand, a positive running of the Planck mass α M in the first branch implies that the mass increases with time, leading to a decrease of the effective Newton constant in the Poisson equation, see eq. (2.21).This implies that a given matter density perturbation leads to a smaller distortion of space through the Poisson equation, thus decreasing the growth of structure and counterbalancing the enhancement arising from from γ c > 0. In the second branch, a negative α B counterbalances a positive γ c .The parameter α B , which encodes a mixing between the scalar and tensor kinetic terms, impacts the kinetic energy of the scalar field, which in turns affects the Poisson equation and therefore the growth of structure.
From figure 2, we notice that the posteriors are highly non-Gaussian, due to the strong degeneracies between the parameters.Such a behavior cannot be captured by a Fisher analysis, which assumes Gaussian posteriors.Therefore, our results show that it is essential to perform an MCMC analysis to fully understand the degeneracies in Horndeski theories and scenarios involving a breaking of the WEP, and correctly forecast the constraining power of the coming generation of surveys.

Including gravitational redshift by adding the dipole
The inclusion of the dipole in the analysis (red contours in figure 2) has a strong impact on the constraints, favoring one of the branches over the other.This is due to the fact that the dipole is directly sensitive to the combination as can be seen from eq. (3.10).From this, we notice that keeping the dipole close to its ΛCDM value requires either a small γ c or a cancellation of the two terms in the square bracket.This tightens the constraints on γ c , as is clearly visible in figure 2.Moreover, the degeneracy between γ c , α B and α M is broken by imposing that (α B − α M + 3γ c w c b c )/c 2 s α is tuned to balance the growth rate f .From figure 2, we see that adding this constraint by including the dipole almost completely removes one of the branches in the parameter space and significantly tightens the other.Since w DE is strongly degenerated with α M , the tightening of the constraints on α M from the dipole also improves the constraints on w DE .
In appendix B, figure 6, we show the constraints arising from the monopole alone and compare them to an analysis based on the combination of monopole and dipole.Also in this case, the dipole efficiently breaks the degeneracy between the parameters and removes one of the branches.The only difference with respect to the analysis with all multipoles is that the monopole is much less constraining on its own than in combination with the quadrupole and hexadecapole.This is due to the fact that the growth rate is degenerated with the galaxy bias using the monopole alone, whereas the inclusion of the other two even multipoles breaks this degeneracy and overall reduces the size of the two branches.
The improvement arising from the dipole is impressive: the constraints on the EFT parameters are tightened by up to 50%, as can be seen from Table 1.This is especially remarkable considering that the signal-to-noise ratio (SNR) of the dipole is significantly smaller than that of the even multipoles.Indeed, for a bias difference ∆b = 1 the cumulative SNR of the dipole in ΛCDM from separations d = 20 − 160 Mpc/h over the redshift range z = 0.15 − 1.55 is 80, while that of the quadrupole is 461, i.e. almost 6 times larger.
The reason why such a small signal can significantly improve the constraints is because deviations from ΛCDM affect it differently from the even multipoles.Hence, the key element to quantify its constraining power is the amplitude of the deviations with respect to ΛCDM.For example, if we choose a point at the edge of the 2 σ contours of figure 2 with α M,0 = 1.3, w DE = −0.8,α B,0 = −0.4 and γ c,0 = 0.17, we find that the SNR of the deviation is 6.2 for the quadrupole.For the dipole, the SNR of the WEP-breaking term (first line in eq.(3.10)) is 3.6, so not even two times smaller.Since this term has a very different dependence on the EFT parameters than the quadrupole, such an SNR is sufficient to significantly tighten the constraints.
The amplitude of the WEP-breaking term in the dipole in eq.(3.10) is directly sensitive to the galaxy bias difference between the two populations of galaxies, which we assume here to be equal to 1.This is what can be typically expected if the sample of galaxies is split into two populations according to their observed flux, for a population of galaxies similar to that of BOSS [89].The situation may be different for the galaxies detected by SKA2, for which the bias difference between two magnitude-limited samples may be smaller, decreasing the constraining power of the dipole.In such a case, it would be important to explore alternative ways of splitting galaxies to achieve the breaking of degeneracies we obtain in our forecasts, for example density splits where the galaxies are divided into two groups according to their environment.In this case, the bias difference between the two populations can become significantly larger [91,115].To mimic this, we have boosted the dipole by a factor 2, 5 and 10, see figure 7 in appendix B. We find that the contours are significantly tightened when the dipole is boosted, completely removing the second branch.This motivates the search for estimators of the dipole that would enhance the gravitational redshift contribution.
Finally, we notice that, while the dipole tightens the constraints on the model parameters, we do not see a significant impact on the nuisance parameters X, b 1,B , b 2,B , b 1,F and b 2,F .Concerning the four parameters specifying the biases of the two galaxy populations, this result matches our expectation that galaxy bias can be well measured from the even multipoles alone.The constraints on the relation X between the background CDM fraction and the average fraction of CDM inside a galaxy are always completely dominated by the physical prior specified below eq.(4.7), which luckily does not spoil the constraints on the model parameters.The corner plot for the full parameter space is included in appendix B, figure 8, showing the constraints from the even multipoles alone and from all multipoles.
We could further extend our analysis by varying the cosmological parameters Ω m , Ω b , Ω r at recombination as well as A s and n s in a joint analysis with CMB data.This would slightly broaden the constraints on our set of parameters.However, we do not expect this to affect the message of our work.A priori, letting the cosmological parameters vary should affect in a similar way the even multipoles and the dipole, such that the percentage improvement brought by the dipole should remain roughly the same.

Changing the parameterization for the dark energy evolution
The equation of state parameter of dark energy, w DE , is defined as the ratio of the dark energy pressure and the dark energy density.From eq. (2.11), we see that the evolution of the dark energy density depends on w DE and also on α M and γ c , where the latter governs α M, 0 This leads to the following expression for the Hubble parameter In figure 3, we show the constraints on the new set of parameters α M,0 , α B,0 , γ c,0 and ŵDE .As in the previous analysis, we see that the inclusion of the dipole significantly alleviates the strong degeneracies among the parameters that would otherwise be present.The direction of the degeneracies is however different with the new parameterization, which is not surprising given that the relation between w DE and ŵDE depends on γ c and α M .In particular, we notice that now an equation of state ŵDE < −1 is degenerated with a positive α M .This follows from the angular diameter distance constraint, which depends on the evolution of H(z).In this new parameterization, the latter depends directly on α M , γ c and ŵDE , as can be seen from eq. (5.6).From this equation we notice that ŵDE < −1 is counterbalanced by α M > 0 and ŵDE > −1 by α M < 0. As before, this leads to two branches in the parameter space when only the even multipoles are considered.We see that the constraints are clearly asymmetric around ΛCDM.This is particularly visible in the α M,0 − ŵDE contours, which are shifted towards positive values of α M,0 .This asymmetry is related to an interplay between the stability conditions and the constraints from the growth rate and the angular diameter distance.Indeed, linearizing the combination c 2 s α around ΛCDM with the new parameterization, we find The degeneracy in the α M,0 − ŵDE plane is governed by the distance to last scattering, which is sensitive to the combination α M + 3(1 + ŵDE )Ω DE through eq.(5.6).Along this degeneracy direction, when α M is positive, the positivity of the combination α M + (1 + ŵDE )Ω DE in eq.(5.7) is automatically maintained.This leads to a positive c 2 s α even when α B = 0, as long as γ c is small enough.Hence, the stability condition and the requirement of having a distance and growth rate close to ΛCDM are simultaneously satisfied.On the other hand, when α M is negative along the degeneracy direction, then α M + (1 + ŵDE )Ω DE < 0. To ensure stability, it is thus necessary to have negative values of α B and γ c .This leads to stronger deviations in the growth rate, which are disfavored.
As with the other parameterization, we see that adding the dipole strongly tightens the constraints on α M,0 by breaking the degeneracies with α B,0 and γ c,0 , which in turns also tightens the constraints on ŵDE .Quantitatively, from Table 2, we see that the 1 σ constraints on α M,0 and ŵDE are improved by 53% and 50% respectively when the dipole is added.On the other hand, the constraints on α B,0 are only mildly improved with this parameterization,  while those on γ c,0 show no improvement.Nevertheless, the red 2-dimensional contours involving γ c,0 in figure 3 show a clear preference of one branch when including gravitational redshift.Overall, comparing Table 1 and 2, we see that the dipole significantly tightens the constraints in both cases, but the way the improvement is distributed among the parameters clearly depends on the parameterization.

Constraints around modified gravity fiducials
As explained in section 4.4, we also explore constraints around two other fiducial models: one that would solve the Hubble tension (MGI) and another one that is closer to ΛCDM (MGII).The results for MGI are plotted in figure 4. We notice that the double branches around MGI are less pronounced in the blue contours (standard RSD analysis), than around ΛCDM (figure 2).This is due to the fact that the MGI fiducial model is well inside one of the branches of the ΛCDM contours.For example, we can see that the point with γ c,0 = 0.15 and α M,0 = −0.1 are well within the vertical branch of the γ c,0 −α M,0 contour in figure 2. As such, when we shift the fiducial model to this region, the degeneracies are already partially broken.While of course new degeneracies can appear around this new fiducial, figure 4 indicates that they are much less pronounced.As a consequence, including gravitational redshift still improves the constraints, but in a less significant way than around ΛCDM.From figure 4, we also notice that if MGI were to be the true model in our Universe, data from SKA2 would not be able to rule out ΛCDM based on the gravity + scalar sector only, as the 2 σ error bars on α M,0 , α B,0 and w DE encompass the corresponding ΛCDM values.However, we see that SKA2 would be able to favor MGI over ΛCDM based on the coupling to dark matter γ c,0 , which is incompatible with zero at more than 2 σ.Note that a model different from MGI might lead to a different conclusion regarding the relation with ΛCDM.
The constraints around MGII are shown in appendix B, figure 9.In this case, the degeneracies are still very pronounced with a standard RSD analysis.Since MGII is only mildly away from ΛCDM, the two branches around this point are quite similar to those around ΛCDM.As a consequence, we see that adding gravitational redshift significantly tightens the constraints, in a similar way as for the ΛCDM fiducial.

Constraints for specific sub-classes of Horndeski theories
As discussed above, the strechangeh of gravitational redshift is to break the degeneracy between the parameters γ c,0 , α B,0 and α M,0 .As a comparison, we investigate specific models where only two of the three parameters are left free.In figure 5 (left panel), we show the contours for the case of Horndeski theories without a breaking of the WEP, i.e. with γ c = 0.In the right panel we show another specific case, where the parameter γ c,0 is left free, but we fix α B = α M /2.This corresponds to various Horndeski theories (cf.Table 1 in ref. [8]), including Brans-Dicke theories [17,123] and some f (R) models [14,[124][125][126].From these figures, we see that the constraints from the RSD analysis are much tighter than before.The degeneracies in the parameter space are indeed much less severe when one parameter is removed.In particular, the strong degeneracy between α M,0 and w DE is reduced, due to the increased constraints on α M,0 from the growth of structure f .In these cases, gravitational redshift does not improve the constraints further.This is consistent with the SNR comparison made above: the statistical power of the dipole term (the full expression in eq.(3.10)) is ∼ 6 times smaller than that of the quadrupole.As such, when there is no strong degeneracy to break, the dipole has almost no impact on the constraints.This is also in agreement with the results of ref. [127], who showed that relativistic effects do not help to constrain Horndeski theories with no breaking of the WEP.
A word of caution is however necessary concerning these results: both the analysis of ref. [127] and our analysis assume a fixed time evolution for the parameters.Even if we have chosen this time evolution to be the same for all parameters in order to minimize its impact on the degeneracies, the simple fact that we assume a known time evolution helps breaking some of the degeneracies.For example, in eq.(5.2), the term w c b c has a different time evolution from γ c and α M and therefore helps breaking the degeneracy between these parameters in the model with α B = α M /2, when we combine multiple redshift bins.If we consider one redshift bin only, i.e. we do not exploit the fact that we have a fixed time evolution for the parameters, we find that the parameters are completely degenerated and cannot be constrained anymore.

Conclusions
In this paper, we show that the effect of gravitational redshift, which will be measurable with the coming generation of galaxy surveys, is crucial to efficiently test models beyond GR.To determine the constraining power of this new observable in a generic way, we have considered the effective theory of interacting dark energy, which encompasses all scalar-tensor theories of gravity and allows for a non-minimal coupling of dark matter.We find that adding gravitational redshift to a standard RSD analysis significantly improves the constraints on the EFT parameters, by up to 50 percent.The precise results depend on the choice of parameterization within the EFT framework, on the time evolution assumed for the EFT functions and on the fiducial cosmological model, but the overall message remains unaffected by these details: gravitational redshift breaks degeneracies between parameters, leading to a significant improvement in the constraints.
To measure this novel observable, no new data are involved: one can simply fit for a dipole in the cross-correlation function between two different populations of galaxies.This requires minimal changes to the analysis pipeline, but leads to a strong impact on the constraints.The improvement is especially remarkable given that the SNR of the dipole is typically six times smaller than that of the standard quadrupole term.However, the specificity of the dipole is that it contains a particular combination of gravitational redshift and Doppler effects, which exactly vanishes if the WEP is valid.This combination is not present in the standard even multipoles of the correlation function and is able to break the degeneracies between the EFT parameters.Note that, since this combination is proportional to the bias difference between the two populations of galaxies, we expect the breaking of degeneracies to scale with the value of this difference.Here, we have assumed a bias difference equal to 1, but estimators that would boost this difference may improve the constraints further, while a lower difference may result in an overall worsening.
In this work we have focused on spectroscopic surveys, which are directly sensitive to the growth rate of structure through RSD.In a future work, we will study whether the inclusion of gravitational lensing can further improve the constraints.Since gravitational lensing is sensitive to the sum of the two gravitational potentials (the spatial and temporal perturbations in the metric), while gravitational redshift only depends on the distortion of time, combining the two may lead to more stringent constraints on the effective theory of interacting dark energy.Possible future work could also include an investigation beyond the quasi-static ap-proximation assumed in our work and the choice of other time parameterizations for the EFT parameters.
We conclude that measuring gravitational redshift provides an efficient way of increasing the constraining power of large-scale structure surveys at minimal cost, and a combination with additional observables could further increase its remarkable constraining power.corner plot in figure 6 shows the constraints obtained on the parameters α M,0 , α B,0 , γ c,0 and w DE when including only the monopole of the galaxy correlation function (yellow contours), compared to the combination with the dipole term, which contains gravitational redshift (green contours).In figure 7, we further illustrate the impact of a dipole boosted by factors 2, 5, and 10 (beige-brown shaded contours), compared to the base case (red contours) without a boost.In figure 8 we show the full corner plot, including all parameters specified in eq.(4.8), for the standard RSD analysis (blue contours) and also including gravitational redshift (red contours).Finally, in figure 9 we show the constraints around the modified gravity fiducial model MGII.The contours are very similar to those obtained around ΛCDM, see figure 2 α fid M,0 = −0.1 , α fid B,0 = −0.5 , γ fid c,0 = 0.15 , w fid DE = −1.41 .(MGI fiducial) (4.5)

15 Figure 1 .
Figure 1.Multipoles of the galaxy correlation function as a function of separation, d, for ΛCDM (green) and the modified gravity fiducial models MGI (red) and MGII (blue), evaluated at the center of the first SKA redshift bin, z = 0.15.

Figure 4 .
Figure 4. Constraints around the modified gravity fiducial MGI.

Figure 5 .
Figure 5. Constraints given a ΛCDM fiducial when one parameter is kept fixed.Left panel : γ c is fixed to zero.Right panel : α B = α M /2. .

Figure 6 .
Figure 6.Constraints around the ΛCDM fiducial when including only the monopole, and combining monopole and dipole.

Figure 7 .
Figure 7. Constraints around the ΛCDM fiducial when including all multipoles and applying various boost factors to the dipole.

Figure 8 .
Figure 8. Constraints around a ΛCDM fiducial for the full parameter space.

Figure 9 .
Figure 9. Constraints around the modified gravity fiducial MGII.The constraints are very similar to those around ΛCDM, see figure 2

Table 2 .
Median values and 68% credible intervals for the ΛCDM fiducial model both for the standard RSD analysis (no dipole) and when adding the dipole that contains gravitational redshift, employing the alternative parameterization with the effective equation of state parameter ŵDE .For the credible interval, we also show the percent reduction obtained by the latter.theinteraction of dark energy with dark matter.Alternatively, we can define an effective equation of state, ŵDE , directly encoding the evolution of dark energy through ρ′ DE = −3(1 + ŵDE )ρ DE .