# Imprints of relativistic effects on the asymmetry of the halo cross-correlation function: from linear to non-linear scales

Michel-Andrès Breton, Yann Rasera, Atsushi Taruya, Osmin Lacombe, Shohei Saga

► To cite this version:

Michel-Andrès Breton, Yann Rasera, Atsushi Taruya, Osmin Lacombe, Shohei Saga. Imprints of relativistic effects on the asymmetry of the halo cross-correlation function: from linear to non-linear scales. Monthly Notices of the Royal Astronomical Society, 2019, 483 (2), pp.2671-2696. 10.1093/mn-ras/sty3206. hal-01758753

HAL Id: hal-01758753

https://hal.science/hal-01758753

Submitted on 1 Mar 2022

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
<figure><img src="image_2.png" /></figure>

# Imprints of relativistic effects on the asymmetry of the halo cross-correlation function: from linear to non-linear scales

Michel-Andrès Breton,<fnref n="1" />★ Yann Rasera,<fnref n="1" />,² Atsushi Taruya,²,³
Osmin Lacombe²,⁴ and Shohei Saga²

¹LUTH, Observatoire de Paris, PSL Research University, CNRS, Université Paris Diderot, Sorbonne Paris Cité 5 place Jules Janssen, F-92195 Meudon, France

²Center for Gravitational Physics, Yukawa Institute for Theoretical Physics, Kyoto University, Kyoto 606-8502, Japan

³Kavli Institute for the Physics and Mathematics of the Universe (WPI), The University of Tokyo Institutes for Advanced Study The University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa, Chiba 277-8583, Japan

⁴Ecole polytechnique, CNRS, F-91128 Palaiseau Cedex, France

Accepted 2018 November 21. Received 2018 November 21; in original form 2018 March 12

## ABSTRACT

The apparent distribution of large-scale structures in the Universe is sensitive to the velocity/potential of the sources as well as the potential along the line of sight through the mapping from real space to redshift space (redshift-space distortions, RSD). Since odd multipoles of the halo cross-correlation function vanish when considering standard Doppler RSD, the dipole is a sensitive probe of relativistic and wide-angle effects. We build a catalogue of ten million haloes (Milky Way size to galaxy-cluster size) from the full-sky light cone of a new 'RayGalGroup-Sims' N-body simulation which covers a volume of (2.625 $h^{-1}$ Gpc)$^3$ with 4096$^3$ particles. Using ray-tracing techniques, we find the null geodesics connecting all the sources to the observer. We then self-consistently derive all the relativistic contributions (in the weak-field approximation) to RSD: Doppler, transverse Doppler, gravitational, lensing and integrated Sachs-Wolfe. It allows us, for the first time, to disentangle all contributions to the dipole from linear to non-linear scales. At large scale, we recover the linear predictions dominated by a contribution from the divergence of neighbouring line of sights. While the linear theory remains a reasonable approximation of the velocity contribution to the dipole at non-linear scales it fails to reproduce the potential contribution below 30–60 $h^{-1}$ Mpc (depending on the halo mass). At scales smaller than ~10 $h^{-1}$ Mpc, the dipole is dominated by the asymmetry caused by the gravitational redshift. The transition between the two regimes is mass dependent as well. We also identify a new non-trivial contribution from the non-linear coupling between potential and velocity terms.

**Key words:** gravitational lensing: weak–methods: numerical–galaxies: distances and redshifts–large-scale structure of Universe.

## 1 INTRODUCTION

Late time structure formation is a non-linear process which is very sensitive to the underlying cosmology. However we do not observe large-scale structures in themselves but rather an image of these objects via messengers (photons, neutrinos, gravitational waves, etc.). Most of our observations come from light, but the information transported by photons is altered during their path from the source to the observer. This leads to several distortions of the image and spectrum of the objects we are interested in. Because of lensing (Schneider, Ehlers & Falco 1992; Bartelmann & Schneider 2001), the angular

position as well as the shape/luminosity of objects can be modified: this is related to the bending of light near local energy overdensities and tidal deformations of light beams respectively. The observed redshift of an object is also perturbed by its proper motion, its gravitational potential, and light propagation in time-varying potentials. As a consequence, the comoving radial distance inferred from redshift (assuming a given homogeneous cosmology) is also perturbed. The apparent distribution of structures is therefore modified by redshift perturbations and lensing: this effect is called redshift-space distortions (RSD<fnref n="1" />; Kaiser 1987; Hamilton 1992). RSD modify the

<footnote n="1">¹We stick to the terminology, RSD, as used widely in the community although we admit that this term is ambiguous. It should be taken as a</footnote>

<footnote>* E-mail: michel-andres.breton@obspm.fr</footnote>
position and properties of objects but also carry relevant cosmological information: studying these effects is a major contemporary challenge. Early RSD studies only took into account peculiar velocities at linear order in the distant-observer approximation. However with the improvement of data precision more subtle effects need to be accounted for: in Szapudi (2004) and Reimberg, Bernardeau & Pitrou (2016) the distant-observer approximation has been relaxed leading to a wide-angle calculation. In Pápai & Szapudi (2008) and Raccanelli et al. (2018) a more sophisticated treatment of Doppler terms has been implemented and, in McDonald (2009) the effect of gravitational redshift has been included leading to an imaginary part of the power spectrum.

However linear theory does not provide a fully satisfactory prediction since small scales are dominated by non-linear clustering and large-scale modes can also be affected by smaller scale modes through mode coupling. In fact, non-linear effects are already visible in real space at $100 h^{-1}$ Mpc scales in Baryon Acoustic Oscillations (BAO; Crocce & Scoccimarro 2008; Taruya et al. 2012; Rasera et al. 2014). Moreover the mapping from real space to redshift space also becomes non-linear. The question of non-linearities has been addressed through semi-analytical approaches: either using different flavours of perturbation theory (Crocce & Scoccimarro 2006; Matsubara 2008; Carlson, White & Padmanabhan 2009; Taruya, Nishimichi & Saito 2010; Taruya et al. 2012; Carlson, Reid & White 2013; Taruya, Nishimichi & Bernardeau 2013), using streaming models (Scoccimarro 2004; Reid & White 2011), Effective Field Theory (Carrasco, Hertzberg & Senatore 2012), or halo model (Tinker et al. 2008). Alternatively, N-body simulations have been performed (Tinker, Weinberg & Zheng 2006; Tinker 2007) to investigate RSD beyond the quasi-linear regime. One important limitation of most of these works is that only standard RSD have been considered (distant-observer approximation and no relativistic effects).

Recently, several authors computed the observed galaxy number count including all relativistic effects at first order in the weak-field approximation (Yoo, Fitzpatrick & Zaldarriaga 2009; Yoo 2010; Bonvin & Durrer 2011; Challinor & Lewis 2011). Given the complexity of the calculation, all the terms were computed using the linear regime approximation and, assuming a linear mapping between real space and redshift space. These works allow to better understand the relative amplitude of all the contributions to the multipoles of the observed galaxy power spectrum (or two-point correlation function) at large linear scales. Measuring the so-called relativistic effects (i.e. beyond standard RSD effects) in a galaxy survey would be exciting since it would provide alternative ways of testing the nature of gravity and of the dark sector. Unfortunately the (usually studied) even multipoles of the observed halo correlation function are dominated by standard RSD and the detection of such effects might be challenging.

While standard RSD generate only even multipoles, the relativistic and wide-angle effects generate an asymmetry in the observed galaxy distribution (i.e. odd multipoles) when cross-correlating two halo populations living in different environments (or more, Bonvin, Hui & Gaztañaga 2016). Using a multipole expansion of the linear cross-correlation function including all relativistic terms, Bonvin, Hui & Gaztañaga (2014) have shown that the dipole is dominated by relativistic terms (which scale as $H/k$, where $H$ is the conformal Hubble constant and $k$ is the wavenumber of interest). Without confirmation by simulation, the validity of these results was however

limited to large linear scales (>100 Mpc). At these scales the dipole generated by the gravitational potential is cancelled out because of velocity terms present in the Euler equation.

On the other hand at halo scale (< few Mpc), the asymmetry of the distribution of galaxies caused by the gravitational potential of galaxy clusters has been investigated through analytical models (Cappi 1995; Kaiser 2013; Zhao, Peacock & Li 2013) and simulations (Cai et al. 2017). In these studies the relative shift between the mean redshift of two galaxy populations was considered instead of the dipole of the galaxy cross-correlation. Wojtak, Hansen & Hjorth (2011) claimed a detection of this effect by stacking galaxy clusters: this has provided an alternative way to test gravity in cluster although the exact interpretation of the measurement is debated (Kaiser 2013; Zhao et al. 2013).

Croft (2013) proposed to use the same concept at larger scales in order to probe the gravitational redshift outside galaxy clusters. They introduce a new estimator since a clear boundary cannot be defined in the universe as is the case in galaxy clusters. The shell estimator measures the relative shift of galaxy's redshift within spherical shells centred on galaxy clusters. In a recent paper, they have measured this estimator from snapshots of N-body simulations (Zhu et al. 2017). Because of the noise related to the limited size of the simulation, they have used an artificial boost factor to increase the signal-to-noise ratio. Interestingly, Alam et al. (2017) have claimed a detection of this effect within the SDSS survey. However, exact predictions for this estimator remain difficult at all scales (Giusarma et al. 2017). While predictions of the dipole at large linear scales are already well established, a clean dipole measurement from linear to non-linear scales within simulations (or observations) is still missing (Alam et al. 2017; Gaztañaga, Bonvin & Hui 2017; Zhu et al. 2017).

In this paper, we directly measure the cross-correlation between two halo populations within the full-sky light cone of a larger and more resolved simulation. We use sophisticated ray-tracing techniques to self-consistently include all relativistic effects at first order in the weak-field approximation. For the first time, we fill in the gap between large-scale linear predictions of the dipole (dominated by the contribution from the divergence of the line of sights) and small-scale non-linear expectations (dominated by the contribution from gravitational redshift). We also decompose all the contributions to the dipole, compare them to linear theory and shed light on to new ones.

The paper is organized as follows. In Section 2 we review the theoretical predictions for dipole in the linear and non-linear regimes. We then present in Section 3 the methodology used to compute halo cross-correlation function from our simulated light cone. In Section 4 we describe our halo catalogues and test our measurements. Finally in Section 5 we show the results of the dipole from linear to non-linear scales.

## 2 THEORY

In this paper, we will consider a well-defined mass-limited collection of haloes within a given cosmological volume. We will not consider observational effects such as selection effects, magnification-bias, absorption/diffusion of light, redshift errors, or the fact that galaxies can be hidden if they are aligned along the line of sight.

### 2.1 Apparent halo overdensity: from real space to redshift space

We consider scalar perturbations of the Friedmann-Lemaître-Robertson-Walker (FLRW) metric in conformal Newtonian gauge.

<footnote>synonymous of 'Observed-Space Distortions' including all distortions induced by the presence of an observer.</footnote>
The metric reads (Ma & Bertschinger 1995)

$$g_{\mu\nu}dx^\mu dx^\nu = a(\eta)^2 \left[ -(1 + 2\psi/c^2)c^2 d\eta^2 + (1 - 2\phi/c^2)\delta_{ij}dx^i dx^j \right], \quad (1)$$

where a is the expansion factor, c is the speed of light, $\delta_{ij}$ is the Kronecker delta, $\psi$ and $\phi$ the Bardeen potentials (Bardeen 1980), x is the comoving position, and $\eta$ the conformal time. Using $k_v k^v = 0$ (where $k_v$ are the components of the wavevector) and the lensing deviation equation we know that the apparent comoving position of a source is (Challinor & Lewis 2011)

$$s = \chi \mathbf{n} + \frac{c}{H} \delta z \mathbf{n} - \mathbf{n} \int_{0}^{\chi} (\phi + \psi) / c^2 d\chi' - \int_{0}^{\chi} (\chi - \chi') \nabla_{\perp} (\phi + \psi) / c^2 d\chi', \quad (2)$$

where $\mathbf{n}$ is a unit vector pointing towards the direction from which the unperturbed photon is coming and $\chi$ is the unperturbed comoving distance of the source. On the right-hand side the first term $\mathbf{x} = \chi \mathbf{n}$ is the unperturbed comoving position of the source, the second term is the distance perturbation along the line of sight due to redshift perturbation $\delta z$. The third term is the (small) Shapiro effect and the last term is the transverse displacement due to lensing.

For the redshift perturbation we will consider the usual first-order terms plus the special relativistic transverse Doppler effect that can be a non-negligible fraction of the gravitational redshift at small scales (Zhao et al. 2013)

$$\delta z = \frac{a_0}{a} \left\{ \frac{\mathbf{v} \cdot \mathbf{n}}{c} - \frac{(\psi - \psi_0)}{c^2} + \frac{1}{2} \left(\frac{v}{c}\right)^2 - \frac{1}{c^2} \int_{\eta_0}^{\eta_1} \frac{\partial(\phi + \psi)}{\partial \eta} d\eta' \right\}, \quad (3)$$

where $\mathbf{v}$ is the velocity and $\phi$ the potential. Quantities with the subscript '0' are evaluated at the observer location today. In the above expression, we have assumed a comoving observer. Assuming mass conservation gives

$$(1 + \delta^{(s)}(s)) \, d^3 s = (1 + \delta^{(x)}(x)) \, d^3 x, \quad (4)$$

where $\delta^{(s)}$ and $\delta^{(x)}$ are the matter density contrast in redshift space and real space, respectively. We have

$$\delta^{(s)} = (1 + \delta^{(x)}) |J|^{-1} - 1, \quad (5)$$

where $J$ is the Jacobian of the transformation from real space to redshift space.

## 2.2 Two-point halo-halo cross-correlation function: linear theory

The halo-halo cross-correlation function between two halo populations $h_1$ and $h_2$ is given by

$$\xi_{h_1 h_2}(\vec{r}_2 - \vec{r}_1) = \langle \delta_{h_1}(\vec{r}_1) \delta_{h_2}(\vec{r}_2) \rangle, \quad (6)$$

with $\delta_{h_i}(\vec{r}_i)$ the overdensity of population $i$ and $\langle \rangle$ the ensemble average. The cross-correlation function is related to the cross-power spectrum through

$$\langle \delta_{h_1}^{(s)}(\mathbf{k}_1) \delta_{h_2}^{(s)*}(\mathbf{k}_2) \rangle = \int d^3 \mathbf{s}_1 d^3 \mathbf{s}_2 \langle e^{-i\mathbf{k}_1 \cdot \mathbf{s}_1} e^{i\mathbf{k}_2 \cdot \mathbf{s}_2} \delta_{h_1}^{(s)}(\mathbf{s}_1) \delta_{h_2}^{(s)}(\mathbf{s}_2) \rangle, \quad (7)$$

where $\delta^{(s)}(\mathbf{k}_i)$ is the Fourier transform of $\delta^{(s)}(\mathbf{s}_i)$. To rewrite this expression in terms of real space quantities we can use equations (4) and (5). Equation (7) is the general formula for the power spectrum in redshift space but this leads to complicated mode couplings

(Zaroubi & Hoffman 1996). In the linear regime, we can linearize the mapping between real and redshift space,

$$\Delta^{(s)} = \Delta^{(r)} + 1 - |J|, \quad (8)$$

where we use $\Delta$ to denote the galaxy number count as an observable thus gauge invariant quantity. Assuming no velocity bias, the observed galaxy number count is given by the sum of the following terms (Bonvin & Durrer 2011; Challinor & Lewis 2011; Bonvin et al. 2014; Tansella et al. 2018)

$$\Delta^{\text{std}} = b\delta - \frac{1}{H} \nabla_r (\mathbf{v} \cdot \mathbf{n}), \quad (9)$$

$$\Delta^{\text{acc}} = \frac{1}{Hc} \dot{\mathbf{v}} \cdot \mathbf{n}, \quad (10)$$

$$\Delta^{\text{q}} = -\frac{\dot{H}}{cH^2} \mathbf{v} \cdot \mathbf{n}, \quad (11)$$

$$\Delta^{\text{div}} = -\frac{2}{H\chi} \mathbf{v} \cdot \mathbf{n}, \quad (12)$$

$$\Delta^{\text{pot, (1)}} = \frac{1}{Hc} \nabla_r \psi \cdot \mathbf{n}, \quad (13)$$

$$\Delta^{\text{pot, (2)}} = \left( \frac{\dot{H}}{H^2} + \frac{2c}{H\chi} \right) \psi/c^2 - \frac{1}{Hc^2} \dot{\psi}, \quad (14)$$

$$\Delta^{\text{shapiro}} = (\phi + \psi)/c^2, \quad (15)$$

$$\Delta^{\text{lens}} = -\frac{1}{c^2} \int_0^\chi \frac{(\chi - \chi')\chi'}{\chi} \nabla_\perp^2 (\phi + \psi) d\chi', \quad (16)$$

$$\Delta^{\text{isw}} = \frac{1}{Hc^2}(\phi + \psi), \quad (17)$$

$$\Delta^{\text{LC}} = \mathbf{v} \cdot \mathbf{n}/c, \quad (18)$$

with $\delta$ the matter density contrast and $b$ a scale-independent bias. $\Delta^{\text{std}}$ is the standard contribution to RSD (Kaiser 1987), $\Delta^{\text{acc}}$ the contribution from the acceleration of sources, $\Delta^{\text{q}}$ the contribution related to the acceleration of the expansion of the universe, $\Delta^{\text{div}}$ the contribution from the divergence of line of sights due to a finite observer, $\Delta^{\text{pot, (1)}}$ the contribution from the gravitational redshift at first order in $H/k$, $\Delta^{\text{pot, (2)}}$ the contribution of the dominant terms in $(H/k)^2$ to the gravitational redshift, $\Delta^{\text{shapiro}}$ the contribution from the Shapiro time delay, $\Delta^{\text{isw}}$ the contribution from the Integrated Sachs-Wolfe effect, $\Delta^{\text{lens}}$ the lensing contribution equal to the lensing convergence as light-beam deformations modify the apparent source distribution, and $\Delta^{\text{LC}}$ the light-cone contribution as the observed position of sources on the light cone is different from their position on constant-time hypersurfaces due to peculiar velocities (Bonvin et al. 2014). A more refined calculation of this effect is given in Kaiser (2013). We neglect the following terms, which are the subdominant $(H/k)^2$ terms:

$$\Delta_{\text{neglect}} = \left( \frac{\dot{H}}{H^2} + \frac{2c}{H\chi} \right) \frac{1}{c^2} \int_{\eta_0}^{\eta_1} \frac{\partial(\phi + \psi)}{\partial \eta} d\eta' \\ + \frac{2}{\chi c^2} \int_0^{\chi} (\phi + \psi) d\chi'. \quad (19)$$
It is straightforward to see that these terms are subdominant: for the dipole the only contribution of integrated terms comes from the integration between $\chi - r/2$ and $\chi + r/2$ where $\chi$ is the pair centre and $r$ the pair separation for the correlation function. We do not consider the term from transverse Doppler effect because it does not contribute to the correlation function in the linear regime. We also neglect higher order terms in redshift (see Ben-Dayan et al. 2012; Umeh, Clarkson & Maartens 2014). The linear correlation function of two different populations of galaxies is given by $\xi = \langle \Delta_1^S \Delta_2^S \rangle$. In the following we will focus on the terms that generate an asymmetry

$$ \xi^A = \sum_i \langle \Delta_1^{\text{std}} \Delta_2^A \rangle + (1) \leftrightarrow (2), \qquad (20) $$

where $\Delta^A$ is given by equations (9)–(18). In the following, we will omit the subscript which indicates the halo population and we will implicitly assume that we perform the summation over the two populations (first and second term in the right-hand side). Odd multipoles come from an asymmetry in the correlation function via exchange of pairs. If the position of each object of a pair is given by $\mathbf{x}_1$ and $\mathbf{x}_2$ we choose the convention $\mathbf{x} = (\mathbf{x}_1 + \mathbf{x}_2)/2$, $\mathbf{r} = \mathbf{x}_2 - \mathbf{x}_1$, and $\mu = \hat{\mathbf{x}} \cdot \hat{\mathbf{r}}$, where a hat denotes a unit vector. The angle defined this way is symmetric under exchange of pairs and we therefore do not need any additional geometrical term for the dipole due to the choice of angle (Reimberg et al. 2016; Gaztanaga et al. 2017).

The contributions to the dipole are the following (Bonvin et al. 2014; Tansella et al. 2018):

$$ \langle \xi^{\text{acc}} \rangle = - \left( f^2 + \frac{\dot{f}}{\mathcal{H}} + \frac{\dot{\mathcal{H}}}{\mathcal{H}^2} f \right) \mathcal{G} \gamma_1^1(r), \qquad (21) $$

$$ \langle \xi^q \rangle = \frac{\dot{\mathcal{H}}}{\mathcal{H}^2} f \mathcal{G} \gamma_1^1(r), \qquad (22) $$

$$ \langle \xi^{\text{div}} \rangle = \frac{2c}{\mathcal{H}\chi} f \mathcal{G} \gamma_1^1(r), \qquad (23) $$

$$ \langle \xi^{\text{pot.}(1)} \rangle = \frac{3}{2} \Omega_m(z) \mathcal{G} \gamma_1^1(r), \qquad (24) $$

$$ \langle \xi^{\text{pot.}(2)} \rangle = -\frac{3}{2} \Omega_m(z) \frac{\mathcal{H}}{\mathcal{H}_0} \mathcal{G} \gamma_1^2(r) \left( \frac{2c}{\mathcal{H}\chi} + \frac{\dot{\mathcal{H}}}{\mathcal{H}^2} - f + 1 \right), \qquad (25) $$

$$ \langle \xi^{\text{shapiro}} \rangle = -3\Omega_m(z) \frac{\mathcal{H}}{\mathcal{H}_0} \mathcal{G} \gamma_1^2(r), \qquad (26) $$

$$ \langle \xi^{\text{isw}} \rangle = 3\Omega_m(z) \frac{\mathcal{H}}{\mathcal{H}_0} \mathcal{G} \gamma_1^2(r) (1 - f), \qquad (27) $$

where we have $f = \partial \ln D_+ / \partial \ln a$ the linear growth rate, $D_+$ the linear growth factor, $\mathcal{G} = (b_1 - b_2) \frac{\mathcal{H}}{\mathcal{H}_0}$, and

$$ \gamma_{\ell}^{m}(r) = \frac{1}{2\pi^2} \left(\frac{\mathcal{H}_0}{c}\right)^{m} \int dk k^{2-m} j_{\ell}(kr) P(k, z). \qquad (28) $$

with $j_\ell$ the $\ell$th spherical Bessel function and $P(k, z)$ the linear power spectrum at redshift $z$. A dot denotes a derivative w.r.t. conformal time The wide-angle term, coming from the fact that two haloes have different line of sights, reads

$$ \langle \xi^{\text{wa}} \rangle = -\frac{2}{5}(b_1 - b_2) f \frac{r}{\chi} \gamma_2^0(r). \qquad (29) $$

The lensing term is given by (Matsubara 2000; Hui, Gaztanaga & Loverde 2007, 2008)

$$ \langle \xi^{\text{lens}} \rangle = \langle \Delta^{\text{std}} \Delta^{\text{lens}} \rangle = -\frac{9}{4} \Omega_m(z) \frac{r\mathcal{H}}{c} \mathcal{G} \sigma(r, z), \qquad (30) $$

with

$$ \begin{aligned} \sigma(r, z) = & \frac{1}{2\pi} \left( \frac{\mathcal{H}_0}{c} \right) \int_{-1}^{1} \mu^2 d\mu \int dk k J_0 \\ & \times \left( kr \sqrt{1 - \mu^2} \right) P(k, z). \end{aligned} \qquad (31) $$

$J_0$ is the Bessel function of the first kind of order 0. The lensing effect is due to the fact that galaxies on the far side of a massive halo look more spread out than the ones in front of the halo due to light bending. This generates an observed underdense region on the far side of the halo which leads to a negative dipole. Following Bonvin et al. (2014), the light-cone term due to peculiar velocities is given by

$$ \langle \xi^{\text{LC}} \rangle = -f \mathcal{G} \gamma_1^1(r). \qquad (32) $$

The bias and growth factor are not constant. This leads to evolution terms

$$ \langle \xi^{\text{evol}} \rangle = \frac{r}{6} \left\{ \left[ (b_1 - b_2) f' - f(b'_1 - b'_2) \right] \left[ \gamma_0^0(r) - \frac{4}{5} \gamma_2^0(r) \right] \right\}, \qquad (33) $$

$$ \langle \xi^{\text{evol2}} \rangle = \frac{r}{2} (b_1 b'_2 - b'_1 b_2) \gamma_0^0(r), \qquad (34) $$

where an apostrophe denotes a derivative w.r.t. comoving distance. We find it convenient to split the evolution term in two parts because they appear in different configurations (depending on the velocity field) as we will see in Section 5. In computing equations (35) and (36), while we use the bias parameter directly measured from our simulation, the comoving distance derivative of the bias is estimated from theoretical model, since the measured comoving distance derivative is basically noisy (see Fig. A1). Here, we specifically use the model given by Sheth & Tormen (1999)

$$ b_{\text{ST}} = 1 + \frac{dv^2 - 1}{\delta_c} + 2 \frac{p/\delta_c}{1 + (dv^2)^p}, \qquad (35) $$

where $p = 0.3$, $d = 0.707$, $\delta_c = 1.673$ for the $\Lambda$CDM cosmology that we used, $v = \delta_c / \sigma(M, z)$, and $\sigma(M, z) = \sigma(M, 0)D_+(z)$. Finally we get

$$ b'_{\text{ST}} = \frac{2\mathcal{H}f}{c\delta_c} \left[ dv^2 - 2 \frac{p^2(dv^2)^p}{(1 + (dv^2)^p)^2} \right]. \qquad (36) $$

We verified that this value is in agreement with the derivative computed from our simulation (see Appendix A).

Note that in order to compute the full correlation function it is also possible to use the pressureless Euler equation $\dot{\mathbf{v}} \cdot \mathbf{n} + \mathcal{H}\mathbf{v} \cdot \mathbf{n} + \nabla_r \psi = 0$.

We sum up in Table 1 the linear regime contributions to the dipole predictions used in this paper.

The dominant term for large separation is $\xi^{\text{div}}$. It is related to a geometrical effect due to the divergence of line of sights for an observer at finite distance. However this is not a 'wide-angle' term. Indeed, even if we consider the pairs to be aligned the effect does not vanish. It comes from the fact that an element of volume seen under a given solid angle is perceived as less dense when receding from us and similarly is perceived as denser when coming towards us. This generates an overall positive dipole while the usual wide-angle term produces a negative dipole.
<table><caption>**Table 1.** Dipole prediction table for the linear regime. This table indicates which terms to consider when predicting the dipole for a specific choice of angle (θ or β) and redshift (zᵢ) which are given by equations (44)–(49). A cross shows if a term should be added to the prediction, while a zero indicates that the term should not be taken into account.</caption><thead><tr><th></th><th>ξ<sup>pot</sup>, (1)</th><th>ξ<sup>pot</sup>, (2)</th><th>ξ<sup>acc</sup></th><th>ξ<sup>q</sup></th><th>ξ<sup>div</sup></th><th>ξ<sup>wa</sup></th><th>ξ<sup>LC</sup></th><th>ξ<sup>lsw</sup></th><th>ξ<sup>lens</sup></th><th>ξ<sup>evo1</sup></th><th>ξ<sup>evo2</sup></th></tr></thead><tbody><tr><td>flrw = (β, z<sub>0</sub>)</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>×</td><td>0</td><td>0</td><td>0</td><td>×</td></tr><tr><td>(β, z<sub>1</sub>)-flrw: Potential only (1)</td><td>×</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>(β, z<sub>1</sub>)-flrw: Potential only (2)</td><td>×</td><td>×</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>(β, z<sub>2</sub>)-flrw: Doppler only</td><td>0</td><td>0</td><td>×</td><td>×</td><td>×</td><td>×</td><td>0</td><td>0</td><td>0</td><td>×</td><td>0</td></tr><tr><td>(β, z<sub>3</sub>)-flrw: Transverse Doppler only</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>(β, z<sub>4</sub>)-flrw: ISW only</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>×</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>(θ, z<sub>0</sub>)-flrw: Lensing only</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>×</td><td>0</td><td>0</td></tr><tr><td>(θ, z<sub>5</sub>)-flrw: All<sup>2</sup></td><td>×</td><td>×</td><td>×</td><td>×</td><td>×</td><td>×</td><td>0</td><td>×</td><td>×</td><td>×</td><td>0</td></tr></tbody></table>

To infer the final prediction, we computed the linear theory prediction at 200 redshifts between the limits of our survey and took the volume average. Indeed, for some of the dipole terms, it is not equivalent to computing the prediction at the volume averaged redshift. Throughout the paper we always consider the case $b_1 > b_2$.

## 2.3 Two-point halo-halo cross-correlation function: non-linear regime

Starting from equation (6), two difficulties arise in the non-linear regime of structure formation. The evolution of the matter fields (density, velocity, and potential) becomes non-linear. Moreover, the mapping from real space to redshift space becomes non-linear too. A vast literature has addressed these questions in the context of the standard RSD with distant observer approximation (i.e. Doppler effect along one fixed direction). A naive perturbation theory expansion provides poor results because of Finger-of-God effect. Taruya et al. (2010) have developed a perturbation-based theory where the damping is characterized by a univariate function with one single free parameter. This model performs well in the quasi-linear regime. It however needs to be extended to include wide-angle effect and other relativistic effects. We plan to work on these aspects in the future. Another approach to the non-linear regime is to rewrite RSD in the context of the streaming model. Following Scoccimarro (2004) the observed (Doppler-only) correlation function $\xi_s^{\text{std}}$ is given (in the distant observer approximation) by

$$1 + \xi_s^{\text{std}}(s_{\perp}, s_{\|}) = \int [1 + \xi(r_{\perp}, r_{\|})] \mathcal{P}(v_{\|}^{12}|r_{\perp}, r_{\|}) dv_{\|}^{12}, \quad (37)$$

where the apparent position in redshift space is decomposed into a component perpendicular to the line of sight ($s_{\perp}$) and along the line of sight ($s_{\|}$), the real-space position is also decomposed as $r_{\perp} = s_{\perp}$ and $r_{\|} = s_{\|} - v_{\|}^2/c$, $v_{\|}^2$ is the pairwise velocity along the line of sight and $\mathcal{P}(v_{\|}^{12}|r_{\perp}, r_{\|})$ is the pairwise-velocity Probability Distribution Function (PDF) at the position ($r_{\perp}, r_{\|}$) in real space. This decomposition is exact even at non-linear scales. However one still needs to predict the PDF of the pairwise velocity using halo model or perturbation theory. Again an extension to wide angles is still missing as it is quite challenging. The detailed analytical modelling of the non-linear velocity PDF which is highly non-Gaussian (Scoccimarro 2004) (plus possible wide-angle effects at intermediate scales) is beyond the scope of this paper and we leave it for future work. However, as we will see the contribution from gravitational redshift dominates the dipole at small non-linear scales: we have therefore focused on this specific contribution. Taking into account the PDF of the gravitational potential at a given pair separation $\mathcal{P}(\phi_{12}|r_{\perp}, r_{\|})$, a general expression to the halo-halo correlation

function with potential-only RSD is given by

$$1 + \xi_s^{\text{pot}}(s_{\perp}, s_{\|}) = \int [1 + \xi(r_{\perp}, r_{\|})] \mathcal{P}(\phi_{12}|r_{\perp}, r_{\|}) dv_{12}, \quad (38)$$

where $r_{\perp} = s_{\perp}$, $r_{\|} = s_{\|} + \phi_{12}/(cH)$ and $\phi_{12}^2$ is the difference of potential between the two halo populations. We assume a simple spherical model to derive the potential difference as a function of radius. Following Croft (2013), the PDF of the potential difference is a single-valued function which depends only on pair separation. The contribution from the potential difference $\phi_{12}$ to the halo-matter correlation function (no velocity), is given by

$$M_{12}(< r) = 4\pi \bar{\rho} \int_{0}^{r} (\xi_{h1m}(x) - \xi_{h2m}(x)) x^2 dx, \quad (39)$$

$$\phi_{12}(R) = -G \int_{0}^{R} \frac{M_{12}(<r)}{r^2} dr, \quad (40)$$

$$\xi_s^{\text{pot,sing}}(s_{\perp}, s_{\|}) = \xi(r_{\perp}, r_{\|}), \quad (41)$$

where $\xi_{hm}$ is the monopole of the halo-matter correlation function, $M_{12}(<r)$ is the enclosed mass, G is the gravitational constant, $R = \sqrt{r_{\perp}^2 + r_{\|}^2}$, $r_{\perp} = s_{\perp}$, and $r_{\|} = s_{\|} + \phi(R)/(cH)$. $\xi_{hm}$ is taken as the maximum between the linear prediction $b\xi_{mm}$ and a spherical NFW profile (Navarro, Frenk & White 1997) with a concentration parameter given by Zhao et al. (2009). This model is an approximation as the distribution of matter is not spherical (Cai et al. 2017), the PDF of the pairwise potential is not single valued and the halo profile can deviate from NFW (Balmès et al. 2014).

Croft (2013) have also taken into account both the standard and (single-valued) potential contributions to the dipole with a simple streaming model. However they have neglected wide-angle effects, other relativistic effects and they have assumed a simple exponential model for the pairwise velocities. To conclude, a full model of the cross-correlation function with all relativistic effects in the non-linear regime is still missing. We will use cosmological simulations to address this regime including all contributions.

## 3 METHODS

The numerical set-up is described in the three following subsections. In the first subsection, we introduce the new large N-body simulation that we have performed. It is part of the RAYGALGROUPSIMS<fnref n="3" /> suite of simulations dedicated to ray-tracing studies (Breton et al.,

<footnote n="3">³Ray-tracing Galaxy Group Simulations, soon available at http://cosmo.ob.spm.fr/</footnote>
in preparation). This simulation gives us Lagrangian and Eulerian
quantities within snapshots (volume at constant times) and light
cones (volume as seen from an observer located at the centre of the
box, a point further away from the centre being seen further back
in time). In the second subsection, we use MAGRATHEA<fnref n="4" /> (Reverdy
2014), a ray-tracing library integrating photon paths using geodesic
equations on the light cones. Searching for the geodesics connecting
a given observer to all detected haloes, we compute the seen angle
and the observed redshift (which are different from the true angle
and comoving redshift provided directly by the simulation). The ray
tracing gives us catalogues of apparent redshift and angular posi-
tions (z, θ, φ). In the last subsection, we introduce the correlation
function estimator.

We would like to point out that in the following we use the Newtonian gauge to interpret the data from *N*-body simulation. In principle this can lead to errors in redshift due to the fact that in this gauge there are relativistic corrections at the horizon scale. To be perfectly rigorous we should interpret the position of particles in the *N*-body gauge (Fidler et al. 2015, 2016). The difference of interpretation comes from the fact that contrary to the Newtonian gauge, the *N*-body gauge leaves the spatial volume unperturbed in a similar way to Newtonian simulations. To account for relativistic corrections in Newtonian gauge we could also apply a time-independent displacement on the particle position (Chisari & Zaldarriaga 2011) which depends on the initial particle distribution. However as shown by Adamek et al. (2016a,b) and Adamek (2017) these gauge effects are small compared to the effects we are interested in, meaning for this work we can safely use the Newtonian gauge to interpret the data from our Newtonian simulation.

**3.1 N-body simulation**

The simulation used in this work consists in a dark matter only sim-
ulation with 4096³ particles within a volume of (2625 h⁻¹ Mpc)³.
The corresponding particle mass-resolution is 1.88 × 10¹⁰h⁻¹ M☉.
The final number of AMR (Adaptive Mesh Refinement) cells is
0.4 trillion and the spatial resolution reaches 5 h⁻¹kpc. We use
CAMB (Lewis, Challinor & Lasenby 2000) to compute the initial
linear matter power spectrum for a ΛCDM cosmology (7-yr WMAP
data, Komatsu et al. 2011) with Hubble parameter h = 0.72, matter
density Ω_m = 0.25733, baryon density Ω_b = 0.04356, radiation
density Ω_r = 8.076 × 10⁻<fnref n="5" />, slope of the primordial power spec-
trum n_s = 0.963, and normalization σ_8 = 0.801. Initial conditions
were generated with a second-order Lagrangian perturbation the-
ory (2LPT) version of MPGRAFIC (Prunet et al. 2008) at a redshift
z_start = 46.

Dark matter particles were evolved with an improved version
of the particle-mesh adaptive-mesh-refinement (PM-AMR) *N*-body
code RAMSES (Teyssier 2002). We borrow the Triangular Shape
Cloud assignment routine from Li et al. (2012) in order to make the
density, potential, and force more isotropic than with the standard
Cloud In Cell assignment.

During the simulation, ~50 snapshots (particles) and three light
cones (particles and gravity cells) with different depths and aper-
tures were written. The light cones were built using the onion-shells
technique (Fosalba et al. 2008, 2015a,b; Teyssier et al. 2009). In this
article, we focus on snapshots between z = 0.5 and z = 0 and the

full-sky light cone up to $z_{\text{max}} = 0.5$. Choosing this maximum redshift ensures that we avoid any replica in the light cone. The light cone consists in ~300 shells (i.e. every coarse time-step): this ensures a good time resolution.

We use pFoF<fnref n="5" /> (Roy, Bouillot & Rasera 2014), a parallel friend-of-
friend algorithm to detect haloes both in snapshots and light cones.
We adopt a standard linking-length of $b = 0.2$ times the mean
interparticle separation and we only pick haloes with more than 100
particles (to guarantee that haloes are sufficiently sampled).

Last, we note that we use a Newtonian simulation and therefore
the two Bardeen potentials are equal.

**3.2 Ray tracing**

The goal of ray tracing is to make a connection between sources and observers. There are many approaches to ray tracing. The most basic approximation (widely used in analytical works) is the Born approximation where deflections and lensing are computed from integral along the undeflected path. Within simulations a commonly used method consists in splitting the universe in several thin lenses orthogonal to the direction of observation (multiple lens formalism, Hilbert et al. 2009). The light is assumed to move in straight lines between lenses while being deflected when crossing the lenses. Lenses are usually flat and plane parallel, thus inducing error for large angles. Moreover the potential is often computed by solving a 2D Poisson equation on the lens (instead of a 3D one). In order to deal with wide-angle effects an onion-shells technique has been implemented (Fosalba et al. 2008, 2015a,b; Teyssier et al. 2009) with spherical lenses. More recently the ray-tracing algorithm has been ported within RAMSES in order to use the accurate 3D potential from the simulation (Barreira et al. 2016). In most of the works, lensing is computed using integral along the light ray while redshift perturbations are computed independently, only considering the Doppler effect.

Our method consists in unifying deflection, lensing, and redshift
calculation by directly solving geodesics equations. While this ap-
proach is usually done in the field of general relativity, it is not
common within cosmological simulations. To our knowledge, this
approach has been introduced by Killedar et al. (2012) using a
fixed grid resolution and limited-size simulation. Here we use MA-
GRATHEA, a hybrid MPI/pthreads C++ 11 ray-tracing library to
propagate photons on null geodesics within the hierarchy of AMR
grids (Reverdy 2014). Using adaptive mesh is crucial to fully re-
solve lensing, potential, and velocity profiles near haloes. We note
that there has been a recent resurgence of interest for a relativis-
tic approach to ray tracing within cosmological simulations (Giblin
et al. 2017 has solved optical scalar equation with a full general rela-
tivistic code at low resolution; Borzyszkowski, Bertacca & Porciani
2017 has solved it at higher resolution but using Born approxima-
tion).

A detailed description of our fast and very accurate ray-tracing li-
brary based on template-meta-programming is available in Reverdy
(2014), we now review some of its specialities. The light cone pro-
vides a regular grid at coarse level with refinement in high-density
areas. In each cell of the grid we have (a, φ, ∂φ/∂x, ∂φ/∂y, ∂φ/∂z), respec-
tively, the scale factor of the shell, the potential and derivatives of
the potential with respect to spatial Cartesian coordinates of the
simulation box. In order to propagate light rays we interpolate these
quantities at the photon space-time location and solve the geodesic

<footnote n="4">⁴https://github.com/vreverdy/magrathea-pathfinder/ The version used for this work is currently very different from the master branch, however they will be synchronized soon</footnote>

<footnote n="5">⁵https://gitlab obsm pr/roy/pFoF</footnote>
equations. In this work, we chose the linear interpolation for space and nearest-neighbour interpolation for time. Higher order interpolation are interesting prospects for future work, but for this paper we will not consider it. Photons are launched backward in time (backward ray tracing) from the light-cone centre at $z = 0$ with a given $k^i$ (spatial part of the wavevector) and initially setting $g_{\mu\nu}dx^\mu dx^\nu = 0$. We then let MAGRATHEA solve the linearized geodesic equations with the metric given by equation (1),

$$ \frac{d^2\eta}{d\lambda^2} = -\frac{2a'}{a} \frac{d\eta}{d\lambda} \frac{d\eta}{d\lambda} - \frac{2}{c^2} \frac{d\phi}{d\lambda} \frac{d\eta}{d\lambda} + 2 \frac{\partial\phi}{\partial\eta} \left(\frac{d\eta}{d\lambda}\right)^2 \quad (42) $$

$$ \frac{d^2 x^i}{d\lambda^2} = -\frac{2a'}{a} \frac{d\eta}{d\lambda} \frac{dx^i}{d\lambda} + \frac{2}{c^2} \frac{d\phi}{d\lambda} \frac{dx^i}{d\lambda} - 2 \frac{\partial\phi}{\partial x^i} \left(\frac{d\eta}{d\lambda}\right)^2 \quad (43) $$

where $\lambda$ is the affine parameter along the photon path.

As we are interested in source-averaged observables rather than direction-averaged ones (Kibble & Lieu 2005), we now describe recent modifications of the solver to build a catalogue of sources including all relativistic effects. To find the null geodesic connecting a source and the observer, we launch several photons from the observer to the tentative directions of observation of the source. Then using a root finder, the geodesics intersecting both the source and observer world lines are identified. In practice, we assume that sources are present at any time (as opposed to an event which corresponds to a specific space-time location). Moreover, since sources are moving, we use a nearest-neighbour interpolation for the time location of the source. Because sources are moving slowly and light-ray deflections are small, the sources lie very close to the null FLRW light cone. A refinement would be to linearly interpolate the position of particles between two light cones at different times. Moreover, we only search for one geodesic for each source since we focus on large scales, dominated by the weak-lensing regime: generalization to strong lensing (i.e. multiple geodesics for one source) is straight-forward with enough resolution. We leave these possible refinements for future work.

Let a halo be at location $(X, Y, Z)$ on the light cone. For an observer at the centre of the simulation, the two components of the true angle $\beta$ are (assuming the same convention as for lensing): $\beta_1 = \arctan(Y/X)$, $\beta_2 = \arccos(Z/R)$ where $R$ is the comoving distance $R = \sqrt{X^2 + Y^2 + Z^2}$. We expect the lensing deviation to be small, we thus launch the photon in the direction $\beta$, but the ray is deflected and does not reach the position $(X, Y, Z)$. As shown in Fig. 1, we iterate on the initial launching conditions using a root-finder method (Newton's method in our case) to find the initial angle that minimizes the angle difference between $\beta$ and the position of the photon at same radius. In practice only one or two iterations are needed to get an angle difference lower than $10^{-2}$ arcsec. With this method we know the true angle $\beta$ and the seen angle $\theta$. We can then directly derive the Jacobian matrix $A_{ij} = \frac{\partial \beta_i}{\partial \theta_j}$, hence the distortion matrix (related to lensing). This way of computing the lensing directly from a beam of light rays (instead of integrating Sachs equation) is similar to the ray-bundle approach (Fluke, Webster & Mortlock 1999; Fluke & Lasky 2011) except that the geodesic equations are directly integrated.

To gain comprehension on the various contributions to the total redshift we decompose it as follows:

$$ z_0 = \frac{a_0}{a} - 1, \quad (44) $$

$$ z_1 = \frac{a_0}{a} \left( 1 + \frac{\phi_o}{c^2} - \frac{\phi_s}{c^2} \right) - 1, \quad (45) $$

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Illustration of the geodesic-finder algorithm. Each tentative photon is designated by $\zeta_n$ with n being the number of iterations. The first photon $\zeta_0$ is launched towards the source with an angle $\beta$. The first ray will generally miss the source, we then iterate using Newton's method in order to get a new initial angle. In this example we iterate twice to find a ray close enough to the source at the same radius, the initial angle of $\zeta_2$ is given by $\theta$ and is interpreted as the seen angle.</figcaption></figure>

$$ z_2 = \frac{a_0}{a} \left( 1 + v_s \cdot n/c \right) - 1, \quad (46) $$

$$ z_3 = \frac{a_0}{a} \left( 1 + \frac{1}{2} \left( \frac{v_s}{c} \right)^2 \right) - 1, \quad (47) $$

$$ z_4 = \frac{a_0}{a} \left( 1 - \frac{2}{c^2} \int_{\eta_s}^{\eta_o} \dot{\phi} d\eta \right) - 1, \quad (48) $$

$$ z_5 = \frac{(g_{\mu\nu} k^\mu u^\nu)_s}{(g_{\mu\nu} k^\mu u^\nu)_o} - 1, \quad (49) $$

with $g_{\mu\nu}k^\mu u^\nu = -ack^0 \left( 1 + \phi/c^2 + v \cdot n/c + \frac{1}{2} \left( \frac{v}{c} \right)^2 \right)$ and the observer velocity set to zero.

Each redshift corresponds to a specific contribution. $z_0$ is the redshift directly inferred from the scale factor. However this scale factor is related to the conformal time computed until arriving at the source, using the geodesic equation (42). It therefore implicitly takes into account time delay. $z_1$ only includes the gravitational redshift perturbation, $z_2$ the Doppler perturbation, $z_3$ the transverse Doppler perturbation, and $z_4$ the ISW perturbation. $z_5$ is the exact general relativity redshift computation. It almost corresponds to $z_0$ plus all redshift perturbations above. The ISW effect is hidden in the $k^0$ term, which comes directly from our geodesic integration.

Finally, ray tracing gives us catalogues with $\beta, \theta, A_{ij}$, various redshifts containing each terms of equation (3) and the number of dark matter particles for each halo. In these catalogues all the relativistic effects have been computed in a self-consistent way. These catalogues<fnref n="6" /> will be described in detail in Breton et al. (in preparation).

### 3.3 Estimation of the correlation function

The halo-halo two-point cross-correlation function $\xi_{h_1h_2}(r) = \langle \delta_{h_1}(x) \delta_{h_2}(x+r) \rangle$ is a measure of the excess of probability relative to a Poisson distribution of finding a pair of haloes separated

<footnote n="6">⁶http://cosmo obspm.fr/raygalgroupsim s-relativistic-halo-catalogs</footnote>
by a vector **r**. For a statistically homogeneous and isotropic field the correlation function can be written as $\xi_{h_1h_2}(r)$ since it only depends on the norm **r** of the separation. However the presence of an observer breaks the isotropy and one needs to specify two components of the separation vector **r**, for instance its norm **r** and projection $\mu$ along the line of sight.

To estimate the correlation function we used a modified version of CUTE (Alonso 2012) (a parallel tree-code pair-counting algorithm). It implements an LS estimator (Landy & Szalay 1993; Kerscher, Szapudi & Szalay 2000), which is one of the most commonly used estimator for the correlation function (since its variance is almost Poisson),

$$ \xi_{LS} = \frac{D_1 D_2 - D_1 R_2 - R_1 D_2 + R_1 R_2}{R_1 R_2}, \quad (50) $$

where $D_1, D_2$ refer to different data sets to be cross-correlated while $R_1$ and $R_2$ are the associated random catalogues. Moreover the pair counts are normalized by the total number of pairs in the catalogues. Since we are interested in correlation function anisotropies, we bin in $(r, \mu)$. Once we compute $\xi(r, \mu)$, we deduce the multipoles as

$$ \xi_{\ell}(r) \approx \frac{2\ell + 1}{2} \sum_{\mu=-1}^{1} \xi(r, \mu) \mathcal{L}_{\ell}(\mu) \Delta\mu, \quad (51) $$

where $\mathcal{L}_{\ell}(\mu)$ is the Legendre polynomial of order $\ell$.

We have cross-checked the results of this direct pair-counting method to a grid method. In this method the halo density is estimated on a thin Cartesian grid using a Cloud-in-Cell assignment scheme. The correlation function is then directly computed from its definition as a function of the overdensities of the halo populations. The two methods give very similar results. With the intention of comparing to linear theory we estimated the linear bias $b_i$ for data set $\text{data.H}_i$ (see Appendix A for more details),

$$ b_{100} \approx \sqrt{\frac{\xi_{h_{100}h_{100}}^{\ell=0}}{\xi_{mm}^{\ell=0}}}, \quad (52) $$

$$ b_i \approx b_{100} \frac{\xi_{h_i h_{100}}^{\ell=0}}{\xi_{h_{100} h_{100}}^{\ell=0}}, \quad (53) $$

where $\xi_{hh}^{\ell=0}$ and $\xi_{mm}^{\ell=0}$ are the halo-halo and matter-matter correlation function monopole, respectively. The bias is estimated by fitting a constant to equation (53) for $r$ between 25 and $75 h^{-1}$ Mpc. Below $25 h^{-1}$ Mpc the number of pairs is too low and the correlation function may fluctuate while above $75 h^{-1}$ Mpc the Poisson noise becomes non-negligible

The last point concerns the estimation of statistical errors. Running again the same heavy simulation being much too time consuming, we estimate the variance using the jackknife method, as it is the internal method that minimizes most of the variance for the linear regime according to Norberg et al. (2009). In this paper we compute the jackknife method with 32 re-samplings. We then estimate the variance of the correlation function as follows,

$$ \sigma_{\ell}^{2}(r) = \frac{N-1}{N} \sum_{k=1}^{N} (\xi_{\ell}^{k}(r) - \bar{\xi}_{\ell}(r))^{2}, \quad (54) $$

with $N$ the number of re-samplings, $k$ the sample number, and $\bar{\xi}_{\ell}(r)$ the mean correlation function given by

$$ \bar{\xi}_{\ell}(r) = \frac{1}{N} \sum_{k=1}^{N} \xi_{\ell}^{k}(r). \quad (55) $$

<table><caption>**Table 2.** Summary of the different data sets used: name, number of haloes, range for the number of particles per halo, mean mass, bias at the volume averaged redshift $z = 0.341$, and estimated mean halo concentration (taken from Zhao et al. 2009).</caption><thead><tr><th>Name</th><th>Number of haloes</th><th>Number of part</th><th>Mass (h<sup>-1</sup> M<sub>☉</sub>)</th><th>Bias</th><th>c<sup>200m</sup><sub>zhao</sub></th></tr></thead><tbody><tr><td>data.H<sub>100</sub></td><td>5.4 × 10<sup>6</sup></td><td>100–200</td><td>2.8 × 10<sup>12</sup></td><td>1.08</td><td>8.2</td></tr><tr><td>data.H<sub>200</sub></td><td>3.4 × 10<sup>6</sup></td><td>200–400</td><td>5.6 × 10<sup>12</sup></td><td>1.22</td><td>7.7</td></tr><tr><td>data.H<sub>400</sub></td><td>1.9 × 10<sup>6</sup></td><td>400–800</td><td>1.1 × 10<sup>13</sup></td><td>1.42</td><td>7.1</td></tr><tr><td>data.H<sub>800</sub></td><td>1.0 × 10<sup>6</sup></td><td>800–1600</td><td>2.2 × 10<sup>13</sup></td><td>1.69</td><td>6.6</td></tr><tr><td>data.H<sub>1600</sub></td><td>4.0 × 10<sup>5</sup></td><td>1600–3200</td><td>4.5 × 10<sup>13</sup></td><td>2.07</td><td>6.1</td></tr><tr><td>data.H<sub>3200</sub></td><td>2.0 × 10<sup>5</sup></td><td>3200–6400</td><td>9.0 × 10<sup>13</sup></td><td>2.59</td><td>5.7</td></tr></tbody></table>

It is important to note that the variance estimated with equation (54) is good enough in the linear regime where the density field is Gaussian but becomes much less accurate for smaller scales, in the non-linear regime. In this regime error bars should be taken with caution.

In the linear regime, we note that the theoretical predictions for the cross-correlation dipole are proportional to the bias difference (except for evolution effects). Therefore, normalizing by this quantity should give the same signal for each pair of populations. We take advantage of this feature by using a weighted average of the normalized dipole for all mass combinations to increase our signal-to-noise ratio. In the linear regime, the mean signal is computed as

$$ \xi_1^{\text{lin}} / \Delta b = \frac{\sum_{ij} \frac{\xi_j^{ij}}{b_i - b_j} \frac{1}{\sigma_{ij}^2}}{\sum_{ij} \frac{1}{\sigma_{ij}^2}}, \quad (56) $$

where $b_i$ and $b_j$ are the bias of different halo populations. $\xi_1^{ij}$ is the dipole of the cross-correlation between two halo populations of bias $b_i$ and $b_j$, and $\sigma_{ij}^2$ its variance normalized by the bias difference. The variance of the weighted average dipole is

$$ \sigma^2 = \frac{1}{\sum_{ij} \frac{1}{\sigma_{ij}^2}}. \quad (57) $$

The error bars are probably underestimated due to the lack of independence.

# 4 DATA AND VALIDATION

We now proceed to the presentation and validation of our data sets. In Section 4.1 we introduce the halo catalogues and in Section 4.2 we validate our two-point correlation measurements on the monopole and quadrupole.

## 4.1 Data sets

We consider haloes between $z_{\text{min}} = 0.05$ and $z_{\text{max}} = 0.465$ (hereafter, when we refer to the full light cone we mean the light cone between these two redshifts), leading to a volume of $8.34 (h^{-1} \text{Gpc})^3$. We choose redshift limits which are not too close to the observer to avoid issues when computing angles, and not too close from the edge of our full-sky light cone to avoid edge effects. We also focus on haloes with mass between $1.9 \times 10^{12} h^{-1} \text{M}_\odot$ and $1.2 \times 10^{14} h^{-1} \text{M}_\odot$. The total number of haloes in this volume is $1.2 \times 10^7$ leading to a mean halo density of $n \approx 5 \times 10^{-4} \text{ Mpc}^{-3}$. We divide the halo catalogue in six logarithmic mass bins as shown in Table 2.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Redshift distribution for each halo data set in Table 2. The distributions from least to most massive haloes population are shown in black, purple, blue, green, orange, and red. The distributions are normalized so that the integral is unity.</figcaption></figure>

data. $H_N$ represents catalogues of haloes sampled by a number of dark matter particles between $N$ and $2N$. We cross-correlate all the data sets, with 15 linear bins in $r$ going from 0 to $150\ h^{-1}\ Mpc$ for large scales and eight linear bins from 0 to $32\ h^{-1}\ Mpc$ for smaller scales. We also use 201 bins in $\mu$. For computations on the full light cone we generate a random catalogue for each data set with more than 10 times the number of haloes. The redshift distribution for random catalogues follow the same distribution as the associated data catalogue. To avoid losing information on the clustering along the line of sight, the redshift distribution is smoothed using eight redshift bins as shown in Fig. 2.

Otherwise, if our computation is done with shell averages we use a random catalogue with $10^8$ object with a uniform distribution.

## 4.2 Monopole and quadrupole validation

### 4.2.1 Matter autocorrelation monopole on the light cone

We first check our measurement of the well-known matter autocorrelation monopole. Out of the ≈30 billion particles in the light cone we randomly pick $10^8$ that we ray trace. We expect this catalogue to be representative of the general distribution of dark matter in the simulation. We compute the monopole on our full light cone using $10^9$ particles for the random catalogue with a uniform distribution (since the mean matter density does not evolve with redshift), and we compare to the real-space prediction at the volume averaged redshift (the light-cone effect being negligible for this multipole). For the prediction we use the emulator COSMICEMU (Heitmann et al. 2016) which agrees with our power spectrum computation for different snapshots at per cent level between roughly $k = 0.02\ h\ Mpc^{-1}$ and $k = 2\ h\ Mpc^{-1}$. The result is shown in Fig. 3. The monopole of the matter autocorrelation is in good agreement (better than two per cent) with the emulator in the range $r = 20–120\ h^{-1}\ Mpc$.

### 4.2.2 Halo autocorrelation monopole and quadrupole in redshift space

In this section we want to validate our monopole and quadrupole measurements in redshift space taking only into account the effects of peculiar velocity e.g. the standard RSD effect. We compute the

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Matter monopole autocorrelation computed on the unperturbed FLRW light cone compared with COSMICEMU emulator (Heitmann et al. 2016). Subplot shows the relative difference.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Monopole autocorrelation in redshift space (Doppler term only). Red full lines give the prediction from RegPT + TNS (Taruya et al. 2010, 2012; Taruya et al. 2013). Black circles give the computation on the full light cone.</figcaption></figure>

correlation on the full light cone and the errors bars are estimated with the jackknife method. For the quadrupole in redshift space we have subtracted the real-space quadrupole.

To predict the monopole (Fig. 4) and quadrupole (Fig. 5) in redshift space we use the RegPT + TNS (Taruya et al. 2010, 2012; Taruya et al. 2013) model with the measured linear bias and the parameter characterizing the damping of small-scale clustering $\sigma_v = 5.204\ h^{-1}\ Mpc$, which is estimated from the linear theory assuming that haloes trace dark matter flow. We compute the prediction at the volume averaged redshift. The prediction is supposed to be accurate in the weakly non-linear regime, which roughly corresponds to the scales larger than $30\ h^{-1}\ Mpc$ for this redshift. Note
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Absolute value of the quadrupole autocorrelation in redshift space (Doppler term only). Red full lines give the prediction from RegPT + TNS (Taruya et al. 2010, 2012; Taruya et al. 2013). Black circles give the computation on the full light cone.</figcaption></figure>

that the validity of the prediction relies on the assumption of linear bias and on the distant observer approximation.

We have a very good agreement on the redshift-space monopole and a good agreement for the quadrupole, although we see a small discrepancy at large scales which might be due to the cosmic variance of the measurement results or alternatively to the presence of a finite distance observer. Our procedure is therefore fully validated.

Finally we have also checked that the relative amplitude of effects beyond the standard Doppler RSD are of order of ~$10^{-2}$ (monopole) or ~$10^{-3}$ (quadrupole). Disentangling these effects from the main contribution is therefore a challenge. We now focus on the dipole where the standard Doppler contribution vanishes in the distant-observer limit.

# 5 RESULTS

We are now interested in the dipole. In order to investigate this multipole we will split it into three parts. The first part is the one generated by statistical fluctuations in real space. This is mostly interesting for the comparison to observations. It can be minimized by increasing the statistics (number and density of haloes). The second part is generated when the halo catalogue is projected from real space to the unperturbed FLRW light cone. This part is a small contribution to the dipole which has to be taken into account when making accurate predictions. The third part of particular interest in this paper is the dipole generated by the perturbation of the FLRW light cone (at first order in the weak-field approximation) related to the formation of large-scale structures in the universe.

In what follows, we compare our results with the linear theory predictions presented in Section 2.2, and summarized in Table 1. In the predictions including evolution terms (i.e. equations 33 and 34), the comoving distance derivative of the bias changes the amplitude

of dipole correlation, and the resultant prediction does not simply scale as $\xi_1 \propto \Delta b$. We thus plot the averaged prediction as well as the maximum and minimum of the predictions among possible combination of haloes, filling up the interval with grey shade (Fig. 8, upper-right panel of Fig. 10, and Fig. 11).

In Section 5.1 and 5.2 we focus on the dipole at large scales. At these scales the theoretical predictions are proportional to the bias difference between two halo populations (except for evolution effects). We take full advantage of all the cross-correlations by using them to compute the weighted dipole normalized by the bias difference (see Section 3.3). Each cross-correlation are shown in Appendix C. In Section 5.3 we investigate the dipole at smaller scales where non-linearities arise. For these scales we show the dipole and its mass dependence.

## 5.1 Statistical fluctuations and light-cone effect

In this section, we measure the dipole of the halo distribution using snapshot information (i.e. distribution of haloes at constant time). We consider seven snapshots in the interval between $z = 0.05$ and $z = 0.465$. For each snapshot, we build a shell of size of order ~$300 h^{-1}$ Mpc centred on the comoving distance corresponding to the snapshot redshift. In this way, we are able to compute the mean dipole in real space. We can also compare the dipole computed at constant time to the one computed on the FLRW unperturbed light cone at the same position (see Fig. 6).

### 5.1.1 Dipole generated by statistical fluctuations

After averaging the snapshot dipole computed in each shell (with a weigh given by the volume of each shell), we obtain the mean snapshot dipole. The standard errors on the mean are computed from the seven snapshots and should therefore be taken with caution because of the lack of independence. For a wide range of radii ($20 < r < 150 h^{-1}$ Mpc) the dipole shown in Fig. 7 is roughly compatible with zero within the statistical error bars (except for several points at 100 and $130 h^{-1}$ Mpc where the error bars are very likely underestimated). Moreover the error bars are limited to ~$2 \times 10^{-4}$.

In principle, for a very large volume and for a very large density of haloes, the dipole tends towards zero. For finite volume surveys with finite number of galaxies per unit of volume such a fluctuation might blur the dipole. However, this noise can be minimized by increasing the size of the surveys and the density of pairs of haloes (smaller haloes, more halo populations, see Bonvin et al. 2016). As we will see for our light cone, the noise is below the signal but it can sometimes reach the same order of magnitude as the signal. Increasing the halo statistics by a factor of ~10 should be enough to boost the signal-to-noise ratio. In our simulation, we have simply subtracted this noisy contribution to extract the physical part of the dipole signal.

### 5.1.2 Light-cone and evolution effects

In this section, we compute the dipole within the shells in the light cone and the corresponding shells in the snapshots. By subtracting the two we can extract the so-called light-cone effect (Kaiser 2013; Bonvin et al. 2014). The main contribution to this effect is related to the peculiar motions of haloes: haloes are not at the same position in the snapshot (constant time) and in the FLRW light cone. Another
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Representation of snapshots and light cone. Left-hand panel : snapshot at redshift $z_I$. Middle panel : snapshot at redshift $z_{II}$. Right-hand panel : full-sky light cone (made of ~300 small shells of size ~4 $h^{-1}$ Mpc not represented here). Here $z_{II} > z_I$ and $d_I$ (respectively $d_{II}$) is the FLRW comoving distance at redshift $z_I$ (respectively $z_{II}$). To estimate light-cone effects, we compare large light-cone shells (of size $2\Delta r \sim 300\ h^{-1}$ Mpc) to the equivalent volume in a snapshot at the corresponding homogeneous redshift.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Real-space (snapshot) dipole of the cross-correlation function, normalized by the bias. We use a weighted average of all the cross-correlations as well as a weighted average on seven redshift shells. The dipole is consistent with zero within the statistical error bars.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Light-cone-effect contribution to the dipole of the cross-correlation function, normalized by the bias. We use a weighted average of all the cross-correlations as well as a weighted average on seven redshift shells. The linear prediction given by the first line of Table 1 is shown by the grey filled curve (as the prediction is not completely proportional to the bias difference) while the averaged prediction is shown with black solid line.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 9.** The dipole of the cross-correlation function normalized by the bias induced by the Shapiro effect. The prediction for this effect is shown in Section 2.2. This contribution is very small and will be neglected afterwards.</figcaption></figure>

contribution comes from the evolution effects: haloes are not exactly the same in the snapshot and light cone as they experience merging.

While the light-cone effect has already been taken into account in simulations (Cai et al. 2017; Zhu et al. 2017), our approach is more sophisticated since we have directly built the light cone on the fly during the simulation (at each coarse time-step). In previous work, the light-cone effect was added as a post-processing procedure on top of the snapshots. This approach usually neglects the variation of the potential during the evolution. Moreover, evolution effect are not easily captured while they are a direct outcome of our approach.

As shown in Fig. 8, the light-cone effect (plus evolution effect) is in agreement with the linear expectation (see first line of Table 1) while error bars remain important and the points at small scales seem to depart from the averaged prediction. The normalized dipole is of order $10^{-4}$ at most. As we will see later it is a subdominant contribution to the full dipole. The linear prediction is broad due to evolution effects which are not proportional to the bias difference. For the cross-correlation of the two most massive halo populations the evolution effect is of the same order as the light-cone effect. This shows the importance of an accurate modelling as well as precise bias measurement to disentangle these two effects.
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Dipole of the cross-correlation function normalized by the bias, at large scales, for different perturbations of the observed halo number count. This leads to: *upper left-hand panel* only the contribution from gravitational potential was taken into account as a source of RSD, in black dashed line we have the prediction when accounting for leading terms in $(H/k)^2$. *Upper right-hand panel* Doppler only, *middle left-hand panel* transverse Doppler only, *middle right-hand panel* ISW/RS only, *bottom left-hand panel* weak lensing only, and finally *bottom right-hand panel* the residual where we subtract all the previous effects to the full dipole taking into account all the effects at once. In black we have the averaged prediction using linear theory at first order in $H/k$.</figcaption></figure>

Last, there is a difference between the particle positions on the FLRW light cone given by the simulation and the particle position on the perturbed FLRW light cone seen by photons due to time delay (Shapiro effect). However this contribution is too small to be detectable (inferior to $10^{-5}$) as seen in Fig. 9 and will therefore be neglected.

## 5.2 The linear regime and its breakdown

From now on, we subtract the effects from statistical fluctuations as well as the light-cone effect described above. In this section we in-

vestigate the dipole at large scales from $20 h^{-1}$ Mpc to $150 h^{-1}$ Mpc. This corresponds to the linear regime (Section 2.2) and the beginning of the quasi-linear regime where, in principle, linear theory does not hold any more. We focus on the weighted average of the normalized dipole. Each cross-correlation is shown in Appendix C.

### 5.2.1 Contributions to the dipole at large scales

At first order in the metric perturbation, RSD are the sum of five contributions: four redshift perturbations (see equations 44–49) plus
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Full dipole of the cross-correlation function normalized by the bias. The dipole is dominated by the Doppler contribution.</figcaption></figure>

the lensing effect which affects the apparent position of haloes (we neglect the small Shapiro term which is already subtracted). To investigate each of these effects we have produced five catalogues where only one effect is present at a time. This allows us to compare our measurements to the linear predictions summed up in Table 1 as shown in Fig. 10.

The upper left-hand panel shows the dipole when we only consider the gravitational potential as a source of RSD. We find a good agreement at large scales with linear theory (even better with higher order terms in $H/k$) down to a radius $r \approx 60–80 h^{-1}$ Mpc where $\xi_1/\Delta b \approx 5 \times 10^{-5}$. The dipole drops sharply at smaller scales and even becomes negative at scales smaller than $30 h^{-1}$ Mpc. Linear theory fails to predict this drop. As we shall see, it is the second most important contribution to the total dipole at large scales. The upper right-hand panel shows the dipole for the Doppler term only. The signal is much larger rising from $\xi_1/\Delta b \approx 2 \times 10^{-4}$ at $100 h^{-1}$ Mpc up to $\xi_1/\Delta b \approx 10^{-3}$ at $20 h^{-1}$ Mpc. Data points and linear prediction are in agreement when looking at the scatter (except for two points at 65 and $75 h^{-1}$ Mpc which are at about $3-\sigma$ above the linear expectation). This is the dominant term and it is mostly related to the divergence of the line of sights at this redshift (this term is inversely proportional to the comoving distance and is therefore fainter at higher redshift). In the middle left-hand panel we only consider the transverse Doppler term which should be null in the linear regime. Our measurement is consistent with zero between $r \approx 60 h^{-1}$ Mpc and $r \approx 150 h^{-1}$ Mpc at a precision better than $2 \times 10^{-6}$. Below $60 h^{-1}$ Mpc the data show a positive dipole. The transverse Doppler yields an overall redshift of the lighter galaxies w.r.t. to the more massive ones which is explained (at smaller scales) by Zhao et al. (2013). In the middle right-hand panel we see the dipole of the ISW/RS effect only. This effect is integrated and suppressed by a factor $H/k$ compared to the other terms and is therefore expected to be small. Here we only check the consistency with zero that is given on the full scales of interest with a precision of order $\lesssim 5 \times 10^{-6}$. Bottom left shows the effect of weak lensing only. This effect is much fainter than the potential and Doppler terms and the data clearly favours a negative dipole of order of a few $10^{-6}$ which follows well the linear prediction.

Lastly, the bottom right-hand panel shows the *residual*, i.e. the full dipole including all effects (full redshift perturbations and lensing) minus the individual contributions mentioned above. It includes all the cross-terms (potential-Doppler, potential-lensing, etc...) as well as the non-linear mapping that was ignored by the linear calculation.

We see that the *residual* is consistent with zero to a good precision beyond $50 h^{-1}$ Mpc. However below this threshold the *residual* gives a negative contribution to the total dipole. We expect this term to be dominated by the potential-Doppler cross-correlation term as they are the two dominant terms. We can see that the departure of the *residual* from zero occurs at approximately the same scale as the departure of the dipole for the potential only term from the linear prediction. It indicates a failure of the linear regime at this scale. In Appendix C we show the behaviour of the *residual* for a wide range of cross-correlations. Depending on the values of bias and bias difference, the amplitude of the signal as well as the scale at which the *residual* departs from zero can change.

### 5.2.2 Total dipole

Now that we have seen all the individual contributions to the dipole, we show the final result that is the sum of all the previous contributions (see Fig. 11). As we expected the dipole is dominated by the Doppler term as the effects of the others terms are small. However the dipole departs from linear theory near $30 h^{-1}$ Mpc due to the non-linear contribution from the potential and residual. In Bonvin et al. (2014) the authors claim that a measurement of the total dipole in the linear regime will allow us to probe velocity field and to test general relativity (through the Euler equation). As we have seen to reach this goal, one has to pay attention to real-space statistical fluctuations of the dipole (therefore huge sample and survey volume are mandatory), evolution effect (as we need to properly model bias evolution while for the moment we are limited to simple phenomenological models) and the bias itself (as the linear prediction is proportional to the bias difference between the population and simple scale-independent bias models are considered).

## 5.3 From quasi-linear scales to non-linear scales

After investigating the linear regime (beyond $50 h^{-1}$ Mpc) and its breakdown (at $40–60 h^{-1}$ Mpc scales for the potential contribution and the *residual*), we now focus on the transition from quasi-linear to non-linear scales between $30 h^{-1}$ Mpc and $5 h^{-1}$ Mpc. We use the conservative lower bound of $5 h^{-1}$ Mpc because it stands well above the coarse grid size ($0.6 h^{-1}$ Mpc), beyond halo's virial radii (~$0.9 h^{-1}$ Mpc for large group-size haloes) and it is larger than the light-cone shell size (~$4 h^{-1}$ Mpc). In this regime, we expect baryonic effects to be negligible although this has to be further investigated with dedicated simulations.

The theoretical predictions are not necessarily proportional to the bias difference, thus showing the weighted average of the normalized dipole would not make much sense. Here we focus on the cross-correlation of $\text{data}_{H_{1600}}$ (halo mass $M_{h1} \approx 4.5 \times 10^{13} h^{-1} M_{\odot}$) with $\text{data}_{H_{100}}$ (halo mass $M_{h2} \approx 2.8 \times 10^{12} h^{-1} M_{\odot}$), which gives a bias difference $\Delta b \approx 1$.

### 5.3.1 Contributions to the dipole at small scales

The various contributions to the dipole are shown in Fig. 12. The dominant contributions at non-linear scales (especially at $r < 10 h^{-1}$ Mpc) are very different from the ones at linear scales. Here the dominant contributions are the potential term (upper left) and the *residual* (bottom right) while in the linear regime the dominant contribution is the Doppler term (upper right). Moreover both contributions are very negative resulting in a negative dipole.
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Dipole of the cross-correlation function between data_H1600 and data_H100, at small scales, for different perturbations of the halo number count. This leads to: *upper left-hand panel* only the contribution from gravitational potential was taken into account as a source of RSD, *upper right-hand panel* Doppler only, *middle left-hand panel* transverse Doppler only, *middle right-hand panel* ISW/RS only, *bottom left-hand panel* weak lensing only, and finally *bottom right-hand panel* the residual for which we subtract all the previous effects to the full dipole taking into account all the effects at once.</figcaption></figure>

The potential contribution (upper left) drops slowly from $\xi_1 \simeq 0$ near $30 \ h^{-1} \ Mpc$ to $\xi_1 \simeq -1 \times 10^{-4}$ at $20 \ h^{-1} \ Mpc$. The fall is then much steeper at smaller separations, from $\xi_1 \simeq -5 \times 10^{-4}$ at $14 \ h^{-1} \ Mpc$ down to $\xi_1 \simeq -6 \times 10^{-3}$ at $6 \ h^{-1} \ Mpc$. We also note that the measurement is very robust since the statistical error bars are very small (although error bars should be taken with caution at these scales). The linear prediction completely fails. The dipole of the halo-halo cross-correlation is a sensitive probe of the gravitational potential up to about ten virial radii for this halo mass.

The velocity contribution (upper right) remains high $\xi_1 \simeq 5-20 \times 10^{-4}$ between 30 and $6 \ h^{-1} \ Mpc$. At smaller scales the error bars increase from $\sigma_\xi \simeq 5 \times 10^{-4}$ to $\sigma_\xi \simeq 5 \times 10^{-3}$ at smaller

scales. Interestingly the Doppler-only dipole remains close to the linear expectation.

The transverse-Doppler contribution to the total redshift in the vicinity of galaxy clusters has originally been highlighted by Zhao et al. (2013). However it was restricted to the region $r < 2 \ R_{vir}$ inside or close to the virial radius $R_{vir} \sim 1-2 \ h^{-1} \ Mpc$ of the clusters. Interestingly, the transverse-Doppler contribution to the dipole (middle-left) is non-zero even at very large radii ($r > 2 \ R_{vir}$). It remains positive of order $\xi_1 \simeq 2-6 \times 10^{-5}$ at radii $14 < r < 30 \ h^{-1} \ Mpc$. At smaller scales there is strong increase from $\xi_1 = 2 \times 10^{-4}$ at $14 \ h^{-1} \ Mpc$ to $\xi_1 = 5 \times 10^{-4}$ at $6 \ h^{-1} \ Mpc$. The ratio to the potential contribution to the dipole is of order -10 at this scale.
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Full dipole of the cross-correlation function between data.H1600 and data.H100. The deviation from linear theory is governed by the potential contribution and the 'residual' (mostly related to the coupling between potential and velocity terms). The dipole is a sensitive probe of the potential well beyond the virial radius of haloes.</figcaption></figure>

The ISW contribution (middle right) and lensing contribution (bottom left) are consistent with zero at small scales. The size of the error bars provide an upper limit for the signal of $\xi_1 < 5 \times 10^{-5}$ for ISW and $\xi_1 < 10^{-4}$ for lensing. It is still in agreement with the linear prediction which is of the same order of magnitude, however the fluctuations are too important to measure the signal.

Surprisingly, the residual (bottom right) is of the same order as the potential contribution (from $\sim -10^{-4}$ at $30 h^{-1}$ Mpc to $\sim -6 \times 10^{-3}$ at $6 h^{-1}$ Mpc). This is an important result of this paper. It means that at these scales and especially below $15 h^{-1}$ Mpc, one cannot add up all the contributions one by one. On the contrary, there are some important contributions involving both potential terms and velocity terms together.

### 5.3.2 Total dipole

The total dipole at non-linear scales is presented in Fig. 13. It remains slightly positive of order $\xi_1 \sim 1 \times 10^{-3}$ above $15 h^{-1}$ Mpc. As shown in the previous section, this is related to the velocity contribution which remains positive in this region. At smaller scales, the potential contribution dominates over the velocity contribution. The total dipole is then falling down quickly to $\xi_1 \sim -1 \times 10^{-2}$ at $6 h^{-1}$ Mpc. Moreover within our simulated survey of 8.34 ($h^{-1}$ Gpc)$^3$, error bars (mostly related to the fluctuations of the velocity field) are smaller than the signal at this scale. The dipole of the group-galaxy cross-correlation function is therefore a good probe of the potential far outside of the group virial radii. Interestingly, deviations from linear theory are mostly governed by the potential and by the *residual*. The interpretation of the dipole is therefore non-trivial because of correlations between potential and velocity terms. However the dipole carries important information about the potential.

### 5.3.3 Mass dependence of the contributions

So far, we have focused on the cross-correlation between haloes of mass $\sim 4.5 \times 10^{13} h^{-1}$ M$_\odot$ and haloes of mass $\sim 2.8 \times 10^{12} h^{-1}$ M$_\odot$. In Fig. 14, we investigate the halo mass dependence of the main dipole contributions (velocity, potential). The mass dependence on the *residual* is shown in Appendix C. We

explore various configurations by cross-correlating all the different halo populations with the lightest halo population. At large linear scales the variation of the dipole is mostly governed by the bias difference between the two halo populations, however at small non-linear scales the evolution of the dipole is less trivial. The velocity contribution to the dipole does not evolve strongly with halo mass. It stays bounded in the range $0 < \xi_1 < 1 \times 10^{-3}$. On the other hand, the potential contribution becomes more negative at larger mass from $\xi_1 \simeq -5 \times 10^{-4}$ to $\xi_1 \simeq -1 \times 10^{-2}$ at $6 h^{-1}$ Mpc. It means that for massive enough haloes the potential contribution dominates over the velocity contribution for a wide range of scales (as seen previously). However for haloes lighter than $\sim 10^{13} h^{-1}$ M$_\odot$ the velocity-contribution dominates. The *residual* also departs from 0 at larger radii for heavier haloes. Interestingly it is mostly following the potential contribution.

The prediction of the potential effect from equation (41) (assuming spherical symmetry) reproduces the trend at a qualitative level. However the potential contribution is overestimated. Taking into account the dispersion around the potential deduced from spherical symmetry as in equation (38) should improve the agreement with the measured dipole (Cai et al. 2017). Note that we have checked (see Appendix B) that our conclusions still hold for a very different halo definition (i.e. linking length $b = 0.1$). The main difference is a slightly better agreement with the spherical predictions for the potential contribution to the dipole.

## 6 CONCLUSIONS

In this work we explored the galaxy clustering asymmetry by looking at the dipole of the cross-correlation function between halo populations of different masses (from Milky Way size to galaxy-cluster size). We took into account all the relevant effects which contribute to the dipole, from lensing to multiple redshift perturbation terms. At large scales we obtain a good agreement between linear theory and our results. At these scales the dipole can be used as a probe of velocity field (and as a probe of gravity through the Euler equation). However one has to consider a large enough survey to overcome important real-space statistical fluctuations. It is also important to take into account the light-cone effect and to accurately model the bias and its evolution.

At smaller scales we have seen deviation from linear theory. Moreover the gravitational redshift effect dominates the dipole below $10 h^{-1}$ Mpc. It is therefore possible to probe the potential outside groups and clusters using the dipole. By subtracting the linear expectation for the Doppler contribution it is in principle possible to probe the potential to even larger radii. This is a path to explore in order to circumvent the disadvantages of standard probes of the potential, usually relying on strong assumptions (such as hydrostatic equilibrium) or being only sensitive to the projected potential (lensing). A simple spherical prediction allows to predict the global trend of the dipole but not the exact value. Moreover as we have seen the *residual* (i.e all the cross terms and non-linearities of the mapping) is of the same order as the gravitational potential contribution and should be taken into account properly. At small scales the pairwise velocity PDF is also highly non-Gaussian, leading to high peculiar velocities and Finger-of-God effect. Coupled to gravitational potential and possibly wide-angle effect we expect this to be a non-negligible contribution to the dipole. To fully understand and probe cosmology or modified theories of gravity at these scales using the cross-correlation dipole we therefore need a perturbation theory or streaming model which takes into account more redshift
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Dipole of the cross-correlation function between different data sets and data_H100 (no autocorrelation). **Left-hand panels:** gravitational potential only, dotted lines gives the spherical prediction computed using equation (41). **Middle panels:** Doppler only. For massive enough halo the negative potential contribution dominates over the positive Doppler contribution. The linear prediction is given by dashed lines. **Right-hand panels:** residual term.</figcaption></figure>

perturbation terms and relaxes the distant observer approximation. This will be the focus of a future paper.

There are multiple possible extensions to this work. At large Gpc scales current analysis are limited by the volume of the simulation as well as gauge effect. At smaller scales the baryons as well as the finite resolution effect might play a role. Extension of this work

in these two directions can open interesting perspectives. When analysing future surveys, it is also important to consider observational effects. One possibility would be to populate haloes with galaxies and to incorporate effects such as magnification bias, absorption by dust, redshift errors, alignment of galaxies, etc. Another straight-forward extension is to explore the influence of cosmology,
dark energy, dark matter and modified gravity on the dipole of the
halo cross-correlation to shed light on the nature of the dark sector
with future large scale surveys.

ACKNOWLEDGEMENTS

MAB and YR acknowledge financial support from the DIM ACAV of the Region Ile-de-France, the Action fédératrice *Cosmologie et structuration de l'univers* as well as the JSPS Grant L16519. AT acknowledges financial support from MEXT/JSPS KAKENHI Grant Numbers JP15H05889 and JP16H03977. SS acknowledge financial support from Grant-in-Aid for JSPS Research Fellow Number 17J10553.

This work was granted access to HPC resources of TGCC/CINES through allocations made by GENCI (Grand Equipement National de Calcul Intensif) under the allocations 2016-042287, 2017-A0010402287 and 2018-A0030402287. We thank Fabrice Roy for important technical support and providing pFoF, and Vincent Reverdy for providing MAGRATHEA. We also thank C. Murray for comments on the catalogues, R. Teyssier for RAMSES, D. Alonso for CUTE and B. Li for sharing his TSC routine. MAB thanks Pier-Stefano Corasaniti, Stefano Anselmi and Paul de Fromont for fruitful discussions. YR thanks Ichihiko Hashimoto and Julian Adamek.

To complete this work, discussions during the workshop, YITP-T-17-03, held at Yukawa Institute for Theoretical Physics (YITP) at Kyoto University were useful. Numerical calculations for the present work have been carried out partly at the Yukawa Institute Computer Facility.

REFERENCES

Adamek J., 2018, Phys. Rev., 97, 021302

Adamek J., Daverio D., Durrer R., Kunz M., 2016a, J. Cosmol. Astropart. Phys., 7, 053

Adamek J., Daverio D., Durrer R., Kunz M., 2016b, Nature Phys., 12, 346

Alam S., Zhu H., Croft R. A. C., Ho S., Giusarma E., Schneider D. P., 2017, MNRAS, 470, 2822

Alonso D., 2012, preprint (arXiv:1210.1833)

Balmès I., Rasera Y., Corasaniti P.-S., Alimi J.-M., 2014, MNRAS, 437, 2328

Bardeen J. M., 1980, Phys. Rev., 22, 1882

Barreira A., Llinares C., Bose S., Li B., 2016, J. Cosmol. Astropart. Phys., 5, 001

Bartelmann M., Schneider P., 2001, Phys. Rep., 340, 291

Ben-Dayan I., Marozzi G., Nugier F., Veneziano G., 2012, J. Cosmol. Astropart. Phys., 11, 045

Bonvin C., Durrer R., 2011, Phys. Rev., 84, 063505

Bonvin C., Hui L., Gaztañaga E., 2014, Phys. Rev., 89, 083535

Bonvin C., Hui L., Gaztañaga E., 2016, J. Cosmol. Astropart. Phys., 8, 021

Borzyszkowski M., Bertacca D., Porciani C., 2017, MNRAS, 471, 3899

Cai Y.-C., Kaiser N., Cole S., Frenk C., 2017, MNRAS, 468, 1981

Cappi A., 1995, A&A, 301, 6

Carlson J., White M., Padmanabhan N., 2009, Phys. Rev., 80, 043531

Carlson J., Reid B., White M., 2013, MNRAS, 429, 1674

Carrasco J. J. M., Hertzberg M. P., Senatore L., 2012, J. High Energy Phys., 9, 82

Challinor A., Lewis A., 2011, Phys. Rev., 84, 043516

Chisari N. E., Zaldarriaga M., 2011, Phys. Rev., 83, 123505

Crocce M., Scoccimarro R., 2006, Phys. Rev., 73, 063519

Crocce M., Scoccimarro R., 2008, Phys. Rev., 77, 023533

Croft R. A. C., 2013, MNRAS, 434, 3008

Fidler C., Rampf C., Tram T., Crittenden R., Koyama K., Wands D., 2015, Phys. Rev., 92, 123517

Fidler C., Tram T., Rampf C., Crittenden R., Koyama K., Wands D., 2016, J. Cosmol. Astropart. Phys., 9, 031

Fluke C. J., Lasky P. D., 2011, MNRAS, 416, 1616

Fluke C. J., Webster R. L., Mortlock D. J., 1999, MNRAS, 306, 567

Fosalba P., Gaztañaga E., Castander F. J., Manera M., 2008, MNRAS, 391, 435

Fosalba P., Gaztañaga E., Castander F. J., Crocce M., 2015a, MNRAS, 447, 1319

Fosalba P., Crocce M., Gaztañaga E., Castander F. J., 2015b, MNRAS, 448, 2987

Gaztanaga E., Bonvin C., Hui L., 2017, J. Cosmol. Astropart. Phys., 1, 032

Giblin J. T., Jr, Mertens J. B., Starkman G. D., Zentner A. R., 2017, Phys. Rev., 96, 103530

Giusarma E., Alam S., Zhu H., Croft R. A. C., Ho S., 2017, preprint (arXiv: 1709.07854)

Hamilton A. J. S., 1992, ApJ, 385, L5

Heitmann K. et al., 2016, ApJ, 820, 108

Hilbert S., Hartlap J., White S. D. M., Schneider P., 2009, A&A, 499, 31

Hui L., Gaztañaga E., Loverde M., 2007, Phys. Rev., 76, 103502

Hui L., Gaztañaga E., Loverde M., 2008, Phys. Rev., 77, 063526

Kaiser N., 1987, MNRAS, 227, 1

Kaiser N., 2013, MNRAS, 435, 1278

Kerscher M., Szapudi I., Szalay A. S., 2000, ApJ, 535, L13

Kibble T. W. B., Lieu R., 2005, ApJ, 632, 718

Killedar M., Lasky P. D., Lewis G. F., Fluke C. J., 2012, MNRAS, 420, 155

Komatsu E. et al., 2011, ApJS, 192, 18

Landy S. D., Szalay A. S., 1993, ApJ, 412, 64

Lewis A., Challinor A., Lasenby A., 2000, ApJ, 538, 473

Li B., Zhao G.-B., Teyssier R., Koyama K., 2012, J. Cosmol. Astropart. Phys., 1, 051

Ma C.-P., Bertschinger E., 1995, ApJ, 455, 7

Matsubara T., 2000, ApJ, 537, L77

Matsubara T., 2008, Phys. Rev. D, 77, 063530

McDonald P., 2009, J. Cosmol. Astropart. Phys., 11, 026

Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493

Norberg P., Baugh C. M., Gaztañaga E., Croton D. J., 2009, MNRAS, 396, 19

Pápai P., Szapudi I., 2008, MNRAS, 389, 292

Prunet S., Pichon C., Aubert D., Pogosyan D., Teyssier R., Gottloeber S., 2008, ApJS, 178, 179

Raccanelli A., Bertacca D., Jeong D., Neyrinck M. C., Szalay A. S., 2018, Phys. Dark Universe, 19, 109

Rasera Y., Corasaniti P.-S., Alimi J.-M., Bouillot V., Reverdy V., Balmès I., 2014, MNRAS, 440, 1420

Reid B. A., White M., 2011, MNRAS, 417, 1913

Reimberg P., Bernardeau F., Pitrou C., 2016, J. Cosmol. Astropart. Phys., 1, 048

Reverdy V., 2014, PhD thesis, Laboratoire Univers et Théories

Roy F., Bouillot V. R., Rasera Y., 2014, A&A, 564, A13

Schneider P., Ehlers J., Falco E. E., 1992, Gravitational Lenses, Springer, Berlin Heidelberg

Scoccimarro R., 2004, Phys. Rev., 70, 083007

Sheth R. K., Tormen G., 1999, MNRAS, 308, 119

Szapudi I., 2004, ApJ, 614, 51

Tansella V., Bonvin C., Durrer R., Ghosh B., Sellentin E., 2018, J. Cosmol. Astropart. Phys., 03, 019

Taruya A., Nishimichi T., Saito S., 2010, Phys. Rev., 82, 063522

Taruya A., Bernardeau F., Nishimichi T., Codis S., 2012, Phys. Rev., 86, 103528

Taruya A., Nishimichi T., Bernardeau F., 2013, Phys. Rev., 87, 083509

Teyssier R. et al., 2009, A&A, 497, 335

Teyssier R., 2002, A&A, 385, 337

Tinker J., Kravtsov A. V., Klypin A., Abazajian K., Warren M., Yepes G., Gottlöber S., Holz D. E., 2008, ApJ, 688, 709

Tinker J. L., 2007, MNRAS, 374, 477

Tinker J. L., Weinberg D. H., Zheng Z., 2006, MNRAS, 368, 85

Umeh O., Clarkson C., Maartens R., 2014, Class. Quantum Gravity, 31, 202001
Wojtak R., Hansen S. H., Hjorth J., 2011, *Nature*, **477**, 567

Yoo J., 2010, *Phys. Rev.*, **82**, 083508

Yoo J., Fitzpatrick A. L., Zaldarriaga M., 2009, *Phys. Rev.*, **80**, 083514

Zaroubi S., Hoffman Y., 1996, *ApJ*, **462**, 25

Zhao H., Peacock J. A., Li B., 2013, *Phys. Rev.*, **88**, 043013

Zhao D. H., Jing Y. P., Mo H. J., Börner G., 2009, *ApJ*, **707**, 354

Zhu H., Alam S., Croft R. A. C., Ho S., Giusarma E., 2017, *MNRAS*, **471**, 2345

## APPENDIX A: BIAS MEASUREMENT AND EVOLUTION

In this section we present the methodology used to compute the bias for each halo population summarized in Table 2. In our case we will consider a scale-independent bias which can be directly fitted with a constant using equation (53). Even if we perform our analysis on the light cone, computing the bias from the full light-cone monopole would give very poor results (see Fig. 4) for haloes with increasing mass. We therefore compute the bias using snapshots with the same comoving volume as our full light cone, and interpolating between snapshots to find the bias of a given population at a given redshift. From a numerical point of view, it would be more accurate to use the halo-matter cross-correlation instead of square root of the halo-halo autocorrelation since we would have more particles than

<figure><img src="image_17.png" /><figcaption>**Figure A1.** Top panel: bias of the halo populations summarized in Table 2 computed on snapshots. Bottom panel: Derivative of bias w.r.t. comoving distance, comparison of our computation (circles, dashed lines) with the prediction from the Sheth & Tormen (1999) model equation (36) (full lines). The associated halo populations are, from top to bottom: data_H3200, data_H1600, data_H800, data_H400, data_H200, data_H100.</figcaption></figure>

haloes and therefore less noise. However this would mean saving all the particles for each snapshot which is unrealistic for our simulation. Instead, we compute $\xi_{mm}$ using COSMICEMU (Heitmann et al. 2016) which is in good agreement which the matter autocorrelation monopole we have on the light cone (see Fig. 3). We also checked that the emulator power spectrum is in agreement with the power spectrum we compute in each snapshot at one per cent level between roughly $k = 0.02$ and $2h$ Mpc$^{-1}$. Finally, we cross-correlate each halo population with the least massive one in order to increase the statistics instead of computing the halo autocorrelation. For the data_H100 bias we indeed use

$$b_{100} \approx \sqrt{\frac{\xi_{hh}^{\ell=0}}{\xi_{mm}^{\ell=0}}}. \qquad (A1)$$

For other populations we then take

$$b_1 = b_{100} \frac{\xi_{h_1 h_2}^{\ell=0}}{\xi_{mm}^{\ell=0}}, \qquad (A2)$$

where $h_1$ denotes the halo population from which we want to know the bias, $h_2$ denotes the data_H100 population with $b_{100}$ its bias. This was done on 12 snapshots which covers well our full-sky light cone. The associated redshift are: $z = 0.123, 0.152, 0.180, 0.208, 0.236, 0.250, 0.265, 0.296, 0.329, 0.364, 0.399$, and 0.428. The results are shown in Fig. A1.

For bias derivative shown in Fig. A1, we use the prediction from Sheth & Tormen (1999), which gives us as a result equation (36). In our case, for each mass bin (increasing in mass) we set $\sigma(M, 0) \equiv \sigma_0 = 1.60, 1.44, 1.27, 1.17, 0.95$, and 0.9.

For the most massive bins the Sheth & Tormen (1999) computation of the bias can be very different from our numerical results, therefore for the dipole predictions we keep the bias from our computation on snapshots but we take the bias derivative from equation (36) as it is smoother.

## APPENDIX B: DIFFERENT HALO DEFINITION

In this appendix we perform the same analysis as in core of the paper but with a different halo definition. In the paper we detected haloes using the friend-of-friend algorithm with $b = 0.2$. To check the sensitivity of our results about the halo definition we here consider a very different linking length, namely $b = 0.1$. While $b = 0.2$ corresponds to an enclosed overdensity very roughly of order of 200 times the mean density of the universe, $b = 0.1$ corresponds to an enclosed overdensity much larger (approximately eight times more). Most of the usual halo definitions lie somewhere in between these two definitions. The data sets are shown in Table. B1, where the bias is computed on the full light cone contrarily to previous

<table><caption>**Table B1.** Summary of the different data sets used: name, number of haloes, range for the number of particles per halo, mean mass, bias at the volume averaged redshift $z = 0.341$.</caption><thead><tr><th>Name</th><th>Number of haloes</th><th>Number of part</th><th>Mass (h<sup>-1</sup> M<sub>☉</sub>)</th><th>Bias</th></tr></thead><tbody><tr><td>data2_H050</td><td>6.9 &times; 10<sup>6</sup></td><td>50–100</td><td>1.4 &times; 10<sup>12</sup></td><td>1.18</td></tr><tr><td>data2_H100</td><td>3.7 &times; 10<sup>6</sup></td><td>100–200</td><td>2.8 &times; 10<sup>12</sup></td><td>1.38</td></tr><tr><td>data2_H200</td><td>2.0 &times; 10<sup>6</sup></td><td>200–400</td><td>5.6 &times; 10<sup>12</sup></td><td>1.54</td></tr><tr><td>data2_H400</td><td>1.0 &times; 10<sup>6</sup></td><td>400–800</td><td>1.1 &times; 10<sup>13</sup></td><td>1.76</td></tr><tr><td>data2_H800</td><td>4.3 &times; 10<sup>5</sup></td><td>800–1600</td><td>2.2 &times; 10<sup>13</sup></td><td>2.13</td></tr><tr><td>data2_H1600</td><td>1.7 &times; 10<sup>5</sup></td><td>1600–3200</td><td>4.5 &times; 10<sup>13</sup></td><td>2.60</td></tr></tbody></table>
<figure><img src="image_18.png" /></figure>

**Figure B1.** Dipole of the cross-correlation function normalized by the bias, at large scales, for different perturbations of the observed halo number count for another halo definition (b = 0.1). This leads to: *upper left-hand panel* only the contribution from gravitational potential was taken into account as a source of RSD, in black dashed line we have the prediction when accounting for leading terms in $(H/k)^2$. *Upper right-hand panel* Doppler only, *middle left-hand panel* transverse Doppler only, *middle right-hand panel* ISW/RS only, *bottom left-hand panel* weak lensing only, and finally *bottom right-hand panel* the residual where we subtract all the previous effects to the full dipole taking into account all the effects at once. In black solid lines we show the averaged prediction using linear theory at first order in $H/k$.
<figure><img src="image_19.png" /><figcaption>**Figure B2.** Full dipole of the cross-correlation function normalized by the bias for another halo definition ($b = 0.1$). The dipole is dominated by the Doppler contribution.</figcaption></figure>

data sets where the bias was estimated by interpolating between snapshots. For previous data sets these two methods agree at the 1 per cent level.

For linear scales in Fig. B1, we see that on a qualitative level the results are similar to the ones using another halo definition (Fig. 10). Quantitatively, it seems that we slightly underestimate the Doppler effect. The remarks on the full dipole shown in Fig. B2 are also similar. More interestingly, for the quasi-linear and non-linear regime in Fig. B3 we can notice three things. First the results are qualitatively similar to the haloes with $b = 0.2$. Second there is a better agreement with $b = 0.1$ with the spherical prediction for the non-linear dipole for the potential only term. This is due to the fact that haloes are now more clustered and sit in deeper potential wells, enhancing the amplitude of the dipole. Still we see that it does not match completely the theoretical prediction. Last, we see that similarly to $b = 0.2$, the point at 6 h⁻¹ Mpc for Doppler only is very negative for the last correlation (most massive halo population with the lightest). This may be due to a coupling between Finger-of-God and wide-angle effects.
<figure><img src="image_20.png" /><figcaption>**Figure B3.** Dipole of the cross-correlation function between different data sets and data2\_H₀₅₀ (i.e. for another halo definition *b* = 0.1). Left-hand panels: gravitational potential only, dotted lines gives the spherical prediction computed using equation (41). Middle panels: Doppler only. For massive enough halo the negative potential contribution dominates over the positive Doppler contribution. The linear prediction is given by dashed lines. Right-hand panels: residual term. In the bottom plot of the middle column, the point at 6 *h*⁻¹ Mpc is at ξ₁ = −0.03.</figcaption></figure>
APPENDIX C: MASS DEPENDENCE OF THE DIPOLE

In Section 5.1 and 5.2, we presented the computation of the dipole normalised by the bias difference. To do so we used all the cross-correlations available with our data sets shown in Table 2. We then performed a sum on the dipoles, weighted by the inverse of their

variance (see Section 3.3). In this Section we show the different cross-correlations for each perturbation effect (potential only in Fig. C1, Doppler only in Fig. C2, transverse Doppler in Fig. C3, residual in Fig. C4 and the full dipole in Fig. C5), and for every combination of populations at large scales. We show the results for the computation of the cross-correlation on the full light cone using jackknife re-sampling.

<figure><img src="image_21.png" /><figcaption>**Figure C1.** Potential only term of the dipole of the cross-correlation function on the full light cone at large scales. The linear predictions at first order in $\mathcal{H}/k$ are shown in dash-dotted lines while the prediction with the dominant $(\mathcal{H}/k)^2$ terms is shown in dashed lines.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>**Figure C2.** Doppler only term of the dipole of the cross-correlation function on the full light cone at large scales. The linear predictions are shown in dashed lines.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>**Figure C3.** Transverse Doppler only term of the dipole of the cross-correlation function on the full light cone at large scales.</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>**Figure C4.** Residual term of the dipole of the cross-correlation function on the full light cone at large scales.</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>**Figure C5.** Full dipole of the cross-correlation function on the full light cone at large scales. The linear predictions are shown in dashed lines.</figcaption></figure>

This paper has been typeset from a TeX/LA TeX file prepared by the author.