# Gravitational Redshift Constraints on the Effective Theory of Interacting Dark Energy

Sveva Castello, Michele Mancarella, Nastassia Grimm, Daniel Sobral Blanco, Isaac Tutusaus, Camille Bonvin

► To cite this version:

Sveva Castello, Michele Mancarella, Nastassia Grimm, Daniel Sobral Blanco, Isaac Tutusaus, et al. Gravitational Redshift Constraints on the Effective Theory of Interacting Dark Energy. Journal of Cosmology and Astroparticle Physics, 2024, 05, pp.003. 10.1088/1475-7516/2024/05/003. hal-04344189

HAL Id: hal-04344189

https://hal.science/hal-04344189v1

Submitted on 17 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
PAPER • OPEN ACCESS

# Gravitational redshift constraints on the effective theory of interacting dark energy

To cite this article: Sveva Castello et al JCAP05(2024)003

View the [article online](#) for updates and enhancements.

## You may also like

- THE EFFECTS OF STELLAR ROTATION, IL A COMPREHENSIVE SET OF STARBURST99 MODELS

Claus Leitherer, Sylvia Ekström, Georges Meynet et al.

- TIGRE: a MATLAB-GPU toolbox for CBCT image reconstruction

Ander Biguri, Manjit Dosanjh, Steven Hancock et al.

- The SST-1M imaging atmospheric Cherenkov telescope for gamma-ray astrophysics

C. Alispach, A. Araudo, M. Balbo et al.
RECEIVED: December 7, 2023

REVISED: February 12, 2024

ACCEPTED: March 18, 2024

PUBLISHED: May 2, 2024

# Gravitational redshift constraints on the effective theory of interacting dark energy

Sveva Castello<sup>D,a,*</sup> Michele Mancarella<sup>D,b,c,*</sup> Nastassia Grimm<sup>D,a</sup>

Daniel Sobral-Blanco<sup>D,a</sup> Isaac Tutusaus<sup>D,d,a</sup> and Camille Bonvin<sup>D,a</sup>

**ABSTRACT:** 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).

**KEYWORDS:** dark energy theory, galaxy clustering, galaxy surveys, redshift surveys

**ARXIV EPRINT:** 2311.14425

<figure><img src="image_2.png" /></figure>

<footnote><sup>a</sup>Département de Physique Théorique and Center for Astroparticle Physics, Université de Genève, Quai E. Ansermet 24, CH-1211 Genève 4, Switzerland</footnote>

<footnote><sup>b</sup>Dipartimento di Fisica “G. Occhialini”, Università degli Studi di Milano-Bicocca, Piazza della Scienza 3, 20126 Milano, Italy</footnote>

<footnote><sup>c</sup>INFN, Sezione di Milano-Bicocca, Piazza della Scienza 3, 20126 Milano, Italy</footnote>

<footnote><sup>d</sup>Institut de Recherche en Astrophysique et Planétologie (IRAP), Université de Toulouse, CNRS, UPS, CNES, 14 Av. Edouard Belin, 31400 Toulouse, France

E-mail: sveva.castello@unige.ch, michele.mancarella@unimib.it, Nastassia.Grimm@unige.ch, Daniel.SobralBlanco@unige.ch, isaac.tutusaus@irap.omp.eu, Camille.Bonvin@unige.ch</footnote>

<footnote>*Corresponding author.</footnote>
# Contents

<table><tr><td>1 Introduction</td><td>1</td></tr><tr><td>2 Effective theory of interacting dark energy</td><td>4</td></tr><tr><td>2.1 Gravitational and matter sector</td><td>4</td></tr><tr><td>2.2 Background</td><td>5</td></tr><tr><td>2.3 Perturbations</td><td>6</td></tr><tr><td>3 Galaxy clustering observables</td><td>8</td></tr><tr><td>4 Survey specifications and numerical implementation</td><td>11</td></tr><tr><td>4.1 Survey specifications</td><td>11</td></tr><tr><td>4.2 Choices of parameterizations and time evolution of the free functions</td><td>12</td></tr><tr><td>4.3 Observables</td><td>13</td></tr><tr><td>4.4 Choice of fiducial models</td><td>14</td></tr><tr><td>4.5 MCMC analysis</td><td>16</td></tr><tr><td>5 Results</td><td>17</td></tr><tr><td>5.1 Constraints around a ΛCDM fiducial model</td><td>17</td></tr><tr><td>5.2 Constraints around modified gravity fiducials</td><td>24</td></tr><tr><td>5.3 Constraints for specific sub-classes of Horndeski theories</td><td>25</td></tr><tr><td>6 Conclusions</td><td>26</td></tr><tr><td>A Relation to the phenomenological parameters</td><td>27</td></tr><tr><td>B Supplementary plots</td><td>28</td></tr></table>

## 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–9], which encompasses a large class of theories in a unified formalism, describing a rich phe-
nomenology beyond GR.<fnref n="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 scalar-tensor 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–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 $\alpha_K$, the braiding $\alpha_B$, the Planck-mass run rate $\alpha_M$ and the tensor speed excess $\alpha_T$.<fnref n="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–37]. However, at the moment only the parameter $\alpha_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–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–41].<fnref n="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. refs. [47–51]), while those on cold dark matter (CDM) could be significant and indeed occur in various theories [52–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 $\gamma_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

<footnote n="1">¹An alternative model-independent framework, which was developed to test gravity with cosmological observations on linear scales, is the Parameterized Post-Friedmann approach, see refs. [10–12].</footnote>

<footnote n="2">²An additional function of time, $\alpha_H$, is required to describe “beyond Horndeski” theories [22, 23], while four further functions are necessary to describe degenerate higher-order theories (DHOST, see refs. [24–26]).</footnote>

<footnote n="3">³Note that other bounds exist for functions describing theories beyond Horndeski, which further motivates our choice of restricting to Horndeski theories here. In particular, imposing the absence of gravitational-wave-induced instabilities in the dark energy sector places a tight bound of $|\alpha_H| \lesssim 10^{-20}$ [36]. Free functions of DHOST models can be constrained by astrophysical observations [42–46].</footnote>
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–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–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–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 RSD makes 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–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 2-dimensional 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\sigma$ error bars by up to $\sim 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.

## 2 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.

### 2.1 Gravitational and matter sector

As in refs. [39, 62], we assume that the gravitational sector is described by a four-dimensional metric $g_{\mu\nu}$ and a scalar field $\phi$, 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<fnref n="4" />

$$ds^2 = -N^2 dt^2 + h_{ij} (dx^i + N^i dt) (dx^j + N^j dt), \quad (2.1)$$

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

$$S_g = \int d^4 x \sqrt{-g} L(N, K_{ij}, R_{ij}, h_{ij}, D_i; t), \quad (2.2)$$

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

$$S_m = \sum_{I}^{N_S} S_I, \quad S_I = \int d^4 x \sqrt{-\check{g}^{(I)}} L_I(\check{g}_{\mu\nu}^{(I)}, \psi_I), \quad (2.3)$$

<footnote n="4">⁴We work in units where the speed of light $c$ is set to 1.</footnote>
with each of the $N_S$ particle species being minimally coupled to the gravitational sector by a generally different effective metric $\tilde{g}_{\mu\nu}^{(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_{\mu\nu}$. In this frame, CDM is allowed to be non-minimally coupled to $g_{\mu\nu}$. We denote baryons, CDM and radiation by the subscripts b, c and r respectively.

## 2.2 Background

We consider a homogeneous and isotropic background described by the Friedmann-Lemaître-Robertson-Walker (FLRW) metric

$$ds^2 = -dt^2 + a^2(t)\delta_{ij}dx^i dx^j, \qquad (2.4)$$

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 $\bar{S}_g + \bar{S}_m$.

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])

$$H^2 = \frac{1}{3M^2} (\bar{\rho}_b + \bar{\rho}_c + \bar{\rho}_r + \bar{\rho}_{\text{DE}}), \qquad (2.5)$$

$$\zeta \equiv \frac{H'}{H} = -\frac{3}{2} - \frac{1}{2M^2H^2} (\bar{\rho}_{\text{DE}} + \bar{\rho}_r). \qquad (2.6)$$

Here, $\rho$ 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 $\alpha_M$, defined by

$$\alpha_M \equiv \frac{d \ln M^2}{d \ln(a)}. \qquad (2.7)$$

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

$$\bar{\rho}'_b = -3 \bar{\rho}_b, \qquad (2.8)$$

$$\bar{\rho}'_c = -3(1 - \gamma_c)\bar{\rho}_c, \qquad (2.9)$$

$$\bar{\rho}'_r = -4\bar{\rho}_r. \qquad (2.10)$$

Hence, baryons and photons evolve in a standard way, as perfect fluids with pressure $p_b = 0$ and $p_r = \rho_r/3$ respectively, while the evolution of CDM is modified by its non-minimal
coupling to gravity, encoded in the free function $\gamma_c$.<fnref n="5" /> Eqs. (2.5)–(2.6) define the energy density and pressure of the dark energy component. The dark energy equation of state is then given by $w_{DE} \equiv \bar{p}_{DE}/\bar{\rho}_{DE}$, and is included as a free parameter in the model. Inserting this into eqs. (2.5)–(2.6), we see that the dark energy density evolves as

$$ \bar{\rho}'_{DE} = -3(1 + w_{DE})\bar{\rho}_{DE} - 3\gamma_c\bar{\rho}_c + \alpha_M(\bar{\rho}_b + \bar{\rho}_c + \bar{\rho}_r + \bar{\rho}_{DE}). \quad (2.11) $$

It is convenient to introduce the dimensionless density parameters $\Omega_X \equiv \frac{\bar{\rho}_X}{3M^2H^2}$, and rewrite the evolution equations for $\bar{\rho}_b$, $\bar{\rho}_c$ and $\bar{\rho}_r$ as

$$ \Omega'_b = -\Omega_b [\alpha_M - 3w_{DE}\Omega_{DE} - \Omega_r], \quad (2.12) $$

$$ \Omega'_c = -\Omega_c [\alpha_M - 3\gamma_c - 3w_{DE}\Omega_{DE} - \Omega_r], \quad (2.13) $$

$$ \Omega'_r = -\Omega_r [1 + \alpha_M - 3w_{DE}\Omega_{DE} - \Omega_r]. \quad (2.14) $$

In the above equations, $\Omega_{DE}$ is defined via eq. (2.5) by $\Omega_{DE} = 1 - \Omega_b - \Omega_c - \Omega_r$. The Hubble parameter can be calculated from the solutions of eq. (2.12) to (2.14) by writing eq. (2.6) as

$$ \zeta = -\frac{3}{2} \left( 1 + \frac{1}{3} \Omega_r + w_{\mathrm{E}} \Omega_{\mathrm{DE}} \right), \quad (2.15) $$

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$),

$$ H(z) = H_0 \exp \left[ \frac{3}{2} \int_{\ln(a(z))}^{0} d\ln a' \left( 1 + \frac{1}{3}\Omega_r + w_{\mathrm{DE}}\Omega_{\mathrm{DE}} \right) \right]. \quad (2.16) $$

## 2.3 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 $\alpha_B$, $\alpha_K$ and $\alpha_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

$$ S_g^{(2)} = \int d^3x \, dt \, a^3 \frac{M^2}{2} \left[ \delta K_j^i \delta K_i^j - \delta K^2 + R \delta N + (1 + \alpha_T) \delta_2 \left( \sqrt{h} R / a^3 \right) \right. \\ \left. + \alpha_K H^2 \delta N^2 + 4 \alpha_B H \delta K \delta N \right]. \quad (2.17) $$

Here, $\delta_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 $\alpha_T$, the kineticity $\alpha_K$, the braiding $\alpha_B$,<fnref n="6" /> and the effective Planck mass $M$, whose evolution is encoded in the parameter $\alpha_M$ defined in eq. (2.7). The parameter $\alpha_T$

<footnote n="5">⁵The function $\gamma_c$ is given by a combination of the two functions governing the conformal and disformal coupling of CDM to gravity, see eqs. (2.24) and (2.16)–(2.17) of ref. [39]. As a consequence, in this framework it is not possible to distinguish between conformal and disformal couplings. These would be distinguishable only for a component with pressure, see section 5 of ref. [62].</footnote>

<footnote n="6">⁶Note that the parameter $\alpha_B$ used here corresponds to the one used in refs. [39, 62], but needs to be divided by a factor -2 to obtain the parameter $\alpha_B$ originally introduced in ref. [8].</footnote>
is directly related to the speed of propagation of tensor modes $c_T \equiv 1 + \alpha_T$, which was constrained to match the speed of light up to $\sim 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 $\alpha_T = 0$ throughout our analysis.<fnref n="7" /> We are thus left with the three free functions $\alpha_M, \alpha_K$ and $\alpha_B$, which take the values $\alpha_K = \alpha_B = \alpha_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,

$$ds^2 = -(1+2\Psi)dt^2 + a^2(t)(1-2\Phi)\delta_{ij}dx^i dx^j. \qquad (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 $\Phi$ and $\Psi$ are interchanged),<fnref n="8" />

$$V'_{b} + V_{b} - \frac{k}{\mathcal{H}}\Psi = 0, \qquad (2.19)$$

$$V'_{c} + V_{c}(1 + 3\gamma_{c}) - \frac{k}{\mathcal{H}}\Psi \left(1 + \frac{6\gamma_{c}}{c_{s}^{2}\alpha\mu}(\alpha_{B} - \alpha_{M} + 3\gamma_{c}\omega_{c}b_{c})\right) = 0. \qquad (2.20)$$

Here, $V_{b,c}$ are the velocity potentials of baryons and CDM in Fourier space, related to the velocity fields by $V_{b,c} = -i\mathbf{k}/k V_{b,c}$, and the quantities $c_s^2\alpha$ and $\mu$ are defined below. The Poisson equation reads (see eq. (4.19) in ref. [62])

$$-\frac{k^2}{a^2}\Psi = 4\pi G\mu\delta_m\bar{\rho}_m, \qquad (2.21)$$

i.e. the Newton constant $G$ is replaced by an effective Newton constant, $G\mu$, which depends on the gravity modifications and the non-minimal coupling $\gamma_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]),

$$\delta''_{b} = -(2 + \zeta) \delta_{b} + \frac{3}{2} \Omega_{m} \left[ 1 + \frac{2}{c_{s}^{2}\alpha} (\alpha_{B} - \alpha_{M}) (\alpha_{B} - \alpha_{M} + 3\gamma_{c}\omega_{c}b_{c}) \right], \qquad (2.22)$$

$$\delta''_{c} = -(2 + \zeta + 3\gamma_{c}) \delta_{c} + \frac{3}{2} \Omega_{m} \left[ 1 + \frac{2}{c_{s}^{2}\alpha} (\alpha_{B} - \alpha_{M} + 3\gamma_{c}) (\alpha_{B} - \alpha_{M} + 3\gamma_{c}\omega_{c}b_{c}) \right]. \qquad (2.23)$$

Here, $\Omega_m$ and $\delta_m$ refer to the total matter contribution,

$$\Omega_m = \Omega_b + \Omega_c, \qquad \delta_m = (\Omega_b\delta_b + \Omega_c\delta_c)/\Omega_m, \qquad (2.24)$$

<footnote n="7">⁷An open possibility is that the propagation speed of GWs exhibits a frequency dependence such that the GR value is recovered in the frequency range where the LIGO-Virgo detectors operate, while deviating from GR at lower frequencies [96]. This effect can be constrained with LISA or multi-band GW observations [97, 98]. In this work, we assume a frequency-independent GW speed.</footnote>

<footnote n="8">⁸The velocity potentials $v_{b,c}$ in ref. [62] are related to the quantities $V_{b,c}$ used here as $v_{b,c} = -aV_{b,c}/k$.</footnote>
and $\omega_c \equiv \Omega_c/\Omega_m$ and $b_c \equiv \delta_c/\delta_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 $\zeta$ given in eq. (2.6), and on the function $\mu$ defined as

$$ \mu \equiv 1 + \frac{2}{c_s^2 \alpha} (\alpha_B - \alpha_M) (\alpha_B - \alpha_M + 3\gamma_c\omega_c b_c). \qquad (2.25) $$

Moreover, the evolution equations contain the quantity $c_s^2\alpha$, involving the speed of sound of scalar perturbations $c_s$ (given by eq. (2.39) of ref. [62]),

$$ c_s^2\alpha = -2 \left[ (1 + \alpha_B)(\zeta - \alpha_M + \alpha_B) + \alpha'_B + \frac{3}{2}\Omega_m + 2\Omega_r \right]. \qquad (2.26) $$

The quantity $\alpha$, defined in eq. (2.36) of ref. [62], contains a combination of the functions $\alpha_K$, $\alpha_B$ and $\alpha_{D,c}$, where the latter arises from a disformal coupling of CDM to the metric $g_{\mu\nu}$. However, we note that the combination $c_s^2\alpha$ in eq. (2.26) does not depend on $\alpha_K$ and $\alpha_{D,c}$. This is a consequence of the quasi-static approximation and implies that the kineticity $\alpha_K$ remains unconstrained when adopting this limit.<fnref n="9" /> Finally, the absence of ghost and gradient instabilities requires $c_s^2\alpha \ge 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 $\alpha_M$, $\alpha_B$, $\gamma_c$ and $w_{DE}$, highlighting the fundamental role played by gravitational redshift.

## 3 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 $\hat{n}$ in the sky. We can use these observations to construct the number count fluctuations observable,

$$ \Delta \equiv \frac{N(z, \hat{n}) - \bar{N}(z)}{\bar{N}(z)}, \qquad (3.1) $$

where $\bar{N}(z)$ is the average number of galaxies at redshift $z$. The quantity $\Delta$ can be computed at linear order in perturbation theory [100–104], yielding<fnref n="10" />

$$ \Delta(z, \hat{n}) = \Delta^{\text{st}}(z, \hat{n}) + \Delta^{\text{rel}}(z, \hat{n}), \qquad (3.2) $$

where

$$ \Delta^{\text{st}}(z, \hat{n}) = b(z) \delta_m(z, \hat{n}) - \frac{1}{\mathcal{H}} \partial_r (\mathbf{V} \cdot \hat{n}), \qquad (3.3) $$

$$ \Delta^{\text{rel}}(z, \hat{n}) = \frac{1}{\mathcal{H}} \partial_r \Psi + \frac{1}{\mathcal{H}} \dot{\mathbf{V}} \cdot \hat{n} + \left[ 1 - \frac{\dot{\mathcal{H}}}{\mathcal{H}^2} - \frac{2}{r\mathcal{H}} - 5s(z) \left( 1 - \frac{1}{r\mathcal{H}} \right) \right] \mathbf{V} \cdot \hat{n}. \qquad (3.4) $$

<footnote n="9">⁹We note that, even when the full set of perturbations equations beyond the quasi-static limit is considered, the kineticity is hardly constrained by CMB, baryon acoustic oscillations, and RSD data, and in general it can be fixed to arbitrary values without affecting the results for the remaining parameters [35].</footnote>

<footnote n="10">¹⁰For a derivation of the relativistic effects at second order, see e.g. [105–107].</footnote>
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_{\mu\nu}$. Eq. (3.3) contains the two dominant contributions to $\Delta$, 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, $\Delta$ 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.<fnref n="11" /> Here, 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].<fnref n="12" />

The standard approach to analyze galaxy clustering data is to apply summary statistics. In this work, we consider the galaxy two-point correlation function,<fnref n="13" /> $\xi \equiv \langle \Delta(\hat{n}, z)\Delta(\hat{n}', z') \rangle$, which can be expanded in multipoles,

$$ \xi_{\ell}(z, d) = \frac{2\ell + 1}{2} \int_{-1}^{1} d\mu \, \xi(z, d, \mu) P_{\ell}(\mu). \qquad (3.5) $$

Here, $\mu$ 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_{\ell}(\mu)$ denotes the Legendre polynomial of order $\ell$.

In the distant-observer regime, the standard contributions (density fluctuations and RSD) generate three even multipoles: a monopole ($\ell = 0$), a quadrupole ($\ell = 2$) and a hexadecapole ($\ell = 4$), which were measured by various surveys, see e.g. ref. [88]. These are given by

$$ \xi_0(z, d) = \left[ b^2(z) + \frac{2}{3}b(z)f(z) + \frac{1}{5}f^2(z) \right] \mu_0(z, d), \qquad (3.6) $$

$$ \xi_2(z, d) = -\left[\frac{4}{3}f(z)b(z) + \frac{4}{7}f^2(z)\right] \mu_2(z, d), \qquad (3.7) $$

$$ \xi_4(z, d) = \frac{8}{35}f^2(z)\mu_4(z, d), \qquad (3.8) $$

<footnote n="11">¹¹Other contributions appearing in the full expression given in refs. [100–102], including gravitational lensing, are subdominant in the redshift range relevant for this work [109, 110].</footnote>

<footnote n="12">¹²The evolution bias $f^{\text{evol}}$ depends on the redshift evolution of the populations of sources. It can be directly measured from the average number of galaxies once data become available (see ref. [112] for a forecast using SKA specifications). However, as it does not affect the terms proportional to $x_c$ encoding the breaking of the WEP, we neglect it here, since we do not expect it to impact the constraining power of the dipole. Note that in general the contribution of $f^{\text{evol}}$ to the dipole is suppressed with respect to the other contributions for two reasons. First, from eq. (8) in ref. [111], we see that the contributions proportional to the magnification bias difference is boosted at low redshift over those proportional to the evolution bias by a factor $(1 - \frac{1}{r^H})$. Secondly, following the argument in section 4.1 of ref. [111], we see that the term depending on the difference $f_B^{\text{evol}} - f_F^{\text{evol}}$ vanishes, further reducing the impact of the evolution bias.</footnote>

<footnote n="13">¹³We choose to work with the two-point correlation function rather than its Fourier-space counterpart, the power spectrum, since this makes it easier to implement the wide-angle correction in eq. (3.10).</footnote>
where $f = d \ln(\delta_m)/d \ln(a)$ is the matter growth rate, and we have defined

$$ \mu_{\ell}(z, d) = \int \frac{dk k^2}{2\pi^2} P_{\delta\delta}(k, z) j_{\ell}(kd). \qquad (3.9) $$

Here, $P_{\delta\delta}(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 $\Delta^{\text{rel}}$ is suppressed by a factor of $(\mathcal{H}/k)^2$ and is therefore negligible on sub-horizon scales. On the other hand, the cross-correlation of $\Delta^{\text{rel}}$ with the standard terms breaks the symmetry of the two-point function $\xi$ [80, 81, 84], generating a dipole moment.<fnref n="14" /> To 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,

$$ \xi_1^{\text{BF}}(z, d) = \frac{\mathcal{H}}{\mathcal{H}_0} \nu_1(d, z) \left\{ \begin{aligned} & 3 x_c (b_B - b_F) \gamma_c \left[ f - (1+z) \left( \frac{\mathcal{H}_0}{\mathcal{H}} \right)^2 \frac{3\Omega_{m,0}}{c_s^2\alpha} (\alpha_B - \alpha_M + 3\gamma_c\omega_c b_c) \right] \\ & - 3(s_B - s_F) f^2 \left(1 - \frac{1}{r\mathcal{H}}\right) + (b_B - b_F) f \left(\frac{2}{r\mathcal{H}} + \frac{\dot{\mathcal{H}}}{\mathcal{H}^2}\right) \\ & + 5f (b_B s_F - b_F s_B) \left(1 - \frac{1}{r\mathcal{H}}\right) \end{aligned} \right\} \\ & - \frac{2}{5} (b_B - b_F) f \frac{d}{r} \mu_2(d, z). \qquad (3.10) $$

Here, $x_c$ is the fraction of CDM mass in a galaxy (which in general can differ from the background CDM fraction $\omega_c$ as well as from the fraction $b_c \equiv \delta_c/\delta_m$ at the level of perturbations), and

$$ \nu_1(z, d) = \int \frac{dk k}{2\pi^2} \mathcal{H}_0 P_{\delta\delta}(k, z) j_{\ell}(kd). \qquad (3.11) $$

The first line in eq. (3.10) arises from the first three terms of $\Delta^{\text{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 $\gamma_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 $\mu_2(d)$, which is enhanced by a factor $k/\mathcal{H}$ with respect to $\nu_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].<fnref n="15" />

<footnote n="14">¹⁴The cross-correlation between standard and relativistic terms also generates an octupole, which is however not sensitive to gravitational redshift [66] and thus does not help in breaking degeneracies in our analysis. We will therefore not consider it in the following.</footnote>

<footnote n="15">¹⁵When considering the galaxy power spectrum instead of the correlation function, the wide-angle term generates a dipole for a single population of galaxies because of the effect of the window function [114].</footnote>
The galaxy velocity $V$ appearing in eqs. (3.3) and (3.4) is a weighted average of the velocities of baryons and CDM [39],

$$V = x_c V_c + (1 - x_c) V_b. \qquad (3.12)$$

The velocity potentials $V_{b,c}$ can be related to the density fields $\delta_{b,c}$ via

$$V_{b,c} = -\frac{\mathcal{H}}{k} f_{b,c}(a) \delta_{b,c}(a, \mathbf{k}), \qquad (3.13)$$

where $f_{b,c}(a) = \frac{d\ln(\delta_{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

$$V = -\frac{\mathcal{H}}{k} [x_c f_c \delta_c + (1 - x_c) f_b \delta_b] = -\frac{\mathcal{H}}{k} f \delta_m, \qquad (3.14)$$

where the growth rate of galaxies $f$ is related to that of baryons and CDM through

$$f = x_c f_c b_c + (1 - x_c) f_b b_b. \qquad (3.15)$$

In the following sections, we forecast the constraints on the effective theory of interacting dark energy arising from measurements of the even multipoles of $\xi$, and we compare them with the constraints obtained when the dipole is added to the analysis.

# 4 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.

## 4.1 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 $\Delta z = 0.1$ and in each of these bins we consider the survey specifications given in table 3 in ref. [94].$^{16}$ We 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],

$$b_P(z) = b_{1,P} \exp(z \cdot b_{2,P}) \pm \Delta b, \qquad (4.1)$$

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$

<footnote>$^{16}$The highest five redshift bins between $z = 1.5$ and $z = 2.0$ stated in table 3 of ref. [94] could as well be included in our analysis, but we omit them for numerical efficiency as they do not provide a significant improvement to our constraints.</footnote>
following ref. [94]. We fix $\Delta 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].

## 4.2 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 $\alpha_B$, the running of the Planck mass $\alpha_M$, the parameter $\gamma_c$ encoding a non-minimal coupling of CDM to the metric $g_{\mu\nu}$, 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 $\alpha_B$, $\alpha_M$ and $\gamma_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 \in \{\alpha_B, \alpha_M, \gamma_c\}$, we write

$$p = p_0 \frac{1 - \Omega_m^{\Lambda\text{CDM}} - \Omega_r^{\Lambda\text{CDM}}}{1 - \Omega_{m,0}^{\Lambda\text{CDM}} - \Omega_{r,0}^{\Lambda\text{CDM}}}, \quad (4.2)$$

where the superscript $\Lambda$CDM denotes quantities evaluated in $\Lambda$CDM, i.e. explicitly $\Omega_m^{\Lambda\text{CDM}} = \Omega_{m,0}^{\Lambda\text{CDM}}(H_0/H^{\Lambda\text{CDM}}(z))^2 \exp(-3 \ln a)$ and $\Omega_r^{\Lambda\text{CDM}} = \Omega_{r,0}^{\Lambda\text{CDM}}(H_0/H^{\Lambda\text{CDM}}(z))^2 \exp(-4 \ln a)$. Note that here we parameterize the time evolution using the fiducial $\Lambda$CDM evolution for simplicity. Given that the background evolution is expected to be close to $\Lambda$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 $\alpha_{B,0}$, $\alpha_{M,0}$ and $\gamma_{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 $\alpha_M$ and $\alpha_B$, but it differs concerning $\gamma_c$. In ref. [39], the authors assumed a constant value for the combination $\beta_\gamma = 3\sqrt{2}\gamma_c/(c_s\alpha^{1/2})$ and provided constraints on this quantity.<fnref n="17" /> Instead, here we assume that also $\gamma_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

<footnote n="17">¹⁷This parameterization was adopted to make an analogy to the standard assumption made for coupled quintessence models [118, 119].</footnote>
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 $\gamma_c$ instead of $\beta_\gamma$ also simplifies the numerical implementation in comparison to ref. [39], as the expression for $c_s^2\alpha$ in eq. (2.26), combined with $\zeta$ in eq. (2.6), is directly given in terms of the free functions. To summarize, the relevant degrees of freedom are $\alpha_{B,0}$, $\alpha_{M,0}$, $\gamma_{c,0}$ and $w_{DE}$, which in the following we refer to as EFT parameters.

## 4.3 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 $\Omega_m$, $\Omega_b$, $\Omega_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 late-time 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,<fnref n="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, $\theta_*$. 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,

$$D_M = \int_0^{z_*} \frac{dz}{H(z)}. \qquad (4.3)$$

<footnote n="18">¹⁸The distance $D_M$ is related to the usual angular diameter distance by $D_M = (1+z)D_A$.</footnote>
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 $\Lambda$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 $\Lambda$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.

## 4.4 Choice of fiducial models

We consider the three following cases:

1. $\Lambda$CDM fiducial with EFT parameters,

$$ \alpha_{M,0}^{\text{fid}} = \alpha_{B,0}^{\text{fid}} = \gamma_{c,0}^{\text{fid}} = 0, \quad w_{\text{DE}}^{\text{fid}} = -1. \quad (\Lambda\text{CDM fiducial}) \quad (4.4) $$

Here we fix $H_0 = 67.66 \text{ km s}^{-1} \text{ Mpc}^{-1}$ according to the Planck measurement [32].

2. Modified gravity with large modifications (MGI). In this case, we assume that $H_0 = 73.04 \text{ km s}^{-1} \text{ 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 $\alpha_{M,0}^{\text{fid}}$, $\gamma_{c,0}^{\text{fid}}$ and $w_{\text{DE}}^{\text{fid}}$ should be different from the $\Lambda$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 $\Lambda$CDM values, unless some fine tuning is present. We thus choose fiducial values $\alpha_{M,0}^{\text{fid}}$, $\alpha_{B,0}^{\text{fid}}$ and $\gamma_{c,0}^{\text{fid}}$, which lie between the 1 and $2\sigma$ constraints of the $\Lambda$CDM fiducial model. Then, the value of $w_{\text{DE}}^{\text{fid}}$ is fixed such that the distance to last scattering matches the Planck measurement when assuming $H_0 = 73.04 \text{ km s}^{-1} \text{ Mpc}^{-1}$. In summary, we choose for the EFT parameters

$$ \alpha_{M,0}^{\text{fid}} = -0.1, \quad \alpha_{B,0}^{\text{fid}} = -0.5, \quad \gamma_{c,0}^{\text{fid}} = 0.15, \quad w_{\text{DE}}^{\text{fid}} = -1.41. \quad (\text{MGI fiducial}) \quad (4.5) $$

Note that $w_{\text{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_s^2\alpha > 0$ (see eq. (5.1) below).
<figure><img src="image_3.png" /><figcaption>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.</figcaption></figure>

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

$$ \alpha_{M,0}^{\text{fid}} = -0.01, \quad \alpha_{B,0}^{\text{fid}} = -0.1, \quad \gamma_{c,0}^{\text{fid}} = 0.01, \quad w_{\text{DE}}^{\text{fid}} = -1.01. \quad (\text{MGII fiducial}) \quad (4.6) $$

In this case, $H_0$ is fixed to the Planck value, $\alpha_{M,0}^{\text{fid}}$, $\alpha_{B,0}^{\text{fid}}$ and $\gamma_{c,0}^{\text{fid}}$ are chosen arbitrarily, and $w_{\text{DE}}^{\text{fid}}$ is determined from the distance to last scattering.

We always fix the values of $A_s$ and $n_s$ and the density parameters $\Omega_m$, $\Omega_b$, $\Omega_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 $\omega_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 $\omega_b$, $\omega_c$ as

$$x_c = X\omega_c, \quad x_b = X\omega_b + 1 - X, \qquad (4.7)$$

so that the condition $x_c + x_b = 1$ is always satisfied. The parameter $X$ can have values $X \in [0, 1 + \omega_b/\omega_c]$, corresponding to the extreme cases $x_c = 1$, $x_b = 0$ for $X = 1 + \omega_b/\omega_c$, and $x_c = 0$, $x_b = 1$ for $X = 0$.<fnref n="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.

## 4.5 MCMC analysis

In summary, our model is characterized by the following set of parameters:

$$(\alpha_{B,0}, \alpha_{M,0}, \gamma_{c,0}, w_{\text{DE}}, X, b_{1,B}, b_{2,B}, b_{1,F}, b_{2,B}). \qquad (4.8)$$

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 $\Omega_m$, $\Omega_c$ and $\Omega_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_s^2\alpha$ from eq. (2.26) with $\zeta$ given by eq. (2.6). Note that the derivative $\alpha'_B$ can be expressed using the time evolution specified in eq. (4.2). If $c_s^2\alpha$ 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

<footnote n="19">¹⁹Note that, if the WEP is satisfied, then $\omega_b$, $\omega_c$ are constant, and the domain of $X$ is fixed by the initial conditions. In our case, however, baryons and CDM can evolve differently, so the upper bound for $X$ depends on the background evolution and in particular on the value of $\gamma_c$. In practice, during the MCMC analysis, at each step we reject the sample if $X$ lies outside the range $[0, 1 + \omega_b/\omega_c]$.</footnote>
for the covariance was thoroughly tested.<fnref n="20" /> For the angular diameter distance to the CMB, we assume a Gaussian likelihood centered on the Planck value $D_M = r_*/\theta_*$, where $r_*$ and $\theta_*$ 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 length of the chains is at least 50 times the autocorrelation time, and further excluding a burn-in phase equal to the correlation length and thinning by the autocorrelation time.

# 5 Results

In this section, we present the constraints obtained on the EFT parameters $\alpha_{M,0}$, $\alpha_{B,0}$, $\gamma_{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.

## 5.1 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).

### 5.1.1 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} - \alpha_{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 $\alpha_M$ through the

<footnote n="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.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Constraints around a ΛCDM fiducial model.</figcaption></figure>

<table><caption>**Table 1.** 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.</caption><thead><tr><th rowspan="2"></th><th colspan="2">Median</th><th colspan="3">Width of the 68% credible interval</th></tr><tr><th>No dipole</th><th>Dipole</th><th>No dipole</th><th>Dipole</th><th>Reduction</th></tr></thead><tbody><tr><td>α<sub>M,0</sub></td><td>0.01</td><td>0.01</td><td>[-0.18, 0.31]</td><td>[-0.11, 0.17]</td><td>43%</td></tr><tr><td>α<sub>B,0</sub></td><td>-0.26</td><td>-0.17</td><td>[-0.37, 0.19]</td><td>[-0.24, 0.13]</td><td>34%</td></tr><tr><td>γ<sub>c,0</sub></td><td>0.09</td><td>0.06</td><td>[-0.05, 0.11]</td><td>[-0.04, 0.09]</td><td>19%</td></tr><tr><td>w<sub>DE</sub></td><td>-1.06</td><td>-1.05</td><td>[-0.27, 0.12]</td><td>[-0.16, 0.08]</td><td>38%</td></tr></tbody></table>

evolution of the radiation density $\Omega_r$, see eq. (2.14). This leads to a strong degeneracy between $\alpha_{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 $\alpha_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 $\alpha_{M,0}$. This leads to the two branches in the contours in figure 2. Note that the direction of the $w_{DE} - \alpha_{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.

* $\alpha_{M,0} - \alpha_{B,0}$ and $w_{DE} - \alpha_{B,0}$ contours: we obtain a first branch with $w_{DE} < -1$ and $\alpha_{M} < 0$, in which case the stability of the perturbations requires $\alpha_{B} < 0$. Indeed, linearizing eq. (2.26) around $\Lambda$CDM and applying eq. (2.15), we find that

$$c_s^2\alpha = 2\alpha_M + 3(1+w_{DE})\Omega_{DE} + (1+\Omega_r - 3\Omega_{DE})\alpha_B - 2\alpha'_B. \quad (5.1)$$

Since $1 + \Omega_r - 3\Omega_{DE} < 0$ at low redshift, we see that $\alpha_B < 0$ is needed to keep $c_s^2\alpha > 0$ when the first two terms are negative (note that $\alpha'_B$ has the same sign as $\alpha_B$ with our chosen time evolution). Moreover, since $\Omega_{DE}$ and $\alpha_M$ both decay at high redshift, while the term in front of $\alpha_B$ increases at high redshift (it becomes less negative), a strongly negative $\alpha_{B,0}$ is necessary to maintain $c_s^2\alpha > 0$ at higher redshift.<fnref n="21" /> This is clearly visible in the $\alpha_{M,0} - \alpha_{B,0}$ and $w_{DE} - \alpha_{B,0}$ contours of figure 2.

On the other hand, we have a second branch with $w_{DE} > -1$ and $\alpha_M > 0$, which both contribute positively to the sound speed. Since the term proportional to $\alpha_B$ in eq. (5.1) dominates over the first two at high redshift, $\alpha_{B,0}$ needs to be negative also in this branch to keep $c_s^2\alpha > 0$ in the past. In this case, however, the absolute value of $\alpha_{B,0}$ does not need to be as large as in the other branch, as is visible from figure 2.

* Contours involving the parameter $\gamma_{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

$$\frac{1}{c_s^2\alpha}(\alpha_B - \alpha_M + 3\gamma_c)(\alpha_B - \alpha_M + 3\gamma_c\omega_c b_c). \quad (5.2)$$

This explains the behavior of $\gamma_{c,0}$ along the aforementioned two branches. Along the first branch, $\alpha_{M,0}$ and $\alpha_{B,0}$ are negative, and $|\alpha_{B,0}| > |\alpha_{M,0}|$ (to ensure stability also in the past). To counterbalance this, and keep the growth rate close to $\Lambda$CDM, it is necessary to have $\gamma_{c,0} > 0$. Along the second branch, we have $\alpha_{M,0} > 0$ and large, while $\alpha_{B,0} < 0$, again requiring $\gamma_{c,0} > 0$ to keep the growth rate consistent with $\Lambda$CDM. This is clearly visible in the $\gamma_{c,0} - \alpha_{M,0}$ and $\gamma_{c,0} - \alpha_{B,0}$ contours of figure 2.

The degeneracy between $\gamma_c$ and $\alpha_M$ can be understood from the Einstein and Euler equations: a positive $\gamma_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 $\alpha_M$ in the first branch implies that the mass increases with time, leading to a decrease of the

<footnote n="21">²¹Note that, when dark energy becomes negligible at very high redshift, $1 + \Omega_r - 3\Omega_{DE}$ changes sign. However, $\alpha_B$ and $\alpha_M$ also approach zero at that point and we recover GR.</footnote>
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 $\gamma_c > 0$. In the second branch, a negative $\alpha_B$ counterbalances a positive $\gamma_c$. The parameter $\alpha_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.

### 5.1.2 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

$$ \gamma_c \left[ f - (1+z) \left( \frac{\mathcal{H}_0}{\mathcal{H}} \right)^2 \frac{3\Omega_{m,0}}{c_s^2\alpha} (\alpha_B - \alpha_M + 3\gamma_c\omega_c b_c) \right], \quad (5.3) $$

as can be seen from eq. (3.10). From this, we notice that keeping the dipole close to its $\Lambda$CDM value requires either a small $\gamma_c$ or a cancellation of the two terms in the square bracket. This tightens the constraints on $\gamma_c$, as is clearly visible in figure 2. Moreover, the degeneracy between $\gamma_c$, $\alpha_B$ and $\alpha_M$ is broken by imposing that $(\alpha_B - \alpha_M + 3\gamma_c w_c b_c)/c_s^2\alpha$ 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 $\alpha_M$, the tightening of the constraints on $\alpha_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 $\Delta b = 1$ the cumulative SNR of the dipole in $\Lambda$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\sigma$ contours of figure 2 with $\alpha_{M,0} = 1.3$, $w_{DE} = -0.8$, $\alpha_{B,0} = -0.4$ and $\gamma_{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 $\Omega_m$, $\Omega_b$, $\Omega_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.

### 5.1.3 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 $\alpha_M$ and $\gamma_c$, where the latter governs the interaction of dark energy with dark matter. Alternatively, we can define an effective
<figure><img src="image_5.png" /><figcaption>Figure 3. Constraints around a ΛCDM fiducial model, employing the alternative parameterization with the effective equation of state parameter $\hat{w}_{DE}$.</figcaption></figure>

equation of state, $\hat{w}_{DE}$, directly encoding the evolution of dark energy through

$$ \bar{\rho}'_{DE} = -3(1 + \hat{w}_{DE})\bar{\rho}_{DE}. \qquad (5.4) $$

Comparing eq. (5.4) with eq. (2.11), we see that the effective equation of state is related to the true equation of state by

$$ \hat{w}_{DE} = w_{DE} - \frac{1}{3}\alpha_{M}\frac{\bar{\rho}_{tot}}{\bar{\rho}_{DE}} + \gamma_{c}\frac{\bar{\rho}_{c}}{\bar{\rho}_{DE}}. \qquad (5.5) $$

This leads to the following expression for the Hubble parameter

$$ H(z) = H_0 \exp \left[ \frac{3}{2} \int_{\ln(a(z))}^{0} d\ln a' \left( 1 + \frac{1}{3}\alpha_M + \frac{1}{3}\Omega_r + \hat{w}_{DE}\Omega_{DE} - \gamma_c\Omega_c \right) \right]. \qquad (5.6) $$

In figure 3, we show the constraints on the new set of parameters $\alpha_{M,0}$, $\alpha_{B,0}$, $\gamma_{c,0}$ and $\hat{w}_{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 $\hat{w}_{DE}$ depends on $\gamma_c$ and $\alpha_M$. In particular, we notice that now an equation of state $\hat{w}_{DE} < -1$ is degenerated with a positive $\alpha_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 $\alpha_M$, $\gamma_c$ and $\hat{w}_{DE}$, as can
<table><caption>**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 $\hat{w}_{DE}$. For the credible interval, we also show the percent reduction obtained by the latter.</caption><thead><tr><th rowspan="2"></th><th colspan="2">Median</th><th colspan="3">Width of the 68% credible interval</th></tr><tr><th>No dipole</th><th>Dipole</th><th>No dipole</th><th>Dipole</th><th>Reduction</th></tr></thead><tbody><tr><td>α<sub>M,0</sub></td><td>0.02</td><td>0.00</td><td>[-0.12, 0.47]</td><td>[-0.08, 0.20]</td><td>53%</td></tr><tr><td>α<sub>B,0</sub></td><td>-0.15</td><td>-0.12</td><td>[-0.20, 0.12]</td><td>[-0.19, 0.10]</td><td>9%</td></tr><tr><td>γ<sub>c,0</sub></td><td>0.06</td><td>0.05</td><td>[-0.03, 0.06]</td><td>[-0.03, 0.06]</td><td>0%</td></tr><tr><td>ŵ<sub>DE</sub></td><td>-1.02</td><td>-1.01</td><td>[-0.06, 0.02]</td><td>[-0.03, 0.01]</td><td>50%</td></tr></tbody></table>

be seen from eq. (5.6). From this equation we notice that $\hat{w}_{DE} < -1$ is counterbalanced by $\alpha_M > 0$ and $\hat{w}_{DE} > -1$ by $\alpha_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 $\alpha_{M,0} - \hat{w}_{DE}$ contours, which are shifted towards positive values of $\alpha_{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_s^2\alpha$ around ΛCDM with the new parameterization, we find

$$c_s^2\alpha = 3\alpha_M + 3(1 + \hat{w}_{DE})\Omega_{DE} - (2 - 3\Omega_m - 4\Omega_r)\alpha_B - 2\alpha'_B - 3\gamma_c\Omega_c. \quad (5.7)$$

The degeneracy in the $\alpha_{M,0} - \hat{w}_{DE}$ plane is governed by the distance to last scattering, which is sensitive to the combination $\alpha_M + 3(1 + \hat{w}_{DE})\Omega_{DE}$ through eq. (5.6). Along this degeneracy direction, when $\alpha_M$ is positive, the positivity of the combination $\alpha_M + (1 + \hat{w}_{DE})\Omega_{DE}$ in eq. (5.7) is automatically maintained. This leads to a positive $c_s^2\alpha$ even when $\alpha_B = 0$, as long as $\gamma_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 $\alpha_M$ is negative along the degeneracy direction, then $\alpha_M + (1 + \hat{w}_{DE})\Omega_{DE} < 0$. To ensure stability, it is thus necessary to have negative values of $\alpha_B$ and $\gamma_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 $\alpha_{M,0}$ by breaking the degeneracies with $\alpha_{B,0}$ and $\gamma_{c,0}$, which in turns also tightens the constraints on $\hat{w}_{DE}$. Quantitatively, from table 2, we see that the $1\sigma$ constraints on $\alpha_{M,0}$ and $\hat{w}_{DE}$ are improved by 53% and 50% respectively when the dipole is added. On the other hand, the constraints on $\alpha_{B,0}$ are only mildly improved with this parameterization, while those on $\gamma_{c,0}$ show no improvement. Nevertheless, the red 2-dimensional contours involving $\gamma_{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.
<figure><img src="image_6.png" /><figcaption>Figure 4. Constraints around the modified gravity fiducial MGI.</figcaption></figure>

## 5.2 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 $\gamma_{c,0} = 0.15$ and $\alpha_{M,0} = -0.1$ are well within the vertical branch of the $\gamma_{c,0} - \alpha_{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\sigma$ error bars on $\alpha_{M,0}$, $\alpha_{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 $\gamma_{c,0}$, which is incompatible with zero at more than $2\sigma$. 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
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Constraints given a ΛCDM fiducial when one parameter is kept fixed. Left panel: $\gamma_c$ is fixed to zero. Right panel: $\alpha_B = \alpha_M/2$.</figcaption></figure>

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.

## 5.3 Constraints for specific sub-classes of Horndeski theories

As discussed above, the strength of gravitational redshift is to break the degeneracy between
the parameters $\gamma_{c,0}$, $\alpha_{B,0}$ and $\alpha_{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 $\gamma_c = 0$. In
the right panel we show another specific case, where the parameter $\gamma_{c,0}$ is left free, but
we fix $\alpha_B = \alpha_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-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 $\alpha_{M,0}$ and $w_{DE}$ is reduced, due to
the increased constraints on $\alpha_{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 $\gamma_c$ and $\alpha_M$ and therefore helps breaking the degeneracy between these parameters in the model with $\alpha_B = \alpha_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.

# 6 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 approximation 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.
**Software availability.** The Python code EF-TIGRE (Effective Field Theory of Interacting dark energy with Gravitational Redshift) developed for this project is publicly available on Github at https://github.com/Mik3M4n/EF-TIGRE. The version used for this paper is available on Zenodo at https://zenodo.org/doi/10.5281/zenodo.10606418.

## Acknowledgments

We thank Lawrence Dam for useful discussions and Goran Jelic-Cizmek for helpful comments concerning the numerical implementation of the covariance. The numerical analysis was performed on the Baobab computer cluster of the University of Geneva. S.C., N.G., I.T. and C.B. acknowledge support from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No. 863929; project title "Testing the law of gravity with novel large-scale structure observables"). D.S.B. and C.B. acknowledge support from the Swiss National Science Foundation. The work of M.M. is supported by European Union's H2020 ERC Starting Grant No. 945155-GWmining, Cariplo Foundation Grant No. 2021-0555, the ICSC National Research Centre funded by NextGenerationEU, and MIUR PRIN Grant No. 2022-Z9X4XS.

## A Relation to the phenomenological parameters

The effective theory of interacting dark energy adopted in this work only contains a limited number of free functions ($\alpha_M$, $\alpha_B$, $\gamma_c$ and $w_{DE}$), which all have a simple theoretical interpretation and enter at the level of the Lagrangian. However, another common approach to investigate the impact of modified gravity and dark energy on cosmological observables consists in introducing some phenomenologically-motivated free functions, entering directly in the Poisson and Euler equations (see e.g. refs. [128-134]). In this appendix, we state the relation between the two sets of parameters arising from these different approaches.

To parameterize modifications in the Poisson equation, the gravitational constant $G$ is usually multiplied by a free function $\mu$, as done in eq. (2.21). Indeed, $\mu$ can be expressed using the free functions $\alpha_M$, $\alpha_B$ and $\gamma_c$, see eq. (2.25). In the phenomenological approach, the difference between the two potentials $\Phi$ and $\Psi$ is expressed by the gravitational slip $\eta$,

$$ \eta = \Phi/\Psi. \qquad (A.1) $$

The gravitational slip is obtained by combining eqs. (4.26) and (4.27) of ref. [62] (note that their $\Phi$ and $\Psi$ are interchanged compared to our notation), and translates into

$$ \eta = \frac{c_s^2\alpha + 2\alpha_B(\alpha_B - \alpha_M + 3\gamma_c\omega_c b_c)}{c_s^2\alpha + 2(\alpha_B - \alpha_M)(\alpha_B - \alpha_M + 3\gamma_c\omega_c b_c)}. \qquad (A.2) $$

Apart from the phenomenologically motivated functions $\mu$ and $\eta$, ref. [92] has introduced the free functions $\Theta$ and $\Gamma$ encompassing a breaking of the WEP for CDM. They appear as a friction and fifth force term in the Euler equation,

$$ V'_{c} + (1 + \Theta) V_{c} - \frac{k}{\mathcal{H}}(1 + \Gamma)\Psi = 0. \qquad (A.3) $$
Comparing this with the Euler equation of the effective theory of interacting dark energy given in eq. (2.20), we obtain

$$ \Theta = 3\gamma_c, \qquad (A.4) $$

and

$$ \Gamma = \frac{6\gamma_c(\alpha_B - \alpha_M + 3\gamma_c\omega_c b_c)}{c_s^2\alpha + 2(\alpha_B - \alpha_M)(\alpha_B - \alpha_M + 3\gamma_c\omega_c b_c)}. \qquad (A.5) $$

Note that these phenomenological modifications enter the perturbation equations as follows (see ref. [64]),

$$ \delta_b'' + (2 + \zeta)\delta_b' - \frac{3}{2}\Omega_m\delta_m\mu = 0, \qquad (A.6) $$

$$ \delta_c'' + (2 + \zeta + 3\gamma_c)\delta_c' - \frac{3}{2}\Omega_m\delta_m\mu(\Gamma + 1) = 0, \qquad (A.7) $$

which, combined with eqs. (2.25), (A.4) and (A.5), indeed recover eqs. (2.22) and (2.23).

Galaxy surveys have provided constraints on the free functions $\mu$ and $\eta$ through gravitational lensing [32, 135, 136] and RSD measurements [88]. The combination of data from the CMB, baryon acoustic oscillations, supernovae, weak lensing and RSD recently allowed for a reconstruction of their redshift evolution [116] as well as their joint redshift and scale dependence [137]. Ref. [64] pointed out that the standard RSD-based constraints on $\mu$ are strongly degenerate with the parameters $\Theta$ and $\Gamma$ describing a breaking of the WEP, and forecasted constraints when taking gravitational redshift into account.

We emphasize that, despite the existence of analytical expressions relating the free functions of the phenomenological and effective theory framework, these two approaches are fundamentally different. First of all, in the phenomenological approach, the free functions $\mu$, $\eta$, $\Gamma$ and $\Theta$ are usually assumed to be independent. However, we see in eqs. (2.25), (A.2), (A.4) and (A.5) that they in fact display dependencies through more fundamental, theory-based free functions. Moreover, while we assume that the theory-based functions evolve according to the dark energy density, see eq. (4.2), their interplay in the equations for the phenomenologically motivated functions leads to a different time dependence. Hence, a time evolution of $\mu$, $\eta$, $\Gamma$ and $\Theta$ proportional to the dark-energy density (as applied e.g. in refs. [32, 64, 88, 135]) cannot be maintained in our effective theory framework.

To summarize, a direct, quantitative comparison between the constraints in our work to those obtained employing the phenomenologically motivated parameters cannot be performed. However, some qualitative results are reflected in both these approaches. Most importantly, the role of gravitational redshift in breaking degeneracies that was pointed out in ref. [64] is also maintained in our effective theory approach, as extensively discussed in section 5.

## B Supplementary plots

Here, we present additional plots, supplementing those presented in section 5. All results in this appendix, except figure 9, are obtained using the $\Lambda$CDM fiducial model. First, the corner plot in figure 6 shows the constraints obtained on the parameters $\alpha_{M,0}$, $\alpha_{B,0}$, $\gamma_{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,
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Constraints around the ΛCDM fiducial when including only the monopole, and combining monopole and dipole.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Constraints around the ΛCDM fiducial when including all multipoles and applying various boost factors to the dipole.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 8. Constraints around a ΛCDM fiducial for the full parameter space.</figcaption></figure>

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.
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Constraints around the modified gravity fiducial MGII. The constraints are very similar to those around ΛCDM, see figure 2.</figcaption></figure>

References

[1] G. Gubitosi, F. Piazza and F. Vernizzi, *The Effective Field Theory of Dark Energy*, JCAP 02 (2013) 032 [arXiv:1210.0201] [INSPIRE].

[2] J.K. Bloomfield, É.É. Flanagan, M. Park and S. Watson, *Dark energy or modified gravity? An effective field theory approach*, JCAP 08 (2013) 010 [arXiv:1211.7054] [INSPIRE].

[3] J. Gleyzes, D. Langlois, F. Piazza and F. Vernizzi, *Essential Building Blocks of Dark Energy*, JCAP 08 (2013) 025 [arXiv:1304.4840] [INSPIRE].

[4] J. Bloomfield, *A Simplified Approach to General Scalar-Tensor Theories*, JCAP 12 (2013) 044 [arXiv:1304.6712] [INSPIRE].

[5] F. Piazza, H. Steigerwald and C. Marinoni, *Phenomenology of dark energy: exploring the space of theories with future redshift surveys*, JCAP 05 (2014) 043 [arXiv:1312.6111] [INSPIRE].

[6] J. Gleyzes, D. Langlois and F. Vernizzi, *A unifying description of dark energy*, Int. J. Mod. Phys. D 23 (2015) 1443010 [arXiv:1411.3712] [INSPIRE].

[7] L.Á. Gergely and S. Tsujikawa, *Effective field theory of modified gravity with two scalar fields: dark energy and dark matter*, Phys. Rev. D 89 (2014) 064059 [arXiv:1402.0553] [INSPIRE].

[8] E. Bellini and I. Sawicki, *Maximal freedom at minimum cost: linear large-scale structure in general modifications of gravity*, JCAP 07 (2014) 050 [arXiv:1404.3713] [INSPIRE].
[9] M. Lagos, T. Baker, P.G. Ferreira and J. Noller, *A general theory of linear cosmological perturbations: scalar-tensor and vector-tensor theories*, JCAP 08 (2016) 007 [arXiv:1604.01396] [INSPIRE].

[10] T. Baker, P.G. Ferreira and C. Skordis, *The Parameterized Post-Friedmann framework for theories of modified gravity: concepts, formalism and examples*, Phys. Rev. D 87 (2013) 024015 [arXiv:1209.2117] [INSPIRE].

[11] P.G. Ferreira, T. Baker and C. Skordis, *Testing general relativity with cosmology: a synopsis of the parametrized post-Friedmann approach*, Gen. Rel. Grav. 46 (2014) 1788 [INSPIRE].

[12] C. Skordis, A. Pourtsidou and E.J. Copeland, *Parametrized post-Friedmannian framework for interacting dark energy theories*, Phys. Rev. D 91 (2015) 083537 [arXiv:1502.07297] [INSPIRE].

[13] G.W. Horndeski, *Second-order scalar-tensor field equations in a four-dimensional space*, Int. J. Theor. Phys. 10 (1974) 363 [INSPIRE].

[14] S.M. Carroll, V. Duvvuri, M. Trodden and M.S. Turner, *Is cosmic speed - up due to new gravitational physics?*, Phys. Rev. D 70 (2004) 043528 [astro-ph/0306438] [INSPIRE].

[15] B. Ratra and P.J.E. Peebles, *Cosmological Consequences of a Rolling Homogeneous Scalar Field*, Phys. Rev. D 37 (1988) 3406 [INSPIRE].

[16] C. Wetterich, *Cosmology and the Fate of Dilatation Symmetry*, Nucl. Phys. B 302 (1988) 668 [arXiv:1711.03844] [INSPIRE].

[17] C. Brans and R.H. Dicke, *Mach's principle and a relativistic theory of gravitation*, Phys. Rev. 124 (1961) 925 [INSPIRE].

[18] C. Deffayet, O. Pujolas, I. Sawicki and A. Vikman, *Imperfect Dark Energy from Kinetic Gravity Braiding*, JCAP 10 (2010) 026 [arXiv:1008.0048] [INSPIRE].

[19] A. Nicolis, R. Rattazzi and E. Trincherini, *The Galileon as a local modification of gravity*, Phys. Rev. D 79 (2009) 064036 [arXiv:0811.2197] [INSPIRE].

[20] C. Deffayet, G. Esposito-Farese and A. Vikman, *Covariant Galileon*, Phys. Rev. D 79 (2009) 084003 [arXiv:0901.1314] [INSPIRE].

[21] C. Deffayet, X. Gao, D.A. Steer and G. Zahariade, *From k-essence to generalised Galileons*, Phys. Rev. D 84 (2011) 064039 [arXiv:1103.3260] [INSPIRE].

[22] J. Gleyzes, D. Langlois, F. Piazza and F. Vernizzi, *Healthy theories beyond Horndeski*, Phys. Rev. Lett. 114 (2015) 211101 [arXiv:1404.6495] [INSPIRE].

[23] J. Gleyzes, D. Langlois, F. Piazza and F. Vernizzi, *Exploring gravitational theories beyond Horndeski*, JCAP 02 (2015) 018 [arXiv:1408.1952] [INSPIRE].

[24] D. Langlois and K. Noui, *Degenerate higher derivative theories beyond Horndeski: evading the Ostrogradski instability*, JCAP 02 (2016) 034 [arXiv:1510.06930] [INSPIRE].

[25] M. Crisostomi, K. Koyama and G. Tasinato, *Extended Scalar-Tensor Theories of Gravity*, JCAP 04 (2016) 044 [arXiv:1602.03119] [INSPIRE].

[26] D. Langlois, M. Mancarella, K. Noui and F. Vernizzi, *Effective Description of Higher-Order Scalar-Tensor Theories*, JCAP 05 (2017) 033 [arXiv:1703.03797] [INSPIRE].

[27] T. Baker et al., *Strong constraints on cosmological gravity from GW170817 and GRB 170817A*, Phys. Rev. Lett. 119 (2017) 251301 [arXiv:1710.06394] [INSPIRE].

[28] P. Creminelli and F. Vernizzi, *Dark Energy after GW170817 and GRB170817A*, Phys. Rev. Lett. 119 (2017) 251302 [arXiv:1710.05877] [INSPIRE].
[29] J. Sakstein and B. Jain, *Implications of the Neutron Star Merger GW170817 for Cosmological Scalar-Tensor Theories*, Phys. Rev. Lett. **119** (2017) 251303 [arXiv:1710.05893] [INSPIRE].

[30] J.M. Ezquiaga and M. Zumalacárregui, *Dark Energy After GW170817: Dead Ends and the Road Ahead*, Phys. Rev. Lett. **119** (2017) 251304 [arXiv:1710.05901] [INSPIRE].

[31] C.D. Kreisch and E. Komatsu, *Cosmological Constraints on Horndeski Gravity in Light of GW170817*, JCAP **12** (2018) 030 [arXiv:1712.02710] [INSPIRE].

[32] PLANCK collaboration, *Planck 2018 results. VI. Cosmological parameters*, Astron. Astrophys. **641** (2020) A6 [Erratum ibid. **652** (2021) C4] [arXiv:1807.06209] [INSPIRE].

[33] J. Noller and A. Nicola, *Cosmological parameter constraints for Horndeski scalar-tensor gravity*, Phys. Rev. D **99** (2019) 103502 [arXiv:1811.12928] [INSPIRE].

[34] J. Noller, *Cosmological constraints on dark energy in light of gravitational wave bounds*, Phys. Rev. D **101** (2020) 063524 [arXiv:2001.05469] [INSPIRE].

[35] E. Bellini, A.J. Cuesta, R. Jimenez and L. Verde, *Constraints on deviations from $\Lambda$CDM within Horndeski gravity*, JCAP **02** (2016) 053 [Erratum ibid. **06** (2016) E01] [arXiv:1509.07816] [INSPIRE].

[36] P. Creminelli, G. Tambalo, F. Vernizzi and V. Yingcharoenrat, *Dark-Energy Instabilities induced by Gravitational Waves*, JCAP **05** (2020) 002 [arXiv:1910.14035] [INSPIRE].

[37] V. Salvatelli, F. Piazza and C. Marinoni, *Constraints on modified gravity from Planck 2015: when the health of your theory makes the difference*, JCAP **09** (2016) 027 [arXiv:1602.08283] [INSPIRE].

[38] LIGO SCIENTIFIC et al. collaborations, *Gravitational Waves and Gamma-rays from a Binary Neutron Star Merger: GW170817 and GRB 170817A*, Astrophys. J. Lett. **848** (2017) L13 [arXiv:1710.05834] [INSPIRE].

[39] J. Gleyzes, D. Langlois, M. Mancarella and F. Vernizzi, *Effective Theory of Dark Energy at Redshift Survey Scales*, JCAP **02** (2016) 056 [arXiv:1509.02191] [INSPIRE].

[40] D. Alonso, E. Bellini, P.G. Ferreira and M. Zumalacárregui, *Observational future of cosmological scalar-tensor theories*, Phys. Rev. D **95** (2017) 063502 [arXiv:1610.09290] [INSPIRE].

[41] A. Spurio Mancini et al., *Testing (modified) gravity with 3D and tomographic cosmic shear*, Mon. Not. Roy. Astron. Soc. **480** (2018) 3725 [arXiv:1801.04251] [INSPIRE].

[42] M. Crisostomi and K. Koyama, *Self-accelerating universe in scalar-tensor theories after GW170817*, Phys. Rev. D **97** (2018) 084004 [arXiv:1712.06556] [INSPIRE].

[43] A. Dima and F. Vernizzi, *Vainshtein Screening in Scalar-Tensor Theories before and after GW170817: Constraints on Theories beyond Horndeski*, Phys. Rev. D **97** (2018) 101302 [arXiv:1712.04731] [INSPIRE].

[44] I.D. Saltas, I. Sawicki and I. Lopes, *White dwarfs and revelations*, JCAP **05** (2018) 028 [arXiv:1803.00541] [INSPIRE].

[45] I.D. Saltas and J. Christensen-Dalsgaard, *Searching for dark energy with the Sun*, Astron. Astrophys. **667** (2022) A115 [arXiv:2205.14134] [INSPIRE].

[46] D. Langlois, R. Saito, D. Yamauchi and K. Noui, *Scalar-tensor theories and modified gravity in the wake of GW170817*, Phys. Rev. D **97** (2018) 061501 [arXiv:1711.07403] [INSPIRE].

[47] T.A. Wagner, S. Schlamminger, J.H. Gundlach and E.G. Adelberger, *Torsion-balance tests of the weak equivalence principle*, Class. Quant. Grav. **29** (2012) 184002 [arXiv:1207.2442] [INSPIRE].
[48] C. van de Bruck, J. Morrice and S. Vu, *Constraints on Nonconformal Couplings from the Properties of the Cosmic Microwave Background Radiation*, Phys. Rev. Lett. **111** (2013) 161302 [arXiv:1303.1773] [INSPIRE].

[49] P. Brax, C. Burrage, A.-C. Davis and G. Gubitosi, *Cosmological Tests of the Disformal Coupling to Radiation*, JCAP **11** (2013) 001 [arXiv:1306.4168] [INSPIRE].

[50] P. Brax and C. Burrage, *Constraining Disformally Coupled Scalar Fields*, Phys. Rev. D **90** (2014) 104009 [arXiv:1407.1861] [INSPIRE].

[51] C.M. Will, *The Confrontation between General Relativity and Experiment*, Living Rev. Rel. **17** (2014) 4 [arXiv:1403.7377] [INSPIRE].

[52] T. Damour, G.W. Gibbons and C. Gundlach, *Dark Matter, Time Varying $G$, and a Dilaton Field*, Phys. Rev. Lett. **64** (1990) 123 [INSPIRE].

[53] J.A. Frieman and B.-A. Gradwohl, *Dark matter and the equivalence principle*, Phys. Rev. Lett. **67** (1991) 2926 [INSPIRE].

[54] M. Archidiacono, E. Castorina, D. Redigolo and E. Salvioni, *Unveiling dark fifth forces with linear cosmology*, JCAP **10** (2022) 074 [arXiv:2204.08484] [INSPIRE].

[55] B.J. Barros, L. Amendola, T. Barreiro and N.J. Nunes, *Coupled quintessence with a $\Lambda$CDM background: removing the $\sigma_8$ tension*, JCAP **01** (2019) 007 [arXiv:1802.09216] [INSPIRE].

[56] H. Desmond and P.G. Ferreira, *Galaxy morphology rules out astrophysically relevant $H$-Sawicki $f(R)$ gravity*, Phys. Rev. D **102** (2020) 104060 [arXiv:2009.08743] [INSPIRE].

[57] B.-A. Gradwohl and J.A. Frieman, *Dark matter, long range forces, and large scale structure*, Astrophys. J. **398** (1992) 407 [INSPIRE].

[58] K. Nordtvedt, *Equivalence Principle for Massive Bodies*. I. Phenomenology*, Phys. Rev. **169** (1968) 1014 [INSPIRE].

[59] L. Hui, A. Nicolis and C. Stubbs, *Equivalence Principle Implications of Modified Gravity Models*, Phys. Rev. D **80** (2009) 104002 [arXiv:0905.2966] [INSPIRE].

[60] S. Clesse and J. García-Bellido, *Seven Hints for Primordial Black Hole Dark Matter*, Phys. Dark Univ. **22** (2018) 137 [arXiv:1711.10458] [INSPIRE].

[61] D. Wands, J. De-Santiago and Y. Wang, *Inhomogeneous vacuum energy*, Class. Quant. Grav. **29** (2012) 145017 [arXiv:1203.6776] [INSPIRE].

[62] J. Gleyzes, D. Langlois, M. Mancarella and F. Vernizzi, *Effective Theory of Interacting Dark Energy*, JCAP **08** (2015) 054 [arXiv:1504.05481] [INSPIRE].

[63] G. D'Amico, Z. Huang, M. Mancarella and F. Vernizzi, *Weakening Gravity on Redshift-Survey Scales with Kinetic Matter Mixing*, JCAP **02** (2017) 014 [arXiv:1609.01272] [INSPIRE].

[64] S. Castello, N. Grimm and C. Bonvin, *Rescuing constraints on modified gravity using gravitational redshift in large-scale structure*, Phys. Rev. D **106** (2022) 083511 [arXiv:2204.11507] [INSPIRE].

[65] C. Bonvin and L. Pogosian, *Modified Einstein versus modified Euler for dark matter*, Nature Astron. **7** (2023) 1127 [arXiv:2209.03614] [INSPIRE].

[66] D. Sobral-Blanco and C. Bonvin, *Measuring anisotropic stress with relativistic effects*, Phys. Rev. D **104** (2021) 063516 [arXiv:2102.05086] [INSPIRE].
[67] D. Sobral-Blanco and C. Bonvin, *Measuring the distortion of time with relativistic effects in large-scale structure*, Mon. Not. Roy. Astron. Soc. **519** (2022) L39 [arXiv:2205.02567] [INSPIRE].

[68] I. Tutusaus, D. Sobral-Blanco and C. Bonvin, *Combining gravitational lensing and gravitational redshift to measure the anisotropic stress with future galaxy surveys*, Phys. Rev. D **107** (2023) 083526 [arXiv:2209.08987] [INSPIRE].

[69] A. Einstein, *The foundation of the general theory of relativity*, Annalen Phys. **49** (1916) 769 [INSPIRE].

[70] R.V. Pound and G.A. Rebka, *Gravitational Red-Shift in Nuclear Resonance*, Phys. Rev. Lett. **3** (1959) 439.

[71] J.C. Lopresto, C. Schrader and A.K. Pierce, *Solar gravitational redshift from the infrared oxygen triplet*, The Astrophysical Journal **376** (1991) 757.

[72] Y.-R. Kim and R.A.C. Croft, *Gravitational redshifts in simulated galaxy clusters*, Astrophys. J. **607** (2004) 164 [astro-ph/0402047] [INSPIRE].

[73] R. Wojtak, S.H. Hansen and J. Hjorth, *Gravitational redshift of galaxies in clusters as predicted by general relativity*, Nature **477** (2011) 567 [arXiv:1109.6571] [INSPIRE].

[74] N. Kaiser, *Measuring Gravitational Redshifts in Galaxy Clusters*, Mon. Not. Roy. Astron. Soc. **435** (2013) 1278 [arXiv:1303.3663] [INSPIRE].

[75] P. Jimeno et al., *Comparing gravitational redshifts of SDSS galaxy clusters with the magnified redshift enhancement of background BOSS galaxies*, Mon. Not. Roy. Astron. Soc. **448** (2015) 1999 [arXiv:1410.6050] [INSPIRE].

[76] Y.-C. Cai, N. Kaiser, S. Cole and C. Frenk, *Gravitational redshift and asymmetric redshift-space distortions for stacked clusters*, Mon. Not. Roy. Astron. Soc. **468** (2017) 1981 [arXiv:1609.04864] [INSPIRE].

[77] H. Zhu et al., *Effects of baryons on the gravitational redshift profile of $\Lambda$CDM halos*, arXiv:1901.05615 [INSPIRE].

[78] S. Alam et al., *Relativistic Effects on Galaxy Redshift Samples due to Target Selection*, Mon. Not. Roy. Astron. Soc. **471** (2017) 2077 [arXiv:1709.07856] [INSPIRE].

[79] S. Alam et al., *Relativistic distortions in the large-scale clustering of SDSS-III BOSS CMASS galaxies*, Mon. Not. Roy. Astron. Soc. **470** (2017) 2822 [arXiv:1709.07855] [INSPIRE].

[80] C. Bonvin, L. Hui and E. Gaztanaga, *Asymmetric galaxy correlation functions*, Phys. Rev. D **89** (2014) 083535 [arXiv:1309.1321] [INSPIRE].

[81] R.A.C. Croft, *Gravitational redshifts from large-scale structure*, Mon. Not. Roy. Astron. Soc. **434** (2013) 3008 [arXiv:1304.4124] [INSPIRE].

[82] C. Bonvin, *Isolating relativistic effects in large-scale structure*, Class. Quant. Grav. **31** (2014) 234002 [arXiv:1409.2224] [INSPIRE].

[83] P. McDonald, *Gravitational redshift and other redshift-space distortions of the imaginary part of the power spectrum*, JCAP **11** (2009) 026 [arXiv:0907.5220] [INSPIRE].

[84] J. Yoo, N. Hamaus, U. Seljak and M. Zaldarriaga, *Going beyond the Kaiser redshift-space distortion formula: a full general relativistic account of the effects and their detectability in galaxy clustering*, Phys. Rev. D **86** (2012) 063514 [arXiv:1206.5809] [INSPIRE].

[85] C. Blake et al., *The WiggleZ Dark Energy Survey: the growth rate of cosmic structure since redshift $z=0.9$*, Mon. Not. Roy. Astron. Soc. **415** (2011) 2876 [arXiv:1104.2948] [INSPIRE].
[86] C. Howlett et al., The clustering of the SDSS main galaxy sample — II. Mock galaxy catalogues and a measurement of the growth of structure from redshift space distortions at z = 0.15, Mon. Not. Roy. Astron. Soc. 449 (2015) 848 [arXiv:1409.3238] [INSPIRE].

[87] BOSS collaboration, The clustering of galaxies in the completed SDSS-III Baryon Oscillation Spectroscopic Survey: cosmological analysis of the DR12 galaxy sample, Mon. Not. Roy. Astron. Soc. 470 (2017) 2617 [arXiv:1607.03155] [INSPIRE].

[88] EBOSS collaboration, Completed SDSS-IV extended Baryon Oscillation Spectroscopic Survey: Cosmological implications from two decades of spectroscopic surveys at the Apache Point Observatory, Phys. Rev. D 103 (2021) 083533 [arXiv:2007.08991] [INSPIRE].

[89] E. Gaztanaga, C. Bonvin and L. Hui, Measurement of the dipole in the cross-correlation function of galaxies, JCAP 01 (2017) 032 [arXiv:1512.03918] [INSPIRE].

[90] C. Bonvin, L. Hui and E. Gaztanaga, Optimising the measurement of relativistic distortions in large-scale structure, JCAP 08 (2016) 021 [arXiv:1512.03566] [INSPIRE].

[91] F. Beutler and E. Di Dio, Modeling relativistic contributions to the halo power spectrum dipole, JCAP 07 (2020) 048 [arXiv:2004.08014] [INSPIRE].

[92] C. Bonvin and P. Fleury, Testing the equivalence principle on cosmological scales, JCAP 05 (2018) 061 [arXiv:1803.02771] [INSPIRE].

[93] S. Saga, A. Taruya, M.-A. Breton and Y. Rasera, Detectability of the gravitational redshift effect from the asymmetric galaxy clustering, Mon. Not. Roy. Astron. Soc. 511 (2022) 2732 [arXiv:2109.06012] [INSPIRE].

[94] P. Bull, Extending cosmological tests of General Relativity with the Square Kilometre Array, Astrophys. J. 817 (2016) 26 [arXiv:1509.07562] [INSPIRE].

[95] R.L. Arnowitt, S. Deser and C.W. Misner, The Dynamics of general relativity, Gen. Rel. Grav. 40 (2008) 1997 [gr-qc/0405109] [INSPIRE].

[96] C. de Rham and S. Melville, Gravitational Rainbows: LIGO and Dark Energy at its Cutoff, Phys. Rev. Lett. 121 (2018) 221101 [arXiv:1806.09417] [INSPIRE].

[97] T. Baker et al., Testing gravitational wave propagation with multiband detections, JCAP 03 (2023) 044 [arXiv:2209.14398] [INSPIRE].

[98] LISA COSMOLOGY WORKING GROUP collaboration, Measuring the propagation speed of gravitational waves with LISA, JCAP 08 (2022) 031 [arXiv:2203.00566] [INSPIRE].

[99] I. Sawicki and E. Bellini, Limits of quasistatic approximation in modified-gravity cosmologies, Phys. Rev. D 92 (2015) 084061 [arXiv:1503.06831] [INSPIRE].

[100] C. Bonvin and R. Durrer, What galaxy surveys really measure, Phys. Rev. D 84 (2011) 063505 [arXiv:1105.5280] [INSPIRE].

[101] A. Challinor and A. Lewis, The linear power spectrum of observed source number counts, Phys. Rev. D 84 (2011) 043516 [arXiv:1105.5292] [INSPIRE].

[102] J. Yoo, A.L. Fitzpatrick and M. Zaldarriaga, A New Perspective on Galaxy Clustering as a Cosmological Probe: General Relativistic Effects, Phys. Rev. D 80 (2009) 083514 [arXiv:0907.0707] [INSPIRE].

[103] D. Jeong, F. Schmidt and C.M. Hirata, Large-scale clustering of galaxies in general relativity, Phys. Rev. D 85 (2012) 023504 [arXiv:1107.5427] [INSPIRE].

[104] J. Yoo, Relativistic Effect in Galaxy Clustering, Class. Quant. Grav. 31 (2014) 234001 [arXiv:1409.3223] [INSPIRE].
[105] J. Yoo and M. Zaldarriaga, *Beyond the Linear-Order Relativistic Effect in Galaxy Clustering: Second-Order Gauge-Invariant Formalism*, Phys. Rev. D 90 (2014) 023513 [arXiv:1406.4140] [INSPIRE].

[106] E. Di Dio, R. Durrer, G. Marozzi and F. Montanari, *Galaxy number counts to second order and their bispectrum*, JCAP 12 (2014) 017 [Erratum ibid. 06 (2015) E01] [arXiv:1407.0376] [INSPIRE].

[107] D. Bertacca, R. Maartens and C. Clarkson, *Observed galaxy number counts on the lightcone up to second order: I. Main result*, JCAP 09 (2014) 037 [arXiv:1405.4403] [INSPIRE].

[108] N. Kaiser, *Clustering in real space and in redshift space*, Mon. Not. Roy. Astron. Soc. 227 (1987) 1 [INSPIRE].

[109] G. Jelic-Cizmek, F. Lepori, C. Bonvin and R. Durrer, *On the importance of lensing for galaxy clustering in photometric and spectroscopic surveys*, JCAP 04 (2021) 055 [arXiv:2004.12981] [INSPIRE].

[110] EUCLID collaboration, *Euclid preparation. XIX. Impact of magnification on photometric galaxy clustering*, Astron. Astrophys. 662 (2022) A93 [arXiv:2110.05435] [INSPIRE].

[111] C. Bonvin et al., *A case study for measuring the relativistic dipole of a galaxy cross-correlation with the Dark Energy Spectroscopic Instrument*, Mon. Not. Roy. Astron. Soc. 525 (2023) 4611 [arXiv:2306.04213] [INSPIRE].

[112] R. Maartens et al., *Magnification and evolution biases in large-scale structure surveys*, JCAP 12 (2021) 009 [arXiv:2107.13401] [INSPIRE].

[113] P.H.F. Reimberg, F. Bernardeau and C. Pitrou, *Redshift-space distortions with wide angular separations*, JCAP 01 (2016) 048 [arXiv:1506.06596] [INSPIRE].

[114] F. Beutler, E. Castorina and P. Zhang, *Interpreting measurements of the anisotropic galaxy power spectrum*, JCAP 03 (2019) 040 [arXiv:1810.05051] [INSPIRE].

[115] E. Paillas, Y.-C. Cai, N. Padilla and A.G. Sánchez, *Redshift-space distortions with split densities*, Mon. Not. Roy. Astron. Soc. 505 (2021) 5731 [arXiv:2101.09854] [INSPIRE].

[116] L. Pogosian et al., *Imprints of cosmological tensions in reconstructed gravity*, Nature Astron. 6 (2022) 1484 [arXiv:2107.12992] [INSPIRE].

[117] Z. Huang, *Observational effects of a running Planck mass*, Phys. Rev. D 93 (2016) 043538 [arXiv:1511.02808] [INSPIRE].

[118] L. Amendola, V. Pettorino, C. Quercellini and A. Vollmer, *Testing coupled dark energy with next-generation large-scale observations*, Phys. Rev. D 85 (2012) 103008 [arXiv:1111.1404] [INSPIRE].

[119] PLANCK collaboration, *Planck 2015 results. XIV. Dark energy and modified gravity*, Astron. Astrophys. 594 (2016) A14 [arXiv:1502.01590] [INSPIRE].

[120] A.G. Riess et al., *A Comprehensive Measurement of the Local Value of the Hubble Constant with 1 km s⁻¹ Mpc⁻¹ Uncertainty from the Hubble Space Telescope and the SH0ES Team*, Astrophys. J. Lett. 934 (2022) L7 [arXiv:2112.04510] [INSPIRE].

[121] M. Karamanis and F. Beutler, *Ensemble slice sampling*, Stat. Comput. 31 (2021) 61 [arXiv:2002.06212] [INSPIRE].

[122] M. Karamanis, F. Beutler and J.A. Peacock, *zeus: a python implementation of ensemble slice sampling for efficient Bayesian parameter inference*, Mon. Not. Roy. Astron. Soc. 508 (2021) 3589 [arXiv:2105.03468] [INSPIRE].
[123] I. Sawicki, I.D. Saltas, L. Amendola and M. Kunz, *Consistent perturbations in an imperfect fluid*, JCAP 01 (2013) 004 [arXiv:1208.4855] [INSPIRE].

[124] Y.-S. Song, W. Hu and I. Sawicki, *The Large Scale Structure of $f(R)$ Gravity*, Phys. Rev. D 75 (2007) 044004 [astro-ph/0610532] [INSPIRE].

[125] S.M. Carroll, I. Sawicki, A. Silvestri and M. Trodden, *Modified-Source Gravity and Cosmological Structure Formation*, New J. Phys. 8 (2006) 323 [astro-ph/0607458] [INSPIRE].

[126] D.N. Vollick, *1/R Curvature corrections as the source of the cosmological acceleration*, Phys. Rev. D 68 (2003) 063510 [astro-ph/0306630] [INSPIRE].

[127] C.S. Lorenz, D. Alonso and P.G. Ferreira, *Impact of relativistic effects on cosmological parameter estimation*, Phys. Rev. D 97 (2018) 023537 [arXiv:1710.02477] [INSPIRE].

[128] L. Amendola, M. Kunz and D. Sapone, *Measuring the dark side (with weak lensing)*, JCAP 04 (2008) 013 [arXiv:0704.2421] [INSPIRE].

[129] E. Bertschinger and P. Zukin, *Distinguishing Modified Gravity from Dark Energy*, Phys. Rev. D 78 (2008) 024015 [arXiv:0801.2431] [INSPIRE].

[130] A. Silvestri, L. Pogosian and R.V. Buniy, *Practical approach to cosmological perturbations in modified gravity*, Phys. Rev. D 87 (2013) 104015 [arXiv:1302.1193] [INSPIRE].

[131] Y.-S. Song et al., *Complementarity of Weak Lensing and Peculiar Velocity Measurements in Testing General Relativity*, Phys. Rev. D 84 (2011) 083523 [arXiv:1011.2106] [INSPIRE].

[132] L. Pogosian, A. Silvestri, K. Koyama and G.-B. Zhao, *How to optimally parametrize deviations from General Relativity in the evolution of cosmological perturbations?*, Phys. Rev. D 81 (2010) 104023 [arXiv:1002.2382] [INSPIRE].

[133] L. Pogosian and A. Silvestri, *What can cosmology tell us about gravity? Constraining Horndeski gravity with $\Sigma$ and $\mu$, Phys. Rev. D 94 (2016) 104014* [arXiv:1606.05339] [INSPIRE].

[134] L. Amendola, D. Bettoni, A.M. Pinho and S. Casas, *Measuring gravity at cosmological scales*, Universe 6 (2020) 20 [arXiv:1902.06978] [INSPIRE].

[135] DES collaboration, *Dark Energy Survey Year 1 Results: Constraints on Extended Cosmological Models from Galaxy Clustering and Weak Lensing*, Phys. Rev. D 99 (2019) 123505 [arXiv:1810.02499] [INSPIRE].

[136] DES collaboration, *Dark Energy Survey Year 3 results: Constraints on extensions to $\Lambda$CDM with weak lensing and galaxy clustering*, Phys. Rev. D 107 (2023) 083504 [arXiv:2207.05766] [INSPIRE].

[137] C. Garcia-Quintero, M. Ishak and O. Ning, *Current constraints on deviations from General Relativity using binning in redshift and scale*, JCAP 12 (2020) 018 [arXiv:2010.12519] [INSPIRE].