# The RayGalGroupSims cosmological simulation suite for the study of relativistic effects: an application to lensing-matter clustering statistics

Y. Rasera, M-A. Breton, P-S. Corasaniti, J. Allingham, F. Roy, V. Reverdy, T. Pellegrin, S. Saga, A. Taruya, S. Agarwal, et al.

► To cite this version:

Y. Rasera, M-A. Breton, P-S. Corasaniti, J. Allingham, F. Roy, et al. The RayGalGroupSims cosmological simulation suite for the study of relativistic effects: an application to lensing-matter clustering statistics. Astronomy and Astrophysics - A&A, 2022, 661, pp.A90. 10.1051/0004-6361/202141908. hal-03463461

HAL Id: hal-03463461

https://hal.science/hal-03463461

Submitted on 10 Feb 2023

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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
# The RayGalGroupSims cosmological simulation suite for the study of relativistic effects: An application to lensing-matter clustering statistics

Y. Rasera<fnref n="1" />, M.-A. Breton<fnref n="2" />, P.-S. Corasaniti<fnref n="3" />,<fnref n="4" />, J. Allingham<fnref n="5" />, F. Roy<fnref n="3" />, V. Reverdy<fnref n="6" />,<fnref n="3" />,<fnref n="7" />, T. Pellegrin<fnref n="8" />,<fnref n="9" />,<fnref n="3" />, S. Saga<fnref n="3" />, A. Taruya<fnref n="10" />, S. Agarwal<fnref n="11" />, and S. Anselmi<fnref n="12" />,<fnref n="3" />

Received 29 July 2021 / Accepted 27 December 2021

## ABSTRACT

**Context.** General relativistic effects on the clustering of matter in the Universe provide a sensitive probe of cosmology and gravity theories that can be tested with the upcoming generation of galaxy surveys. These will require the availability of accurate model predictions, from large linear scales to small non-linear ones.

**Aims.** Here, we present a suite of large-volume high-resolution *N*-body simulations specifically designed to generate light-cone data for the study of relativistic effects on lensing-matter observables without the use of simplifying approximations. As a case study application of these data, we perform an analysis of the relativistic contributions to the lensing-matter power spectra and cross-power spectra.

**Methods.** The RayGalGroupSims suite (RAYGAL for short) consists of two *N*-body simulations of (2625 $h^{-1}$ Mpc)$^3$ volume with 4096$^3$ particles of a standard flat $\Lambda$CDM model and a non-standard $w$CDM phantom dark energy model with a constant equation of state. Light-cone data from the simulations have been generated using a parallel ray-tracing algorithm that has integrated more than 1 billion geodesic equations without the use of the flat-sky or Born approximation.

**Results.** Catalogues and maps with relativistic weak lensing that include post-Born effects, magnification bias (MB), and redshift-space distortions (RSDs) due to gravitational redshift, Doppler, transverse Doppler, and integrated Sachs-Wolfe-Rees-Sciama effects are publicly released. Using this dataset, we are able to reproduce the linear and quasi-linear predictions from the CLASS relativistic code for the ten power spectra and cross-spectra (3 × 2 points) of the matter-density fluctuation field and the gravitational convergence at $z = 0.7$ and $z = 1.8$. We find a 1–30% level contribution from both MB and RSDs to the matter power spectrum, while the fingers-of-God effect is visible at lower redshift in the non-linear regime. Magnification bias also contributes at the 10–30% level to the convergence power spectrum, leading to a deviation between the shear power spectrum and the convergence power spectrum. Magnification bias also plays a significant role in the galaxy-galaxy lensing by decreasing the density-convergence spectra by 20% and coupling non-trivial configurations (such as the configuration with the convergence at the same redshift as the density, or at even lower redshifts).

**Conclusions.** The cosmological analysis shows that the relativistic 3×2 points approach is a powerful cosmological probe. Our unified approach to relativistic effects is an ideal framework for the investigation of gravitational effects in galaxy studies (e.g., clustering and weak lensing) as well as their effects in galaxy cluster, group, and void studies (e.g., gravitational redshifts and weak lensing) and cosmic microwave background studies (e.g., integrated Sachs-Wolfe-Rees-Sciama and weak lensing).

**Key words.** cosmology: theory – large-scale structure of Universe – dark energy – dark matter – gravitational lensing: weak – methods: numerical

## 1. Introduction

Observations of the distribution of galaxies in the Universe carry a wealth of cosmological information pertaining to the cosmic

matter content, the state of expansion, the amplitude of matter-density fluctuations, and the nature of dark energy and dark matter. They also enable general relativity (GR) to be tested on large cosmic scales. This is because, in addition to the cosmological

<footnote n="1">¹ Laboratoire Univers et Théories, Université de Paris, Observatoire de Paris, Université PSL, CNRS, 92190 Meudon, France
e-mail: yann.rasera@obspm.fr</footnote>

<footnote n="2">² Aix Marseille Univ, CNRS, CNES, LAM, Marseille, France
e-mail: michel-andres.breton@lam.fr</footnote>

<footnote n="3">³ Laboratoire Univers et Théories, Observatoire de Paris, Université PSL, Université de Paris, CNRS, 92190 Meudon, France</footnote>

<footnote n="4">⁴ Sorbonne Université, CNRS, UMR 7095, Institut d'Astrophysique de Paris, 98 bis bd Arago, 75014 Paris, France</footnote>

<footnote n="5">⁵ Sydney Institute for Astronomy, School of Physics, A28, The University of Sydney, NSW 2006, Australia</footnote>

<footnote n="6">⁶ Department of Mathematics (DMA), École Normale Supérieure, 45 rue d'Ulm, 75005 Paris, France</footnote>

<footnote n="7">⁷ National Center for Supercomputing Applications, University of Illinois at Urbana-Champaign, Illinois, USA</footnote>

<footnote n="8">⁸ LESIA, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université, Université de Paris, Meudon, France</footnote>

<footnote n="9">⁹ Laboratoire de Physique des Plasmas (LPP), École Polytechnique, IP Paris, Sorbonne Université, CNRS, Observatoire de Paris, Université PSL, Université Paris Saclay, Paris, France</footnote>

<footnote n="10">¹⁰ Center for Gravitational Physics, Yukawa Institute for Theoretical Physics, Kyoto University, Kyoto 606-8502, Japan</footnote>

<footnote n="11">¹¹ African Institute for Mathematical Sciences, 6 Melrose Road, Muizenberg 7945, Cape Town, South Africa</footnote>

<footnote n="12">¹² INFN, Sezione di Padova, via Marzolo 8, 35131 Padova, Italy</footnote>
imprint arising from the dynamical processes that lead to the formation of cosmic structures, observables of the clustering of galaxies also carry the signature of relativistic effects that alter the path of photons propagating between emitting galaxies and the observer.

Relativistic effects modify the apparent distribution of galaxies. They leave distinct imprints on galaxy and lensing two-points statistics through two processes: (i) the weak-lensing effect, which alters the apparent angular position, size, and shape of galaxies in the sky; and (ii) the redshift-space distortions (RSDs) due to the proper motion of the emitting galaxies, the local gravitational potentials, and the propagation of light rays through time-varying potentials.

The current as well as upcoming generation of galaxy surveys, such as eBOSS, HSC, DES, KIDS, DESI, LSST, Euclid, and SKA, will provide a precise mapping of the distribution of galaxies across a wide range of scales with measurements that are sensitive to such effects. These observational programmes will increase the number of galaxies with measured redshifts and/or ellipticities by several orders of magnitude, thus enabling accurate analyses of weak-lensing and clustering data to be performed (as well as their cross-spectra, the so-called 3 × 2 points analysis). This will allow strong constraints on the cosmological model parameters to be inferred and alternative models of dark energy and modified gravity theories to be tested, provided accurate theoretical model predictions are available.

The effects of GR on the galaxy two-points statistics beyond the RSD signal caused by the peculiar velocity of galaxies have been investigated in the linear regime in numerous works (see e.g., Yoo et al. 2009; Yoo 2010; Bonvin & Durrer 2011; Challinor & Lewis 2011), and extensions to the quasi-linear and non-linear regime have been primarily carried out using approximate analytical approaches (Di Dio & Seljak 2019; Taruya et al. 2020; Saga et al. 2020; Di Dio & Beutler 2020; Beutler & Di Dio 2020). Similarly, studies of GR effects on the weak-lensing shear have been investigated under various approximations. As an example, GR corrections up to second order in the gravitational potential of the lensing shear were computed in Bernardeau et al. (2010), and the imprint on the shear power spectrum in the linear regime was evaluated under the Born approximation in Di Dio et al. (2013).

Analytical and numerical studies of the GR contributions to the gravitational potentials of matter-density fluctuations have found that they leave a small imprint on the clustering of matter (Chisari & Zaldarriaga 2011; Adamek et al. 2016a). That is why Newtonian N-body simulations are still a valid tool for investigating relativistic kinematics effects (Fidler et al. 2015, 2016), which, on the other hand, leave observational imprints that cannot be neglected (see e.g., Ghosh et al. 2018; Breton et al. 2019).

Several general relativistic N-body and/or ray-tracing codes have been developed to study these effects in the non-linear regime of the clustering of matter (Adamek et al. 2016a; Borzyszkowski et al. 2017; Giblin et al. 2017; Sgier et al. 2021). Nevertheless, the use of approximations at different stages of the numerical computation prevents the possibility of inferring cosmological model predictions that are of general validity. For instance, a standard approach to generating light-cone data from high-resolution N-body simulations consists in assuming the Born approximation (Borzyszkowski et al. 2017) while possibly building the light cone from snapshots with the replica method (Sgier et al. 2021). A ray-tracing method based on the calculation of the lensing distortions of a ray bundle via the solving of the geodesic equations was presented in Killedar et al. (2012). However, the ray-tracing algorithm is limited to a fixed grid

and thus does not enable the high spatial resolution that can be achieved by N-body codes with adaptive mesh refinement (AMR) to be exploited.

Giblin et al. (2017) performed fully non-linear general relativistic simulations for an ideal fluid; however, the resolution remains limited. Adamek et al. (2016b) developed a numerical GR N-body code in the weak-field limit without AMR. Adamek et al. (2019) and Lepori et al. (2020) implemented a ray-tracing algorithm similar to that of Reverdy (2014) and Breton et al. (2019), which solves the geodesic equations without relying on the Born approximation; it was used to investigate several relativistic effects on the weak-lensing convergence power spectrum (Lepori et al. 2020), the matter power spectrum multipoles (Guandalin et al. 2021), and the cross-correlation power spectrum between the matter-density and lensing convergence (Coates et al. 2021).

Here, we describe and publicly release the numerical data generated from the RayGalGroupSims (which stands for Ray-tracing Galaxy Group Simulations, hereafter RAYGAL) simulation suite<fnref n="1" />, a set of simulations specifically designed to study the relativistic effects on the clustering of matter observables. In particular, the light-cone data have been generated using a ray-tracing algorithm that solves the geodesic equations with the gravitational potential evaluated along the different AMR levels of the N-body solver (Reverdy 2014; Breton et al. 2019; Breton & Fleury 2021). Thus, these numerical datasets benefit from a large dynamical range of simulations with a large volume and a high spatial resolution without the need to rely on the Born approximation or other limiting assumptions. In Breton et al. (2019), the light-cone data from the RAYGAL simulation of a Λ-cold dark matter model (hereafter ΛCDM) were used to perform a study of the asymmetry of the halo cross-correlation function due to the relativistic effects.

In this work, we introduce the complete simulation suite, including data from the simulation of a non-standard dark energy model, and present a case study application of the dataset. In particular, we present the evaluation of relativistic effects on the angular matter-density power spectrum, lensing convergence power spectrum, and their cross-spectra. We find a good agreement in the linear and quasi-linear regimes with the analytical predictions from the relativistic code CLASS, while deviations occur in the non-linear regime. Furthermore, we find a non-negligible imprint of RSDs and magnification bias (MB) on the harmonic matter power spectrum (1–30%). At low redshift, the damping due to the fingers-of-God effect is clearly visible. The MB contributes to the lensing convergence power spectrum to the 10–30% level. As the shear power spectrum is largely insensitive to the MB, this results in important differences with respect to the convergence power. Magnification bias also decreases the matter overdensity–gravitational convergence spectra (i.e. galaxy-galaxy lensing) by 20%. By coupling the signal between different redshift shells, relativistic effects boost the correlations between lensing-matter clustering observables. This is the case for the correlation between the matter-density field and the convergence evaluated at the same redshift as well as the correlation between the matter-density field and the convergence measured at a redshift smaller than that of the density field. The cosmological evolution shows that the relativistic contributions to the 3×2 points statistics make it promising probe of dark energy and modified gravity.

The paper is organised as follows: We describe the RAYGAL simulations and the associated numerical datasets in Sect. 2;

<footnote n="1">¹ https://cosmo.obspm.fr/public-datasets/</footnote>
<table><caption>**Table 1.** Cosmological parameter values of the ΛCDM and wCDM simulated models, respectively.</caption><thead><tr><th>Model</th><th>Ω<sub>m</sub></th><th>σ<sub>8</sub></th><th>w</th></tr></thead><tbody><tr><td>ΛCDM</td><td>0.25733</td><td>0.80101</td><td>-1.0</td></tr><tr><td>wCDM</td><td>0.27508</td><td>0.85205</td><td>-1.2</td></tr></tbody></table>

**Notes.** The columns from left to right give the values of the cosmic matter density, $\Omega_m$, the linear root-mean-square fluctuations smoothed on the 8 h$^{-1}$ Mpc scale, $\sigma_8$, and the redshift-independent equation of state, $w$. In all models the baryon density is set to $\Omega_b = 0.04356$, the density of the relativistic species to $\Omega_r \approx 0.00008$, the curvature to $\Omega_k = 0$, the reduced Hubble constant to $h = 0.72$, and the slope of the primordial spectrum to $n_s = 0.963$.

we present a case study of the imprint of relativistic effects on weak-lensing observables using the RAYGAL data in Sect. 3; and finally we discuss the conclusions in Sect. 4.

## 2. Methods: Simulations, ray tracing, and data

### 2.1. Cosmological models

The nature of dark energy still remains unknown. In order to test the cosmological dependence of relativistic effects on the clustering of matter observables, we simulated two different spatially flat cosmological models: a standard ΛCDM model and a ‘phantom’ dark energy wCDM model with constant equation of state $w = -1.2$. The cosmological parameters of these models have been calibrated against the cosmic microwave background (CMB) anisotropy power spectra from the WMAP 7 year data (Komatsu et al. 2011), these are summarised in Table 1. The values of $\Omega_m$ and $\sigma_8$ of the wCDM model have been chosen so as to be within the $1\sigma$ confidence region of the WMAP inferred constraints in the $\Omega_m - \sigma_8$ plane and along the degeneracy line of $\sigma_8 - w$ contours (see beige filled contours in Fig. 1). Henceforth, this model is statistically indistinguishable at $1\sigma$ from ΛCDM best-fit model to the WMAP 7 year data. The simulated models are the same as those selected for the DEUS-Full Universe Runs (Alimi et al. 2012; Rasera et al. 2014; Bouillot et al. 2015). These are also consistent with the constraints from the latest Planck primary CMB analysis within $2\sigma$ (TT, TE, and EE, the temperature and polarisation spectra Planck Collaboration VI 2020) (see red-yellow filled contours in Fig. 1). We notice that values of the equation of state varying in the range $-1.3 \le w \le -1$, being compatible with Planck data, are often considered in the literature (e.g., Alestas et al. 2020 and reference therein), for instance in the context of the ‘$H_0$ tension’ (Riess et al. 2019). Since these models are hard to discriminate at the homogeneous and linear perturbation level, it is important to find new cosmological probes, for instance in the non-linear regime of cosmic structure formation.

### 2.2. N-body simulations

The RAYGAL suite was run at the French TGCC and CINES super-computing centres. It consists of N-body simulations of the ΛCDM and wCDM models, respectively, spanning a volume of $(2625 \text{ h}^{-1} \text{ Mpc})^3$ sampled by $4096^3$ dark-matter particles. The optimal compromise between mass resolution (of order $\sim 2 \times 10^{10} \text{ h}^{-1} \text{ M}_\odot$ depending on the simulated cosmology) and large statistics makes them ideal to investigate the cosmology-dependent distribution and properties of halos from Milky Way

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** 1 and $2\sigma$ credible regions from the WMAP 7 yr data analysis (beige filled contours) and the Planck primary CMB analysis of the TT, TE, and EE anisotropy spectra (Planck Collaboration VI 2020) (red and yellow filled contours). The parameter values of the ΛCDM and wCDM models correspond to the dot and the cross markers, respectively.</figcaption></figure>

size to galaxy-cluster size (including galaxy-group size halos, which are the main deflectors in weak-lensing studies). The main novelty of these simulations is the relativistic ray tracing within the gravity light cones, which allows us to build realistic catalogues of dark-matter particles and halos that contain relativistic effects.

Initial conditions were generated with a modified version of MPGRAFIC (Prunet et al. 2008). In this code, the density field is generated as a Gaussian random field, while the associated displacement field is computed (in the version we developed) using second-order Lagrangian perturbation theory (2LPT). The starting redshift was chosen so as to ensure that the maximum displacement (among all particles) is of order one coarse cell (i.e. $\sim 0.64 \text{ h}^{-1} \text{ Mpc}$). The use of 2LPT minimises the effects of transients (Crocce et al. 2006), while assuming a late starting redshift minimises discreteness errors (Michaux et al. 2021).
<table><caption>**Table 2.** RAYGAL ΛCDM simulation parameters (wCDM simulation parameters are given in parentheses when they differ).</caption><thead><tr><th>L<sub>box</sub> (h<sup>-1</sup> Mpc)</th><th>n<sub>part</sub></th><th>n<sub>x</sub></th><th>n<sub>ref</sub></th><th>m<sub>ref</sub></th><th>n<sub>cell</sub></th><th>m<sub>p</sub> (h<sup>-1</sup> M<sub>☉</sub>)</th><th>Δ<sub>x</sub> (h<sup>-1</sup> kpc)</th><th>z<sub>i</sub></th><th>C<sub>dr</sub></th><th>P<sub>lin</sub>(k)</th></tr></thead><tbody><tr><td>2625</td><td>6.87 · 10<sup>10</sup></td><td>4096</td><td>7</td><td>8</td><td>4.41(4.46) · 10<sup>11</sup></td><td>1.88(2.01) · 10<sup>10</sup></td><td>5</td><td>46(51)</td><td>0.5</td><td>CAMB</td></tr></tbody></table>

**Notes.** $L_{\text{box}}$ is the box length, $n_{\text{part}}$ the number of particles, $n_x$ the grid size, $n_{\text{ref}}$ the number of refinements, $m_{\text{ref}}$ the refinement threshold, $n_{\text{cell}}$ the final number of AMR cells, $m_p$ the particle mass, $\Delta_x$ the spatial resolution, $z_i$ the starting redshift, $C_{\text{dr}}$ the Courant-like factor that determines the size of the integration time step (with respect to the largest stable time step), and $P_{\text{lin}}(k)$ specifies the source of the linear power spectrum used to generate initial conditions. All calculations have been performed in double precision.

Simulations were performed with (an optimised version of) the parallel AMR N-body code RAMSES<fnref n="2" /> (Teyssier 2002), which makes use of a multi-grid Poisson solver (Guillet & Teyssier 2011). A triangular shaped cloud (TSC) assignment scheme was used (instead of the default cloud-in-cell assignment) so as to increase the isotropy of the gravitational field. The characteristics of the simulations are summarised in Table 2.

While storing all particles and gravity cells at all time steps would be ideal, in practice it is unfeasible since it would overflow all the available storage devices (each snapshot occupies approximately 25 TB of memory). Therefore, we opted for a hybrid strategy and stored from 40 to 50 redshift snapshots (depending on the simulated cosmological model). On the other hand, particles and gravity cells are also stored at every coarse time step (thus ensuring good time resolution) in concentric shells located at the light-travel distance from an observer located at the centre of the box at $z = 0$: these are the background light-cone data. We saved three light cones with different depths (with $z_{\text{max}} = 0.5, 2$, and 10, respectively) and apertures (full sky, 2500, and 400 deg<fnref n="2" />, respectively) using the onion-shells techniques of Fosalba et al. (2008), Teyssier et al. (2009). Unlike many other studies, the maximum redshift $z_{\text{max}}$ of the cones was chosen in a conservative way to make sure that the comoving volume of the cone is smaller than the comoving volume of the simulation: this is required to minimise the effects of replica. The full-sky light cone does not rely on any replica techniques, while the narrow light cones are tilted with respect to the x-axis of the simulations (by 17 or 25 degrees along the two spherical coordinates angles depending on the light cones) to minimise replica effects. Rather than storing the AMR cells that are the closest to the past null Friedmann-Lemaître-Robertson-Walker (hereafter FLRW) light cone of the observer, we stored two cells at each spatial location: one slightly ahead of time and one slightly behind in time. Thanks to this ‘double-layer’ strategy, it is then possible to explore different approaches for building the light cones: by taking the closest cell (in time), by averaging between the two cells, or by interpolating between the cells. It is also possible to compute time derivatives, such as the derivative of the gravitational potential $\phi$, the key quantity of integrated Sachs-Wolfe-Rees-Sciama (ISW-RS) effect. It is also worth noticing that, in principle, light propagates on the null real light cone and not on the null FLRW light cone (i.e. at a given spatial location the gravitational field can be different for different light rays because of the variation in the potential during the time delay between light-ray arrivals). While in this article we neglect this effect in our ray tracing, very advanced studies could in principle account for this very subtle effect by interpolating between the two cells for each individual light ray.

In RAMSES, each task from the MPI library writes a binary file that corresponds to a segment of the Peano-Hilbert space filling curve. Since these numerous raw binary files are not user

friendly, they are rewritten as a small number of HDF5 files (with the tools provided by PFOF<fnref n="3" />). For snapshots, each HDF5 file corresponds to a cubic sub-volume of the whole simulation. Particle files contain information about the dark matter particles. For light cones, each HDF5 file corresponds to a given shell. The particle files contain information about particles while the gravity files contain information about the AMR cells. Within each file, particles and cells are reorganised in small cubes (corresponding to an HDF5 group) to make further selection faster. Finally, the HDF5 files are self-documented to facilitate their usage.

Halos were detected in all the snapshots using a parallel friend-of-friends algorithm (as implemented in pFoF as described in Roy et al. 2014) with various linking lengths ranging from $b = 0.05$ to $b = 0.4$. A variety of classical halo properties were computed and the constituent halo particles were also stored for each halo in the catalogues. Since the volume spanned by halos defined with $b = 0.4$ is much larger than the volume of most commonly defined halos, users are free to re-compute their favourite halo definition from halo particles and to evaluate any physical quantity of interest (with the parallelisation on a per halo basis being straightforward). We also provide halo properties for halos detected with a spherical overdensity halo finder using a range of overdensity thresholds (with respect to the mean matter density) from $\Delta_m = 50$ to $\Delta_m = 12800$ (as well as other popular definitions based on critical density or virial overdensity). Halos are also detected on the light cone with pFoF for some selected linking lengths, while the rest of the detections are currently in progress.

## 2.3. Relativistic ray tracing and catalogue creation

To produce realistic observables, we use the ray-tracing framework MAGRATHEA-PATHFINDER (Breton & Reverdy 2022) based on the AMR library MAGRATHEA (Reverdy 2014), which computes the propagation of light rays within the 3D AMR structure of light cones (such as the ones described in Sect. 2.2).

We consider a weakly perturbed FLRW metric given by

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

with $\phi$ the gravitational potential, $\eta$ the conformal time, $x^i$ the comoving position and $c$ the speed of light. The integration is performed backwards in time, starting from the observer at $z = 0$. Photons are initialised with $k^0 = 1$ and $k_\nu k^\nu = 0$, where $k^\alpha = (d\eta/d\lambda, dx^i/d\lambda)$ and $\lambda$ is the affine parameter. We interpolate the gravity information from the light cone to the photon position using an inverse TSC scheme (so that to be consistent with the interpolation procedure used by the N-body solver) at the finest refinement level that contains the photon. If there are not enough neighbours at the same level, we perform an interpolation on a

<footnote n="2">² https://bitbucket.org/rteyssie/ramses/src/master/</footnote>

<footnote n="3">³ https://gitlab.obspm.fr/roy/pFoF</footnote>
coarser grid and eventually repeat this procedure until we find
enough neighbours to compute the interpolation. If we cannot
interpolate at the coarse level, the propagation stops (it means
that the photon is out of the numerical light cone). Once all the
gravity information at the photon location has been evaluated,
we solve the geodesic equations:

$$
\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 (2)
$$

$$
\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 (3)
$$

with $\phi$ the gravitational potential and $a' = da/d\eta$ the time derivative of the scale factor. In particular, we use the precise value of $\partial\phi/\partial\eta$ estimated from the double-layer strategy adopted to build the light cones, thus ensuring that the ISW-RS effect is correctly implemented even in the strongly non-linear regime.

Our method fully takes advantage of the adaptive mesh. In
particular, we use an adaptive step equal to one-fourth of the
finest cell at the photon location. This allows for a precise light
propagation in high-density regions.

To emulate what we really observe, we found the null
geodesics connecting the observer to each individual source (the
dark matter particle, halo, etc.). To do so, we used the methodol-
ogy described in Breton et al. (2019), which consists of the fol-
lowing steps. First, we launch a light ray towards the comoving
direction of the source. However, due to deflections induced by
the matter-density field, the photon will probably not hit the
source. If that is the case, we then iterate on the initial conditions
using a Newton-like method until the angle subtended by the
source’s comoving position and the photon at the same distance
is smaller than a given threshold (in our case, set to 10 milliarc-
second). Once the null geodesic connecting the observer and a
source is found, we estimate all the contributions to the observed
redshift at first order in metric perturbations, accounting for local
and integrated terms using either a linear decomposition of the
redshift or its exact definition. We also provide several (partial)
redshift definitions that progressively incorporate part of the rel-
ativistic effects in order to disentangle the various contributions<fnref n="4" />:

$$
z_0 = \frac{a_0}{a} - 1, \tag{4}
$$

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

$$
z_2 = \frac{a_0}{a} \left( 1 + \frac{\mathbf{v}_s \cdot \mathbf{n}}{c} + \frac{\phi_o - \phi_s}{c^2} \right) - 1, \quad (6)
$$

$$
z_3 = \frac{a_0}{a} \left( 1 + \frac{\mathbf{v}_s \cdot \mathbf{n}}{c} + \frac{\phi_o - \phi_s}{c^2} + \frac{1}{2} \left(\frac{v_s}{c}\right)^2 \right) - 1, \quad (7)
$$

$$
z_4 = \frac{a_0}{a} \left( 1 + \frac{\mathbf{v}_s \cdot \mathbf{n}}{c} + \frac{\phi_o - \phi_s}{c^2} + \frac{1}{2} \left(\frac{v_s}{c}\right)^2 - \frac{2}{c^2} \int_{\eta_s}^{\eta_o} \dot{\phi} d\eta \right) - 1, \quad (8)
$$

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

where $g_{\mu\nu}k^{\mu}u^{\nu} = -ack^0\left(1 + \phi/c^2 + \mathbf{v} \cdot \mathbf{n}/c + \frac{1}{2}\left(\frac{\mathbf{v}}{c}\right)^2\right)$, the subscripts 'o' and 's' denote an evaluation at the observer and the source, respectively, and $z_n$ gives the different redshifts provided for each single source. $z_0$ is the comoving redshift (from background expansion), $z_1$ also includes gravitational redshift, $z_2$

includes the Doppler effect and $z_3$ adds the transverse Doppler effect. The $z_4$ includes the contribution to the final redshift from all the effects, including the ISW-RS effect. The $z_5$ is the exact redshift definition and should match $z_4$ (except for some small higher-order contributions).

Furthermore, given the fact that the observed angular posi-
tion of the source, θ, is known, as is the comoving angular posi-
tion, β, we can compute the deflection angle and evaluate the
lensing distortion matrix, 𝓅, along the true photon trajectory.
This reads as

$$
\mathcal{A} \equiv \frac{\partial \beta}{\partial \theta} = \begin{pmatrix} 1 - \kappa - \gamma_1 & -\gamma_2 + \omega \\ -\gamma_2 - \omega & 1 - \kappa + \gamma_1 \end{pmatrix}, \quad (10)
$$

where $\kappa$ and $\gamma = (\gamma_1, \gamma_2)$ are the weak-lensing quantities related, respectively, to the convergence and shear, while $\omega$ is the (small) image rotation. This matrix can be estimated using either an infinitesimal-beam or finite-beam methods (both are described in Breton & Fleury 2021). For the present release, we use the finite-beam method with a bundle semi-aperture of $10^{-5}$ radians (which corresponds for instance to the apparent angle of a ruler of comoving size $\sim 10 h^{-1} \text{ kpc}$ located at a comoving distance of $\sim 1 h^{-1} \text{ Gpc}$ in a flat universe). It is important to stress that there are two possible ways to define the beam cross-sectional area: perpendicular to the radial direction or perpendicular to the light-ray wave vector. Here we consider the second option (i.e. the so-called Sachs basis), which is more closely related to observable quantities: we account for the tilt of the cross-sectional area.

One has to be careful since the link between these lensing
quantities and actual observables is not necessarily straightfor-
ward. The magnification µ is directly given by the inverse of
the determinant of the distortion matrix. For a comoving source,
this gives the amplification of the flux by weak lensing, where
the flux is indeed an observable. In virtue of the distance-duality
relation, µ also provides the fractional change in the apparent
area of the source where the area is also an observable. How-
ever, for a non-comoving source, one should also account for
redshift perturbations to evaluate the change in flux and area:
this is called Doppler convergence (Bonvin 2008). The relative
variation in the flux is given by the relative variation (induced
by the redshift perturbations) of the square of the luminosity
distance. These effects play a subdominant role at high redshift
(z > 0.5) but may play a role at lower redshift. These subtleties
are well described in Breton & Fleury (2021). Another observ-
able in weak-lensing studies is the apparent shape of the sources
(i.e. their ellipticity e). The average ellipticity is related to the
reduced shear g = γ/(1 − κ) as ⟨e⟩ = ⟨g⟩ (if κ ≪ 1 one finds
⟨e⟩ ≈ ⟨γ⟩). The role of redshift perturbations on the shape is
negligible. Textbooks often assume that the shear power spec-
trum and the gravitational convergence power spectrum are the
same (thus highlighting the central role of the convergence in
weak-lensing studies): while this is true at first order, we discuss
this point in this article. Finally, the orientation of the source
is an observable that is a probe of the rotation. However, it is
a very subtle effect and it is also very challenging to know the
original orientation of the sources. There are however some pro-
posals using the polarisation of radio-galaxies (Francfort et al.
2021). Here we provide both the fundamental lensing quanti-
ties and the detailed redshift perturbations so that the user can
reconstruct any physical or observable quantities of interest. In
any case, the gravitational convergence κ play an important role
since κ > |γ| > ω and the Doppler corrections are subdominant at
z > 0.5, this is why we focus on this key quantity in this article.

It should be noted that in the current implementation, we
only detect one image per source. This is a valid approximation

<footnote n="4">⁴ Although we implemented the possibility to give any peculiar veloc-
ity to the observer $v_0$, for simplicity we decided to set $v_0 = 0$ in this
release.</footnote>
in the weak-lensing regime (although convergence and shear are not necessarily small). A possible extension of our work would be the use a multiple-root finder to account for multiple images. However, there are multiple reasons as to why the weak-lensing regime should still be a reasonable approximation for RAYGAL. First, for most practical cases the number of strong lensing events for galaxy-size lenses is small. Even if the probability to have strong-lensing events increase with redshift, for a realistic radial selection function the number of observable sources heavily decreases beyond $z \approx 2-3$. Second, strong lensing is usually relevant at very small angular scales, which is not necessary the focus of the present simulation. Third, when we perform source-averaging, the tail of the magnification's probability distribution function is damped with respect to the directional-averaging case (see e.g., Takahashi et al. 2011; Breton & Fleury 2021, and references therein). Moreover, this effect is further enhanced when using the finite-beam method (Fleury et al. 2017, 2019; Breton & Fleury 2021). Finally, even if strong lensing events occur, our geodesic finder would most likely find the principal image of the source.

## 2.4. Maps

Maps were built by projecting the matter distribution and its properties onto a sphere. We discretised the maps by considering the HEALPIX (Hierarchical Equal Area isoLatitude PIXelation) pixelation scheme (Górski et al. 2005).

The matter distribution can be weighed by radial selection functions of different type. Here, we generated two type of maps: (i) Dirac radial selection maps, HEALPIX maps generated assuming a Dirac selection function, and (ii) top-hat radial selection maps, HEALPIX maps generated assuming a top-hat selection function.

Dirac radial selection maps are computed by stopping all the light rays at a fixed comoving distance from the observer. The quantity of interest (for instance the distortion matrix) is then interpolated from the closest values on the light cone. Top-hat radial selection maps are built from the relativistic catalogues. Particles (or halos) are selected within a shell of mean redshift $z_{\text{mean}}$ and (half-)width $\Delta z$. Particles (or halos) are then sorted according to Healpix while keeping track of the empty pixels (no particles) and the masked pixels (outside of the cones). The surface overdensity is evaluated by computing the excess of surface density in pixels with respect to the mean surface density of the shell in the unmasked region. The distortion matrix is evaluated as an average over all particles. Empty cells are filled with the default undefined value provided by Healpix.

## 2.5. RAYGAL data: Snapshots, light cones, relativistic catalogues, and maps

RAYGAL data are divided into four categories: snapshots, FLRW light cones, relativistic catalogues, and maps. Most of the RAYGAL data are available on the website of the COS team at LUTH<fnref n="5" />.

Snapshot data, corresponding to particles at a given time, consist of: (i) snapshots of all particles at a given time (position, velocity, and identity), (ii) all particles within halos (position, velocity, and identity), and (iii) the properties of all the halos (identity, number of particles, most-bounded particle position,

centre of mass position, mean velocity, maximum radius, moment of inertia, circular velocity, velocity dispersion, angular momentum, kinetic energy, self-binding energy, density profile, and circular velocity profile).

The FLRW light-cone data (corresponding to particles and cells near the past FLRW light cone) consist of: (i) the particle light cone, all particles just before or just after the past null light cone (position, velocity, redshift) and (ii) the gravity light cone, all gravity cells just before and just after the past null light cone (position, potential, gravitational field, redshift). These snapshot and light-cone data are built from the simulations and post-processing as described in Sect. 2.2.

Relativistic catalogues, corresponding to sources as seen by the observer and accounting for relativistic effects such as weak lensing and RSDs, consist of: (i) the halo catalogue (un-lensed and lensed angular positions, background, and perturbed redshifts (from Eq. (4) to Eq. (9))), weak-lensing distortion matrix, and halo mass), and (ii) the particle catalogue, a random subset of dark matter particles (same properties as for the halo catalogue, except mass). They represent a population of sources with a bias $b = 1$. These catalogues are built from the ray-tracing techniques described in Sect. 2.3.

Maps, corresponding to the projection of the matter distribution and its properties, consist of Dirac radial selection maps ($\kappa$, $\gamma_1$, $\gamma_2$, and $1/\mu$, respectively the convergence, shear components, and inverse magnification).

We provide Dirac radial selection maps for several redshifts within the range of our light cones. Top-hat radial selection maps (for $\kappa$, $\gamma_1$, $\gamma_2$ and $\delta$) can be computed from the released catalogues using HEALPIX tools. For the narrow cones, pixels outside the footprint are masked. In that case, we only provide the relevant pixels within the angular selection function, therefore producing lighter files.

The characteristics of RAYGAL data are summarised in Tables 3 and 4. An illustration of the apparent matter distribution and the various relativistic effects within the light cone 'Narrow 2' are shown in Figs. 2 and 3⁶. In this cone, the source angular density (~40 arcmin⁻²) is similar to the one expected from present and future photometric surveys such as DES, KIDS, HSC, *Euclid*, and LSST while the halo density (~10⁻³ h³ Mpc⁻³) is also similar to the one expected in spectroscopic surveys such as BOSS, eBOSS, DESI, *Euclid*, and SKA.

## 2.6. Case study: Configuration, definition, angular cross-spectrum evaluation, and analytical predictions

As already mentioned, in this work we present a case study application of the RAYGAL light-cone datasets to the analysis of relativistic effects on the angular cross-spectra. In particular, we correlate the key quantity of clustering (i.e. the matter overdensity, $\delta$) and that of lensing (i.e. the gravitational convergence, $\kappa$) using HEALPIX maps. It is worth noting that these are not directly observable quantities (such as the galaxy overdensities or the galaxy ellipticities), nonetheless they are the central quantities in clustering and lensing studies from which one can deduce the correlation between observables. The goal here is not to focus on specific observables for a given survey, rather to provide general results about the imprint of relativistic effects on the density and lensing angular power spectra and cross-spectra.

<footnote n="5">⁵ https://cosmo.obspm.fr/public-datasets/</footnote>

<footnote n="6">⁶ We made use of Py-SPHViewer (Benitez-Llambay 2015) for the realisation of these images.</footnote>
<table><caption>**Table 3.** Data available for RAYGAL snapshots in ΛCDM and wCDM cosmologies.</caption><thead><tr><th>Data | z</th><th>z<sub>i</sub></th><th>3.00</th><th>2.33</th><th>2.00</th><th>1.50</th><th>1.00</th><th>0.67</th><th>0.50</th><th>0.43</th><th>0.25</th><th>0.11</th><th>0.00</th><th>Other z</th></tr></thead><tbody><tr><td>Part.</td><td>x</td><td></td><td></td><td>x</td><td></td><td>x</td><td></td><td>x</td><td></td><td></td><td></td><td></td><td>x</td></tr><tr><td>Halos part.</td><td></td><td>x</td><td>x</td><td>x</td><td>x</td><td>x</td><td>x</td><td>x</td><td>x</td><td>x</td><td>x</td><td>x</td><td>x</td></tr><tr><td>Halos prop.</td><td>x</td><td>x</td><td>x</td><td>x</td><td>x</td><td>x</td><td>x</td><td>x</td><td>x</td><td>x</td><td>x</td><td>x</td><td>x</td></tr></tbody></table>

**Notes.** We note that halo properties are available for about 30–40 ‘other z’ per cosmology, but the redshifts may differ from one cosmology to another. *z* is the redshift, *z*<sub>i</sub> the starting redshift, ‘Part.’ stands for particles, and ‘prop.’ stands for properties.

<table><caption>**Table 4.** Data available for RAYGAL light cones in ΛCDM and wCDM cosmologies (wCDM simulation parameters are given in parentheses when they differ).</caption><thead><tr><th>Type | charac.</th><th>z<sub>max</sub></th><th>Aperture (deg<sup>2</sup>)</th><th>Gravity</th><th>part.</th><th>rel. cat. part(#)</th><th>rel. cat. halo (#)</th><th>Double layer</th></tr></thead><tbody><tr><td>Narrow 1</td><td>10.00</td><td>400</td><td>x</td><td>x</td><td>10<sup>8</sup></td><td>2.5 × 10<sup>6</sup> (on going)</td><td>x</td></tr><tr><td>Narrow 2</td><td>2.00</td><td>2500</td><td>x</td><td>x</td><td>3.5 × 10<sup>8</sup></td><td>1.2 × 10<sup>7</sup> (1.3 × 10<sup>7</sup>)</td><td>x</td></tr><tr><td>Full sky</td><td>0.48</td><td>41253</td><td>x</td><td>x</td><td>2 × 10<sup>8</sup></td><td>1.3 × 10<sup>7</sup> (1.5 × 10<sup>7</sup>)</td><td>x</td></tr></tbody></table>

**Notes.** ‘Charac.’ stands for characteristics, *z*<sub>max</sub> is the maximum redshift, aperture is the angular aperture, gravity indicates the presence of gravity information, ‘part.’ indicates the presence of particle information, ‘rel. cat. part(#)’ indicates the number of particles in the relativistic catalogues, and ‘rel. cat. halo (#)’ indicates the number of halos. ‘Double-layer’ indicates the presence of two light cones so as to compute time derivatives or perform time interpolation. The particle density in cone Narrow 2 is similar to the galaxy density of present and future photometric surveys, while the halo density is similar to the galaxy density of present and future spectroscopic surveys.

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Apparent matter overdensity within the RAYGAL intermediate light cone called ‘Narrow 2’. This includes ‘all’ weak-field relativistic contributions.</figcaption></figure>

### 2.6.1. Two-shell configuration

The matter overdensity and the gravitational convergence are averaged within 3D pixels defined by the radial extension of the shell and the HEALPIX pixelation scheme for the angular extension. The number of pixels is equal to $12 \times N_{side}^2$ and we chose $N_{side} = 4096$ for top-hat radial selection map. Given the number of sources in each shell, using a larger $N_{side}$ increases the

number of empty pixels and the gain in precision is limited. For test maps under Born approximation, we first computed a map with $N_{side} = 8192$ (resulting in a total number of light rays of the same order of magnitude as the number of sources in catalogues shells), which was then degraded (i.e. each pixel is computed as an average over the 4 smaller pixels within it) down to $N_{side} = 4096$ so as to use the same resolution in all maps.
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Various contributions to the apparent matter overdensity within the RAYGAL intermediate light cone: Doppler effect (top left), weak lensing (top right), gravitational redshift (bottom left), and ISW-RS effect (bottom right).</figcaption></figure>

We considered two shells (1 and 2) extracted from the $\Lambda$CDM 2500 deg$^2$ light cone. We isolated two possible redshifts for each shell: $z = 0.7 \pm 0.2$ and $z = 1.8 \pm 0.1$. We chose $z = 1.8$ since it is one of the largest possible redshift in the cone while the comoving distance to $z = 0.7$ is about half the comoving distance to $z = 1.8$, thus corresponding to a maximum of the lensing kernel for a source located at $z = 1.8$. Given that relativistic effects correlate distant shells, we investigated all the ten possible auto- and cross-power spectra for these two shells. We first investigated the clustering (three density-density spectra), then the weak lensing (three convergence-convergence spectra), and finally the galaxy-galaxy lensing (four density-convergence spectra).

### 2.6.2. Relativistic effect definition

Weak lensing as a relativistic effect has already been widely studied in the literature, especially under the Born approximation (i.e. by computing lensing quantities and in particular the convergence as an integral of the transverse Laplacian along the un-deflected light-path). Henceforth, in the context of weak lensing, we use as a reference case the convergence computed using the Born approximation ($\kappa_{\text{Born}}$). The mean Born convergence for a population of source galaxies (see Kilbinger 2015) is given by

$$ \kappa_{\text{Born}}(\boldsymbol{\theta}) = \int_{0}^{\chi_{\text{lim}}} n(\chi_s) d\chi_s \int_{0}^{\chi_s} d\chi \frac{\chi_s - \chi}{\chi_s} \frac{\Delta_{\perp} \phi(\chi_{\theta}, a(\chi))}{c^2} \chi, \quad (11) $$

with $\chi$ the comoving distance (in a flat universe), $\chi_s$ the comoving distance of the source, $\chi_{\text{lim}}$ the limiting comoving distance of the galaxy sample, $n(\chi_s)$ the source probability distribution, and $\Delta_{\perp} \phi$ the Laplacian perpendicular to the line of sight. We considered a top-hat radial selection function (i.e. a shell). For a shell of volume $V_{\text{shell}}$ and constant density, $n(\chi_s) \propto 4\pi\chi_s^2/V_{\text{shell}}$. It is worth noting that a common additional approximation is to assume

$$ \Delta_{\perp} \phi = \frac{3H_0^2 \Omega_m \delta}{2a}, \quad (12) $$

with $H_0$ the Hubble constant. We do not make this assumption in our Born calculation. Moreover in this expression $\delta$ is often estimated within shells of typical size between 10–100 $h^{-1}$ Mpc. Here we estimate $\Delta_{\perp} \phi$ at the resolution of the simulation (5–600 $h^{-1}$ kpc). We checked that using such approximations leads to similar results in our article. A detailed study of the impact of these two approximations on various statistics would deserve a dedicated paper. Our reference is therefore an accurate Born convergence calculation along the un-deflected path.

In the context of clustering or number count, our reference is the comoving matter overdensity $\delta_{\text{com}}$. The overdensity is estimated as the local excess of particles (in real space) with respect to the mean matter density,

$$ \delta_{\text{com}}(\boldsymbol{\theta}) = \frac{N_{\text{com}}(\boldsymbol{\theta})}{\langle N_{\text{com}}(\boldsymbol{\theta}) \rangle} - 1, \quad (13) $$
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** HEALPIX maps from RAYGAL. Left: overdensity map at $z = 0.7$ including relativistic effects. Right: convergence map at $z = 1.8$ including relativistic effects. The size of the field extracted from cone Narrow 2 is about 2500 deg$^2$. A gnomonic projection has been used, and the maps have been smoothed with a Gaussian beam of 10 arcmin full width at half maximum for representation purposes. We note that some large overdensities at $z = 0.7$ can be seen in the convergence map at $z = 1.8$.</figcaption></figure>

where $N_{\text{com}}(\theta)$ is the number density of source. Here again, light rays are assumed to propagate in a straight line. The Doppler effect can be considered as one of the standard effect in RSD studies (involving 3D correlation). However, the Doppler effect is often neglected in the context of angular correlation (see Grasshorn Gebhardt & Jeong 2020) as well as in weak-lensing and galaxy-galaxy lensing studies (see Ghosh et al. 2018).

We therefore call ‘relativistic effects’ all the non-trivial GR effects beyond these two references ($\delta_{\text{com}}$ and $\kappa_{\text{Born}}$) cases. To gain further understanding on the nature of relativistic effects, we further decomposed them into two pieces.

The first, MB effects, are related to the decrease or increase in the apparent source overdensity due to weak-lensing magnification or de-magnification as introduced by Turner et al. (1984). The apparent density is related to the comoving density by $1 + \delta = (1 + \delta_{\text{com}})\mu^{-1}\mu^{2.5s}$ with $\mu$ the magnification and $s$ the logarithmic slope of the luminosity function. The first term $\mu^{-1}$ is called the dilution term and is always present. The second term is only present for flux limited sample. Because we consider a volume and mass-limited sample of sources, we are only sensitive to the dilution effect (corresponding to a flux limited sample with a logarithmic slope of the luminosity function of $s = 0$). This also biases the estimate of the convergence of a sample of galaxies as highly magnified region (with large convergence) are poorly sampled while de-magnified region are very well sampled. Moreover, widely used source-weighted estimator are even more biased by this effect. We also include in this category other (subdominant) weak-lensing effects, such as post-Born effects related to the propagation of light, which does not happen in a straight line.

The second, RSDs, are related to redshift perturbations, which cause a change in the apparent density. This includes the Doppler effect, gravitational redshift effect, transverse Doppler

effect, and ISW-RS effect. We note that lensing deflections effect are already included in MB.

All these effects do not exactly add up linearly. We therefore did not include them one by one (e.g., none, MB alone, RSDs alone) but progressively (e.g., none, MB, and MB+RSDs). We note that even though the way to decompose the full signal in various components depends on which phenomenon we are interested in, the observed signal including all contributions is unique and well defined.

An illustration of the overdensity and convergence within the two shells and the importance of relativistic effects is shown Figs. 4 and 5.

### 2.6.3. Angular cross-spectrum estimator

The angular cross-spectra between two Healpix maps is estimated with POLSPICE<fnref n="7" /> (Szapudi et al. 2001; Chon et al. 2004). The spectra are computed with fast spherical harmonic transforms. They are corrected from angular selection function effects (i.e. mask effects) as well as pixel effect (i.e. deconvolution from the pixel window function)<fnref n="8" />. The density auto-spectra are polluted by shot noise at large $\ell$. To correct this effect, we generated a noise map with the same number of randomly placed particles as in the data from which we compute the shot noise spectra. We then subtract these shot noise spectra from the measured spectra. This is more accurate than using the analytical formula since it removes the possible non trivial bias from the estimator in presence of a mask. In the context of weak-lensing and galaxy-galaxy lensing studies, we follow Schmidt et al. (2009).

<footnote n="7">⁷ http://www2.iap.fr/users/hivon/software/PolSpice/</footnote>

<footnote n="8">⁸ In practice, for the Narrow 2 cone we used the options: -apodizetype 0 -apodizesigma 25. -thetamax 50. -pixelfile YES.</footnote>
<figure><img src="image_6.png" /><figcaption>Fig. 5. Relativistic contributions in the maps shown in Fig. 4. Top left: relative difference between the overdensity maps at z = 0.7 with and without MB. Top right: relative difference between the overdensity maps at z = 0.7 with and without RSDs. Bottom left: relative difference between the convergence maps at z = 1.8 with and without the MB effect. As expected, the MB effect is anti-correlated with the value of the convergence. Bottom right: (Small) relative difference between the convergence maps at z = 1.8 with and without RSDs. Relativistic effects play a non-trivial role at the ten percent level in the density and convergence maps.</figcaption></figure>

Since we use a pixel-based estimator we weight each κ pixel with the density 1 + δ (this is equivalent to an inverse variance weight). As shown in Appendix B, such a weight allows a good match with the direct pair-based estimator.

### 2.6.4. Angular cross-spectrum predictions (with CLASS)

Analytical predictions usually rely on the linear mapping assumption. Especially in the case of the Doppler contribution at small scale this is an approximation. However, non-linear corrections are less visible than in RSD studies because here we

consider the angular power spectrum (which is indirectly sensitive to radial motions in a way that depends on the thickness of the shell). Therefore, under this assumption one can decompose the overdensity as

$$ \delta_i \approx \delta_i^{\text{com}} + \delta_i^{\text{MB}} + \delta_i^{\text{RSD}}, \quad (14) $$

where $\delta_i^{\text{com}}$ is the comoving overdensity, $\delta_i^{\text{MB}} \approx -2\kappa_i$ is the MB contribution dominated by the dilution term ($\kappa_i$ is the gravitational convergence and we assumed a logarithmic slope $s = 0$ for the source count in the context of our mass selected sample), and $\delta_i^{\text{RSD}}$ is the contribution from redshift perturbations. In the
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** RAYGAL real-space matter power spectra. Left: 3D matter power spectra at $z=0$ (top lines in black and red) and $z=1$ (bottom lines in green and blue) for $\Lambda$CDM (continuous lines) and wCDM (dashed lines) cosmologies. RAYGAL power spectra are shown in black ($z=0$) and green ($z=1$), and emulator power spectra are shown in red ($z=0$) and blue ($z=1$). Right: relative deviations between simulations and emulator at $z=0$ (black) and $z=1$ (green) in $\Lambda$CDM (continuous lines) and wCDM (dashed lines) cosmologies. Percent level agreement is shown up to $k=1-2h\,\text{Mpc}^{-1}$, where resolution effects damp the power spectrum. Baryonic effects are also expected to damp the power spectrum near this wavenumber.</figcaption></figure>

## 3.2. Clustering: Density-density angular spectra

In this section we focus on the angular matter clustering, $C_{\delta_1\delta_2}(\ell, z_1, z_2)$, for shells 1 and 2 at redshift $z_1 \in \{0.7, 1.8\}$ and $z_2 \in \{0.7, 1.8\}$, respectively.

### 3.2.1. Comoving density

The comoving matter-density angular auto- and cross-spectra are shown Fig. 7, left panel. The spectra are binned with 30 bins in the multipole interval $30 \le \ell \le 5000$. Symbols represent the mean value of the spectra in the bin, while the error bars represent the standard deviation of the spectra within the same bin. The continuous lines show the analytical predictions. The auto-spectra $C_{\delta_1\delta_2}(\ell, z_1 = 0.7, z_2 = 0.7)$ and $C_{\delta_1\delta_2}(\ell, z_1 = 1.8, z_2 = 1.8)$ show an excellent agreement with the CLASS analytical prediction across the entire multipole interval. At large scale ($\ell < 50$), the fluctuations we see are very likely to be related to the finite size of the light-cone angular aperture (sample variance). The good agreement in the linear regime between $\ell = 50$ and $\ell = 500$ validates our simulation measurement. The non-linear regime begins near $\ell \sim 500$, corresponding approximately to $k \sim 0.2h\,\text{Mpc}^{-1}$ at $z \sim 1$, where $P(k)$ is known to deviate from linear theory at the ~10% level (see e.g., Rasera et al. 2014). As we calibrated the $P(k)$ used in the analytical predictions to that from the RAYGAL simulation, the agreement still holds in the non-linear regime between $\ell = 500$ and $\ell = 5000$. We halt at $\ell = 5000$, which corresponds approximately to $k \sim 2h\,\text{Mpc}^{-1}$ at $z = 1$, where finite resolution effects in our simulation (and possibly baryonic effects in the Universe) can no longer be neglected. Since the maximum $\ell$ we consider is of order of $N_{\text{side}}$, we also expect the aliasing effect to be small. The green curve shows the cross-spectrum $C_{\delta_1\delta_2}(\ell, z_1 = 0.7, z_2 = 1.8)$, which is completely negligible in real space both in numerical data and analytical prediction. This is because the shells are separated by about $1750\,h^{-1}\,\text{Mpc}$. The right panel shows the relative deviation from CLASS. It illustrates the agreement with the analytical prediction to be better than the 5% level and within error bars for the auto-spectra. The cross-spectrum lies outside the plot because of a division by a nearly zero CLASS prediction. Overall, this is a good cross-validation of CLASS and RAYGAL's random subset matter angular spectra.

### 3.2.2. Relativistic effects

In Fig. 8, we plot the matter overdensity auto- and cross-spectra in the presence of the relativistic effects (upper left panel). We can see a noticeable change compared to the case without relativistic effects. shown in Fig. 7. The agreement between RAYGAL measurements and CLASS prediction (including RSDs and MB) is as good as before from the linear to the non-linear regime. The cross-spectrum between the $z=0.7$ and $z=1.8$ shells is now very different from zero: it reaches about 10% of the auto-spectra, which is not negligible at all. This is due to weak lensing that couples the clustering signal between two density shells: the density in the low-redshift shell is a lens for the density of the high-redshift shell. Here again there is a good agreement between predictions and simulations except below $\ell = 50$ where the simulation spectrum strongly underestimates the analytical predictions. This is likely related to the limited aperture of the narrow light cone, which affects this subtle correlation.

The other panels allow us to investigate relativistic effects in more details: they show the relativistic contributions normalised by the CLASS spectra including all relativistic effects (we chose this normalisation since this is a smooth and non-zero spectrum in all our plots). The upper right and bottom-left panels show the relativistic contributions to the auto-spectra. The MB contribution is mostly constant of order +3% at $z=0.7$ and -3% at $z=1.8$. This signature of lensing on the matter distribution is a subtle effect, well captured by RAYGAL, in good agreement with CLASS expectation from $\ell=50$ to $\ell=5000$. This effect has been investigated in several works (e.g., Fosalba et al. 2015 and reference therein). It is related to the so-called MB, where the apparent density becomes (at first order) $\delta_i = \delta_i^{\text{com}} + (5s-2)\kappa_i$ (with $s$ the logarithmic slope of the source number count, which is $s=0$ for our mass limited sample; in this peculiar case it is also called dilution bias). As a consequence, there is an additional lensing contribution to the matter-density angular power spectrum (see the second line of Eq. (A.2)) that can be simplified as

$$C_{\delta_1\delta_2}(l) \approx C_{\delta_i^{\text{com}}\delta_j^{\text{com}}}(l) - 2C_{\delta_i\kappa_j}(l) - 2C_{\kappa_i\kappa_j}(l) + 4C_{\kappa_i\kappa_j}(l). \quad (16)$$

The RSD contributions to the multipoles of the 3D galaxy-galaxy two-point correlation function is important in the linear regime (the Kaiser effect, which corresponds to an amplification of the clustering due to coherent motions plus a quadrupole
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Comoving matter-density angular auto-spectrum in a single shell at $z = 0.7$ (red) and $z = 1.8$ (pink) in ΛCDM cosmology. The cross-spectrum between two shells at $z = 0.7$ and $z = 1.8$ is shown in green. **Left:** spectrum from the RAYGAL 2500 deg$^2$ light cone (diamonds with error bars) and the spectrum from CLASS (dashed line). We note that the CLASS prediction for the cross-spectrum is not visible since it is oscillating around zero with a small amplitude, reaching at maximum $10^{-10}$ at $\ell = 30$. The RAYGAL cross-spectrum is also oscillating around zero but with a larger amplitude (very likely due to sample variance), thus explaining the missing points, which are negative (or smaller than $10^{-9}$). **Right:** relative deviation from the CLASS prediction. The cross-spectrum is nearly zero in real space, and the relative deviation is therefore outside the graph, with very large error bars (that are omitted for readability). Overall, we find a good agreement for the auto-spectrum, except at large scales (small $\ell$) due to the finite area of the light cone.</figcaption></figure>

and hexadecapole component) and in the non-linear regime (the finger-of-God effect, which corresponds to a damping of the small-scale clustering due to incoherent motions). These contributions are however much smaller when using large spherical shells (>100 $h^{-1}$ Mpc) centred around the observer. Redshift perturbations mostly modify the apparent radial distance of the sources and therefore do not have a strong effect on the angular power spectra. There is however, a subtle and non negligible amplification (see the first term of the third line of Eq. (A.2)) due to the Kaiser effect visible at large scale below $\ell \sim 100$. For the shell thickness we use, it can reach 5% at $\ell = 20$ (for the $z = 0.7 \pm 0.2$ spectra) and 20% at $\ell = 90$ (for the $z = 1.8 \pm 0.1$ spectra). While the geometry makes the effect non-trivial, it is well reproduced by the simulation on agreement with CLASS predictions. The MB effect is also well captured up to $\ell \approx 1000$ using our MB estimate based on Born inverse magnification weight (we discuss the origin of the small discrepancy at large $\ell$ and low redshift in the lensing section). For both shells, the RSD effect is too small at large $\ell$ and we cannot see the damping behaviour of the spectra. In Sect. 3.5.1, we performed again the same exercise but at lower redshift, for narrower shells using the full-sky data. The damping of the angular spectra (related to the fingers-of-God effect) now becomes apparent and is not well reproduced by CLASS.

The bottom-right panel shows the relativistic contribution to the density cross-spectrum between two shells. Because the real-space spectrum is nearly zero, the relativistic contribution is compatible with 100% for both CLASS and RAYGAL (between $\ell = 50$ and $\ell = 5000$). We note that we normalised the spectrum variation by the CLASS spectrum (even for RAYGAL spectrum). The good agreement between CLASS and RAYGAL means that the amplitudes of the total spectra including relativistic effects are close to each other. The amplitude is also in agreement with our MB estimate based on Born inverse magnification weight. Here again the relativistic effect can be understood following Eq. (16) where now the contribution of $C_{\delta\kappa}(l)$ becomes the dominant one compared to $C_{\delta\text{com},\delta\text{com}}(l)$ (and $C_{\kappa\kappa}(l)$). It means that the matter cross-spectra between two widely separated shell is a good estimate of (-2 times) the density-convergence spectrum.

We now investigate these two spectra ($C_{\kappa\kappa}(l)$ and $C_{\delta\kappa}(l)$) in detail in the next subsections.

Overall, we conclude that it is important to consider both RSD and weak-lensing effects including all the corrections from the linear to non-linear regime: our simulations provide a good ‘laboratory’ to test these effects by naturally including all these effects in a natural common framework.

### 3.3. Lensing: Convergence-convergence angular spectra

In this section we focus on the weak-lensing spectra, $C_{\kappa\kappa_2}(\ell, z_1, z_2)$, for shells 1 and 2 at redshifts $z_1 \in \{0.7, 1.8\}$ and $z_2 \in \{0.7, 1.8\}$. We follow exactly the same methodology as for the clustering. We remind that we focus on the convergence defined from Eq. (10).

#### 3.3.1. Born approximation

We investigate the gravitational convergence angular spectra under the Born approximation in Fig. 9. The spectra decrease with the multipole but now increase with redshifts (unlike the overdensity spectra). This is because the convergence is an integral of the overdensity (weighed by the lensing kernel) along the line of sight. A farther away source is therefore more distorted than a closer one. The cross-spectrum behaviour is similar to the auto-spectra one since $\kappa$ is not a local quantity but an integrated one: all the shells are therefore correlated. Here again we find a good agreement between CLASS prediction and RAYGAL measurements in the linear and non-linear regime. The right panel showing the relative deviation between the two indicates that they are compatible within the error bars between $\ell = 50$ and $\ell = 2000$. The deviation below $\ell = 50$ is related to the limited angular size of the cone. Between $\ell = 2000$ and $\ell = 5000$ there is a ~5–10% deviation at low redshift, the simulation measurement is slightly higher than the prediction from CLASS. We are not sure of the origin of the deviation, which might be related to shot noise. Despite such a discrepancy at large $\ell$ that seems to vanish at large redshift, this provides a cross-validation of RAYGAL Born convergence maps and the CLASS predictions.
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Measurements from the RAYGAL 2500 deg² light cone (diamonds) and the CLASS predictions (dashed lines). **Top left:** matter-density angular auto-spectrum with relativistic corrections (MB+RSDs) in a single shell at *z* = 0.7 (blue) and *z* = 1.8 (orange) in ΛCDM cosmology. The opposite of the cross-spectrum between two shells at *z* = 0.7 and *z* = 1.8 is shown in purple. **Top right:** relative deviation from comoving matter-density angular auto-spectrum at *z* = 0.7 in ΛCDM cosmology due to relativistic effects. The effects of the dilution term of MB (i.e. *s* = 0) on the observed matter-density spectrum is shown in green, and the effect of the dilution term of MB plus RSDs is shown in dark blue. An estimate of the MB effect using an inverse magnification (|μ<sub>Born</sub>|⁻¹) weight is shown as light blue triangles (see text for details). **Bottom left:** same but for *z* = 1.8, with the MB effect in light green, MB+RSDs in orange, and the |μ<sub>Born</sub>|⁻¹ weight MB estimate in green. **Bottom right:** same but for the *z* = 0.7 – *z* = 1.8 cross-spectrum with MB in light blue, MB+RSDs in purple, and the |μ<sub>Born</sub>|⁻¹ weight MB estimate in pink. Overall, there is an impressive agreement with CLASS even though there is a combination of MB and RSD effects. For the cross-spectrum, relativistic effects largely dominate. These matter cross-spectra become sensitive probes of the density-convergence spectra.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Gravitational convergence angular auto-spectrum in a single shell at *z* = 0.7 (red) and *z* = 1.8 (pink) in ΛCDM cosmology assuming the Born approximation. The cross-spectrum between two shells at *z* = 0.7 and *z* = 1.8 is shown in green. **Left:** spectrum from the RAYGAL 2500 deg² light cone (diamonds with error bars) and spectrum from CLASS (dashed line). **Right:** relative deviation from the CLASS prediction. Overall, the measured power spectra are in agreement with CLASS predictions within the error bars except at large scales (small *l*) due to the finite area of the light cone. We also note a shift at the ~5–10% level at large *l* and low redshifts, which vanishes at larger redshifts.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Measurements from the RAYGAL 2500 deg² light cone (diamonds) and the CLASS predictions (dashed lines). Top left: gravitational convergence auto-spectrum $C_{\kappa\kappa}^l$ with relativistic corrections (MB and RSDs) in a single shell at $z = 0.7$ (blue) and $z = 1.8$ (orange) in ΛCDM cosmology. The cross-spectrum between two shells at $z = 0.7$ and $z = 1.8$ is shown in purple. Top right: relative deviation from gravitational convergence angular auto-spectrum (under the Born approximation) at $z = 0.7$ in ΛCDM cosmology due to non-trivial relativistic effects. The full effects of MB (with $s = 0$) are shown in green, while the effect of RSDs(+MB) is shown in blue. An estimate of the MB effect using an inverse magnification ($|\mu_{\text{Born}}|^{-1}$) weight is shown as light blue triangles (see text for details). Bottom left: Same but for $z = 1.8$ with MB in light green, RSDs(+MB) in orange, and the $|\mu_{\text{Born}}|^{-1}$ weight MB estimate in green. Bottom right: Same but for the $z = 0.7 - z = 1.8$ cross-spectrum with MB in light blue, RSDs(+MB) in purple, and the $|\mu_{\text{Born}}|^{-1}$ weight MB estimate in pink. MB effects are not accounted for in CLASS and cause ~10–30% deviations. This confirms that the dominant effect is the dilution term of MB related to source averaging (as opposed to angular averaging).</figcaption></figure>

### 3.3.2. Relativistic effects

We now move to Fig. 10, which includes relativistic effects. The upper left plot shows a similar trend as with the Born approximation. However, a careful inspection shows that the points have moved from above CLASS predictions to below. We now focus on the three other panels (upper right, bottom left and bottom right), which show the relative deviations from the Born convergence spectra due to relativistic effects. It is interesting to note that RSDs do not seem to play any role in both CLASS predictions and RAYGAL measurements. This is not trivial since redshift perturbations modify the apparent distribution of sources (and thus the observed densities). However, these changes do not affect the lensing itself since they are mostly radial. Instead, MB plays an important role on the convergence spectra, which is not predicted by CLASS (CLASS relies on Born approximation, it can only account for the effect of MB on the overdensity but not on the convergence). The first contribution one can think of is the post-Born effect. These contributions are at the percent level and play a role mostly beyond $\ell = 2000$ (see e.g., Petri et al. 2017). The second (and main) contribution comes from (the dilution term of) MB: magnified regions have fewer sources than de-magnified ones. The decrease in the convergence spectra (computed from the observed position of sources) is not

negligible at all since it grows from 5% to 25% between $\ell = 100$ and $\ell = 5000$. The decrease is also more important at larger redshifts. This is because the convergence is weighed by the density, which itself is de-magnified by lensing. At first order (see Eq. (15)), $\kappa_i = \kappa_i^{\text{Born}}(1 + (5s - 2)\kappa_i^{\text{Born}})$ (with $s = 0$ in our mass-limited sample), and as a consequence, the convergence two-point spectrum (see Eq. (A.5)) is approximated by

$$C_{\kappa_i \kappa_j}(l) \approx C_{\kappa_i^{\text{Born}} \kappa_j^{\text{Born}}}(l) - 2C_{\kappa_i \kappa_j}(l) - 2C_{\kappa_i^2 \kappa_j}(l) + 4C_{\kappa_i^2 \kappa_j^2}(l). \quad (17)$$

The convergence spectrum depends on some $\langle \kappa_i^2 \rangle$ terms that are not necessarily small, especially at high redshifts.

Alternatively, it is possible to emulate such an effect using the Born approximation and weighting the convergence by the inverse magnification for $s = 0$. This is what we call Born + Weight in Fig. 10. The trend is the same as with the catalogues, with small discrepancies at large $\ell > 1000$. These most likely come from the fact that the value of the magnification (computed with the Born approximation) is not exactly the same as the density ratio in a pixel between the lensed and un- lensed matter distribution. However, this is a cross-check of the large effect due to lensing on convergence angular spectra.

Such a lensing bias has been investigated by Schmidt et al. (2009), who focused on the shear power spectrum and found
similar trend (in their notation, $q = 2$ corresponds to the oppo-
site of our $s = 0$ correction). They found the correction to
increase with $\ell$ and with redshift. For a source at $z = 1$ the
(norm of the) correction increase from 2% at $\ell = 500$ to 7%
at $\ell = 5000$ in their article, while we find a correction about
four times larger (at $z = 0.7$). More recently the same trend was
found but with smaller amplitude in Deshpande et al. (2020)$^{12}$.
Our result is not in contradiction with these findings since the
effect of MB is much stronger on the convergence than on the
shear (see Appendix B). The origin of the lensing convergence
bias is the same as the lensing shear bias: the magnified regions
are less sampled than the de-magnified ones leading to a bias.
However, the magnification is directly related at first order to
$\kappa$ (and not to $\gamma$) and the effect on the convergence power spec-
trum is therefore larger. The bi-spectrum $\langle\gamma\gamma\kappa\kappa\rangle$ is indeed much
smaller than the bi-spectrum $\langle\kappa\kappa\kappa\kappa\rangle$. An important consequence
is therefore that the convergence power spectrum (or correlation
function) can differ from the shear power spectrum (or correla-
tion function) by up to 10–30% due to MB.

To conclude, we find that we reproduce basic weak-lensing
results under the Born approximation and we capture non-trivial
effects, such as MB, which turn out to play an important role in
the resulting spectra.

3.4. Galaxy-galaxy lensing: Density-convergence angular cross-spectra

In this section we focus on the galaxy-galaxy lensing spectra, $C_{\delta_1\kappa_2}(\ell, z_1, z_2)$, for shells 1 and 2 at redshift ($z_1, z_2 \in \{0.7, 1.8\}$)$^2$. We follow exactly the same methodology as for the clustering and the weak-lensing parts. Investigating galaxy-galaxy lensing in the non-linear regime is challenging since it involves non-trivial relativistic effects on both the overdensity and the convergence. To our knowledge an accurate study of such effects down to the non-linear regime is still lacking and we now try to fill in the gap.

3.4.1. Comoving density and Born approximation

The cross-spectra between the gravitational convergence (computed with Born approximation) and the (comoving) matter overdensity are presented in Fig. 11. Unlike for lensing and clustering there are now four cross-spectra since $C_{\delta_1\kappa_2}(\ell, z_1 = 0.7, z_2 = 1.8)$ is not equal to $C_{\delta_1\kappa_2}(\ell, z_1 = 1.8, z_2 = 0.7)$. The dominant spectrum is the one with the overdensity at low redshift, which acts as a lens for the convergence at high redshift, $C_{\delta_1\kappa_2}(\ell, z_1 = 0.7, z_2 = 1.8)$. This spectrum is in good agreement with the CLASS expectation. The cross-spectra with the density and the convergence at the same redshift $C_{\delta_1\kappa_2}(\ell, z_1 = z_2)$ are smaller but non-zero. This is because our Born calculation accounts for the shell thickness. As a consequence, the density in the inner part of the shell acts as a lens for the convergence in the outer part of the shell. This contribution vanishes in the limit of infinitesimal shells (Dirac radial selection function). While the signal is clearly visible for the shell at $z = 0.7$, it becomes smaller and very noisy for the shell at $z = 1.8$ because the shell thickness is much smaller than the distance to the observer. Finally, the non-trivial spectrum correlating the convergence at $z = 0.7$ with the density at $z = 1.8$, $C_{\delta_1\kappa_2}(\ell, z_1 = 1.8, z_2 = 0.7)$, is nearly zero both in RAYGAL and CLASS expectations. The right panel shows the relative deviation from CLASS predictions.

Again, we observe a good agreement between the analytical pre-
diction and the numerical data points from $\ell = 50$ to $\ell = 5000$
for the main spectrum $C_{\delta_1\kappa_2}(\ell, z_1 = 0.7, z_2 = 1.8)$. The devi-
ations at lower $\ell$ are related to sample variance. The spectrum
$C_{\delta_1\kappa_2}(\ell, z_1 = 0.7, z_2 = 0.7)$ is noisy but compatible with CLASS
within the error bars. The spectra $C_{\delta_1\kappa_2}(\ell, z_1 = 1.8, z_2 = 1.8)$
seems roughly in agreement with the expectation between $\ell =
100$ and $\ell = 5000$ but the oscillations are extremely large and can
reach almost 100%. Finally, the spectrum $C_{\delta_1\kappa_2}(\ell, z_1 = 1.8, z_2 =
0.7)$ is outside the range of the graph because of a division by
nearly zero. Overall, there is a good agreement between CLASS
predictions and RAYGAL measurements.

3.4.2. Relativistic effects

We include relativistic effects in Fig. 12. The upper left panel shows that the 4 cross-spectra are now different from zero and they are even of the same order of magnitude (i.e. within a factor of ten in absolute value). Moreover, the order has changed since the second largest spectrum is now $C_{\delta_1\kappa_2}(\ell, z_1 = 1.8, z_2 = 1.8)$, which has been seriously boosted by relativistic effects. The most striking spectrum is $C_{\delta_1\kappa_2}(\ell, z_1 = 1.8, z_2 = 0.7)$, which is now of the same order as $C_{\delta_1\kappa_2}(\ell, z_1 = 0.7, z_2 = 0.7)$ although with a minus sign. Even though the changes are important, all the RAYGAL spectra seem at first sight in agreement with CLASS predictions (which include MB and RSDs) within a few $\sigma$ as we will see later. This is a good cross-check of the results. We now investigate the amplitude of relativistic effects to better understand the origin of the changes. The relativistic contributions in galaxy-galaxy lensing are presented in Ghosh et al. (2018) (using analytical predictions from CLASS). In their work and in all our panels, the dominant effect is the MB effect, which modifies the density (and, in a subtler way, the convergence, as shown previously). While RSDs are negligible here, they cannot always be neglected since they play a (small) role at low redshifts and low multipoles for some cross-spectra within thin shells (e.g., we checked that they are non-negligible but noisy contributions to $C_{\delta_1\kappa_2}(\ell, z_1 = 0.45, z_2 = 0.45)$ of order 10% for $\ell < 80$). The dominant contribution is again MB, where the apparent density becomes (at first order) $\delta_i = \delta_i^{\text{com}} + (5s-2)\kappa_i$ and the convergence becomes $\kappa_i = \kappa_i^{\text{Born}}(1+(5s-2)\kappa_i^{\text{Born}})$ (with $s=0$ here). As a consequence, there is an additional lensing contribution in the cross-spectrum (see Eq. (A.6)):

$$
C_{\delta_1 \kappa_j}(l) \approx C_{\delta_i^{\text{com}} \kappa_i^{\text{com}}}(l) - 2C_{\kappa_i \kappa_j}(l) - 2C_{\delta_1 \kappa_2^2}(l) + 4C_{\kappa_i \kappa_2^2}(l). \quad (18)
$$

The first two terms (as well as the RSD contribution) are included in CLASS but not the two last ones (nor the post-Born effects).

The relativistic contribution to $C_{\delta_1\kappa_2}(\ell, z_1 = 0.7, z_2 = 1.8)$ is presented in the bottom-right panel. The dominant relativistic contribution is MB with -15% between $\ell = 50$ and $\ell = 500$ and -20% between $\ell = 500$ and $\ell = 5000$. CLASS predicts a negative effect of order 10% beyond $\ell = 50$. Given the size of the error bars, the deviation is of order $3\sigma$ below $\ell = 500$ and $5\sigma$ above. The discrepancy is likely to be related to the effect of MB on the convergence (i.e. the two last terms of Eq. (18)), which also decrease the $\kappa - \kappa$ cross-spectra (as seen in Sect. 3.3) and is not accounted for in CLASS. This is confirmed by our MB estimate from Born inverse magnification weight. In this configuration, the relativistic contributions to galaxy-galaxy lensing cannot be neglected. The other galaxy-galaxy lensing spectra are dominated by relativistic effects. The upper right panel shows the relativistic contributions to $C_{\delta_1\kappa_2}(\ell, z_1 = 0.7, z_2 = 0.7)$. This negative effect ranges from -50% to -100%. The shape is highly

<footnote n="12">12 At first order, all the results can be compared by re-normalising the correction so that $5s - 2 \rightarrow -2$.</footnote>
<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Comoving matter density–gravitational convergence angular cross-spectrum in a single shell at $z = 0.7$ (red) and $z = 1.8$ (pink) in ΛCDM cosmology assuming the Born approximation. The cross-spectrum between two shells at $z = 0.7$ (density) and $z = 1.8$ (convergence) is shown in green, and the non-trivial cross-spectrum between two shells at $z = 0.7$ (convergence) and $z = 1.8$ (density) is shown in black. **Left:** spectrum from the RAYGAL 2500 deg$^2$ light cone (diamonds with error bars) and the spectrum from CLASS (dashed line). **Right:** relative deviation from the CLASS prediction. The cross-spectrum for a single shell at $z = 1.8$ being very noisy, the number of bins has been reduced and the large error bars have been omitted (for readability). The non-trivial cross-spectrum is nearly zero in real space, so the data points are outside the plot. Overall, the measured power spectra are in agreement with CLASS predictions within the error bars, except at large scales (small $l$) due to the finite area of the light cone.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Measurements from the RAYGAL 2500 deg$^2$ light cone (diamonds) and the CLASS predictions (dashed lines). **Top left:** matter density–gravitational convergence cross-spectra with relativistic corrections (MB+RSDs) in a single shell at $z = 0.7$ (blue) and $z = 1.8$ (orange) in ΛCDM cosmology. The (usual) cross-spectrum between two shells at $z = 0.7$ (density) and $z = 1.8$ (convergence) is shown in purple, and the (non-trivial) reverse one (cross-spectra between convergence at $z = 0.7$ and density at $z = 1.8$) is shown in green. **Top right:** relative deviation from comoving matter density–Born gravitational convergence cross-spectrum at $z = 0.7$ in ΛCDM cosmology due to relativistic effects. MB is shown in green, while the effect of RSDs(+MB) is shown in blue. An estimate of the MB effect using an inverse magnification ($|\mu_{\text{Born}}|^{-1}$) weight is shown as light blue triangles (see text for details) **Bottom left:** same but for $z = 1.8$ with MB in green, RSDs(+MB) in orange, and the $|\mu_{\text{Born}}|^{-1}$ weight MB estimate in green. **Bottom right:** same but for the $z = 0.7$ (density)-$z = 1.8$ (convergence) cross-spectrum with MB in light blue, RSDs(+MB) in purple, and the $|\mu_{\text{Born}}|^{-1}$ weight MB estimate in pink. Relativistic corrections are between 50 and 100% for the cross-spectrum at a single redshift. Overall, there is a good agreement with CLASS at the 20% level. Here again the overestimation of the power spectrum is related to the dilution term of the MB effect on the convergence.</figcaption></figure>
<figure><img src="image_14.png" /></figure>

**Fig. 13.** Non-trivial reverse configuration of gravitational convergence at $z = 0.7$ and matter density at $z = 1.8$, similar to the bottom-right panel of Fig. 12. MB is in grey, MB+RSDs in cyan, and the $|\mu_{\text{Born}}|^{-1}$ weight MB estimate in light green. Relativistic effects almost reach 100%, in agreement with CLASS. This configuration turns out to be a sensitive probe of the lensing convergence spectrum.

non-trivial with a bump near $\ell = 100$, a minimum near $\ell = 1000$ and another bump near $\ell = 5000$. All these effects are well reproduced by CLASS (with deviations at the $2\sigma$ level) and our $|\mu^{-1}|$ weight MB spectra below $\ell = 1000$. The bottom-left panel shows the contributions to $C_{\delta_1\kappa_2}(\ell, z_1 = 1.8, z_2 = 1.8)$, which are about +100% since $C_{\delta_{\text{com}}\kappa_{\text{Born}}}(\ell, z_1 = 1.8, z_2 = 1.8)$ was small because of the small thickness of the shell. The relativistic contribution agrees with the CLASS prediction. Finally, as expected, the cross-spectrum $C_{\delta_1\kappa_2}(\ell, z_1 = 1.8, z_2 = 0.7)$ is 100% explained by relativistic effects, and is compatible with the CLASS prediction (see Fig. 13). This non-trivial contribution can be understood using Eq. (18). It is dominated by $-2C_{\kappa_1\kappa_2}(\ell, z_1 = 1.8, z_2 = 0.7)$. The corresponding $\kappa-\kappa$ cross-spectrum presented in Fig. 10 is indeed very similar (when multiplied by a factor -2). For thin shells or when the density shell is beyond the convergence shell, it is therefore possible to probe lensing spectra from galaxy-galaxy lensing spectra. To conclude, the relativistic contributions to galaxy-galaxy lensing play an important role that cannot be neglected. The contributions of the MB effect (on the density) are well captured by CLASS while we found additional contributions from the MB (on the convergence). The actual impact of the MB depends on the slope $s$, and the shell configurations, but it can give a 20–100% effect.

## 3.5. Evolution and cosmological dependence

### 3.5.1. Redshift evolution

We now investigate the redshift evolution of relativistic effects. We therefore perform the same analysis with two shells at low redshifts $z_1 = 0.225 \pm 0.045$ and $z_2 = 0.45 \pm 0.025$ extracted from the full-sky catalogue. The resolution of the HEALPIX maps built from the catalogue is $N_{\text{side}} = 1024$ while the resolution of the Born map is $N_{\text{side}} = 2048$. These maps are then degraded to $N_{\text{side}} = 1024$. The estimation of the power spectrum is straightforward since the cone is full-sky (i.e. no need of apodisation to account for large masked region). We find that the global trend is similar to the one of high redshift measurements. There is also a good agreement with CLASS in most cases, while the deviations are well understood.

We now highlight the main differences between the high- and low-redshift cross-spectra – we illustrate some of them in

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

**Fig. 14.** Relativistic effects at low redshift on matter power spectrum $C_{\delta_1\delta_2}(\ell, z_1 = 0.225, z_2 = 0.225)$ (symbols are RAYGAL measurements, and lines are CLASS predictions). The MB effect (green) and RSDs(+MB) effect (blue) are shown. The trends are similar compared to higher redshifts, but the RSD effect plays a dominant role and the MB effects are smaller. Finger-of-God effects (ignored by CLASS) are also present.

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

**Fig. 15.** Relativistic effects at low redshift on convergence power spectrum $C_{\kappa_1\kappa_2}(\ell, z_1 = 0.45, z_2 = 0.45)$ (symbols are RAYGAL measurements, and lines are CLASS predictions). MB effects are in green and RSDs(+MB) in orange. The trends are similar compared to at higher redshifts, but MB effects are smaller.

Figs. 14–16. It shows the impact of relativistic effects on three spectra. Figure 14 shows the impact of relativistic effects on $C_{\delta_1\delta_2}(\ell, z_1 = 0.225, z_2 = 0.225)$. We notice that the influence of the MB is now negligible, while RSDs play a more important role. This is due to the weak-lensing effect, which decreases with redshift. CLASS performs a good prediction at large scales; however, beyond $\ell \sim 100$ (corresponding roughly to $k \sim 0.15 \, \text{h Mpc}^{-1}$), the effects of RSDs on the power spectrum become negative. As expected, this effect is not well reproduced by CLASS because of the assumption of linear mapping between real-space and redshift-space. Analytical predictions could be made using perturbation theory. A recent work by Grasshorn Gebhardt & Jeong (2020) about non-linear RSDs in the harmonic space galaxy power spectrum finds similar trend thus reinforcing our claim of the damping of the angular matter power spectrum at large $\ell$ due to the fingers-of-God effect.

We now turn to the lensing power spectrum. We verified that the agreement between the power spectrum of Born convergence map and CLASS expectation is better than 5%. The impact of relativistic effects on $C_{\kappa_1\kappa_2}(\ell, z_1 = 0.45, z_2 = 0.45)$ is illustrated in Fig. 15. The RSD effects are still negligible. The effect of MB
<figure><img src="image_17.png" /><figcaption>**Fig. 16.** Relativistic effects at low redshift on density-convergence spectrum $C_{\delta_1\kappa_3}(\ell, z_1 = 0.225, z_2 = 0.45)$ (symbols are RAYGAL measurements, and lines are CLASS predictions). The MB effect is in blue and RSDs(+MB) in purple. The trends are similar compared to at higher redshifts, but the MB effects are smaller.</figcaption></figure>

is also much smaller, being of order -2% over a wide range of multipoles from $\ell \sim 10$ to $\ell \sim 400$. We also noted that the MB effect remains large ($5$ to $15\%$) for lensing cross-spectra.

We then investigated galaxy-galaxy lensing at low redshifts. We checked that trends are similar and that the spectra are noisier (and some sign changes occurred, such as for $C_{\delta_1\kappa_2}(\ell, z_1 = 0.45, z_2 = 0.45)$). The impact of relativistic effects on the overdensity-convergence power spectrum $C_{\delta_1\kappa_2}(\ell, z_1 = 0.225, z_2 = 0.45)$ is illustrated in Fig. 16. The MB effect is negative of order -2% to -5% between $\ell = 100$ and $\ell = 1000$. The decrease is smaller than at higher redshift but still larger than the CLASS expectation (which ignores the MB effect on the convergence). Though RSD effects are small, we noted that, according to CLASS and our (noisy) measurement, they contribute to ~10% of $C_{\delta_1\kappa_2}(\ell, z_1 = 0.45, z_2 = 0.45)$ for $\ell < 80$. While RSD seem negligible in most galaxy-galaxy lensing studies, they might play a role in some configurations: it is therefore important to consider all relativistic effects and to verify their relative amplitude before neglecting some of these effects. To conclude, at lower redshift the relative importance of RSDs increases with respect to those of MB effects. It is also interesting to note that fingers-of-God effects are also visible within the angular power spectra.

### 3.5.2. Cosmological dependence

A detailed study of the full cosmological dependence of all cross-spectra and all relativistic effects is beyond the scope of this paper. Here we would like to broadly address the questions of whether the wCDM cross-spectra are valid; what the relative variation in the ten cross-spectra (including all relativistic effects) is when changing cosmology from ΛCDM to wCDM; if CLASS can predict this relative variation in the non-linear regime; and whether these spectra allow us to discriminate between these two competitive scenarios, both of which are compatible with CMB constraints at the $2\sigma$ level.

The analysis of the wCDM simulation is similar to the one of the ΛCDM. We consider two shells at the same redshift and we build maps with the same resolution. The ΛCDM and wCDM cosmologies do not differ strongly since they are both calibrated on CMB data. As a consequence the relativistic effects behave qualitatively in a similar way and we therefore did not plot again the same figures as in ΛCDM cosmology.

<figure><img src="image_18.png" /><figcaption>**Fig. 17.** Relative deviations of the matter-density spectra between wCDM and ΛCDM spectra (including all relativistic effects) considering one shell at $z = 0.7$ (blue), one shell at $z = 1.8$ (orange), and two shells at $z = 0.7$ and $z = 1.8$. Diamonds are measurements from the RAYGAL 2500 deg$^2$ light cone, and dashed lines are CLASS predictions. This plot shows the interest of considering all cross-spectra since they are more or less sensitive to cosmology. Even though the two models were calibrated on CMB data, the differences can reach 20%. The agreement with CLASS is good.</figcaption></figure>

In Figs. 17–19 we investigate the relative variation in the ten density-convergence auto- and cross-spectra (including relativistic effects) when moving from ΛCDM to wCDM. The first global remark is that there is a good agreement below $\ell = 5000$ between our results and CLASS calculation, which is an encouraging validation of our work. In Fig. 17, the matter-density auto-spectrum in wCDM cosmology differs from the one of ΛCDM cosmology by ~-10% near $\ell \sim 30$, ~0% near $\ell \sim 250$ and ~+10% near $\ell \sim 3000$. The same trend can be found in the 3D matter power spectrum at $z = 1$ in Fig. 6. The relative variation in the cross-spectrum (which is dominated by relativistic effect) reach larger values of order +20% near $\ell \sim 5000$. All the dependence (including the non-trivial cross-spectrum) are well reproduced by CLASS (since the 3D power spectra were calibrated to RAYGAL simulations). In Fig. 18, the convergence power spectrum shows a larger relative difference of 15–25% at $\ell \sim 30, 30–40\%$ at $\ell \sim 1000$ and 30% at $\ell \sim 5000$. There is a good agreement with CLASS expectation up to $\ell \sim 500$. Beyond this multipole, the relative variation is overestimated by CLASS. This is likely to be related to the MB effect, which is not included in the analytical prediction. The difference remains reasonable of order 10%. In Fig. 19, the density-convergence cross-spectra are presented. The cosmological dependence of the stronger signal $C_{\delta_1\kappa_2}(\ell, z_1 = 0.7, z_2 = 1.8)$ (with smaller error bars in purple) rises from 0% near $\ell \sim 40$ up to ~20% at $\ell \sim 1000$ and beyond. The agreement with CLASS is good at the $2-\sigma$ level. The spectrum $C_{\delta_1\kappa_2}(\ell, z_1 = 0.7, z_2 = 0.7)$ shows a similar trend rising from 0% near $\ell \sim 100$ up to ~20% at $\ell \sim 2000$ and beyond. Finally, the spectra $C_{\delta_1\kappa_2}(\ell, z_1 = 1.8, z_2 = 1.8)$ and $C_{\delta_1\kappa_2}(\ell, z_1 = 1.8, z_2 = 0.7)$ dominated by relativistic effects show a larger variation rising from 20% at $\ell \sim 30$ to 40% near $\ell \sim 1000$ and beyond. It indicates that even though these spectra are noisier, they carry a strong cosmological signal. The cosmological dependence is well captured by CLASS for overdensity-convergence spectra.

To conclude, if one assumes that the cross-spectra are accurately known in ΛCDM, the cosmological dependence for the models in the vicinity of ΛCDM can be captured at first order by CLASS. To go further the modelling of MB and non-linear RSD
<figure><img src="image_19.png" /><figcaption>**Fig. 18.** Relative deviations of the convergence spectra between wCDM and ΛCDM spectra (including all relativistic effects) considering one shell at *z* = 0.7 (blue), one shell at *z* = 1.8 (orange), and two shells at *z* = 0.7 and *z* = 1.8. Diamonds are measurements from the RAYGAL 2500 deg² light cone, and dashed lines are CLASS predictions. This plot shows the interest of considering all cross-spectra since they are more or less sensitive to cosmology. Even though the two models were calibrated on CMB data, the differences can reach 40%. The agreement with CLASS is good except for *ℓ* > 500, which shows a discrepancy of ~10%.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Fig. 19.** Relative deviations of the density-convergence spectra between wCDM and ΛCDM spectra (including all relativistic effects) considering one shell at *z* = 0.7 (blue), one shell at *z* = 1.8 (orange), and two shells at *z* = 0.7 and *z* = 1.8 (purple for the standard configuration and green for the reverse one). Diamonds are measurements from RAYGAL 2500 deg² light cone, and dashed lines are CLASS predictions. This plot shows the interest of considering all cross-spectra since they are more or less sensitive to cosmology. Even though the two models were calibrated on CMB data, the differences can reach 40%. The agreement with CLASS is good.</figcaption></figure>

should be improved. In any case, the cross-spectra carry complementary cosmological information (the cosmological densities are different). For the studied wCDM and ΛCDM cosmologies, the weak-lensing power spectra and the galaxy-galaxy spectra dominated by relativistic effect seems the most sensitive probes followed by the relativistic matter power spectrum and the standard density-convergence spectrum. Relativistic 3×2 points spectra are therefore a promising cosmological probe.

## 4. Conclusion

In this work we have introduced a simulation suite for the study of relativistic effects in cosmology and illustrated its use with an application to lensing-matter clustering statistics.

We have described the RAYGAL suite. It consists of two large and well-resolved N-body simulations of two different cosmological models (ΛCDM and wCDM) calibrated on CMB data. For each simulation, several snapshots have been stored in addition to (wide, intermediate, or deep) particles and gravity light cones. Halos have been detected in snapshots and light cones according to a wide range of possible definitions. The specificity of this work lies in the relativistic ray-tracing technique used to generate the data: light rays are propagated from the observer to the sources according to weak-field geodesic equations of GR. Several hundred million particles and tens of millions of halos have been replaced in redshift space. Within our intermediate light cone (2500 deg² up to *z* = 2), the source (or halo) angular densities are similar to those expected from present and future photometric (or spectroscopic) surveys, such as DES, KIDS, HSC, *Euclid*, and LSST (or BOSS, eBOSS, DESI, *Euclid*, and SKA). These relativistic particle and halo catalogues contain information about the lensing deflection angle, the weak-lensing distortion matrix, and the apparent redshift, including all contributions at first order in metric perturbation (comoving, Doppler, potential, transverse Doppler, and ISW-RS). Relativistic HEALPIX maps were also produced. All these catalogues and maps are publicly available<fnref n="13" />.

The catalogues and maps are very generic: they allow a wide range of applications. We illustrate the strength of the approach with an application to the study of lensing-matter clustering statistics, which is a hot topic in modern cosmology. Since the matter overdensity and the gravitational convergence are the fundamental quantities of clustering and weak lensing, we focus on the position-position, position-convergence, and convergence-convergence 3×2 points correlations. The extension to more observable quantities, such as the usual 3×2 points involving the shear instead of the convergence (see e.g., DES Collaboration 2022), can be investigated in a future work. We have considered two shells extracted from our intermediate light cone: a low-redshift shell ($z_1$ = 0.7 ± 0.2) and a high-redshift shell ($z_2$ = 1.8 ± 0.1). We computed the matter overdensity, $\delta$, and gravitational convergence, $k$. We then evaluated the ten auto- and cross-spectra between overdensity and convergence at these two redshifts (including the non-trivial ones) up to $\ell$ = 5000. Finally, we investigated the relative amplitude of all the contributions beyond the standard one (evaluated from comoving overdensity, $\delta^{\text{com}}$, and Born convergence, $k^{\text{Born}}$): MB and RSDs. These contributions were compared to CLASS predictions. Our main results are the following:

- **Validation:** There is an agreement at the few percent level between the standard cross-spectra (without relativistic effects) and CLASS predictions once the 3D real-space matter power spectra are properly calibrated to the spectrum of RAYGAL simulations (which itself is in agreement at the percent level with that of COSMICEMU).

- **Clustering:** Both RSDs and MB contribute to angular density spectra. Magnification bias is a nearly scale-independent contribution at the 2–5% level, while the RSD contribution is small for $\ell$ > 100–500 (depending on redshift) but can reach 20–40% at small $\ell$ < 20–100. Magnification bias also correlates distant shells (~100% contribution to the cross-spectra). In this configuration, the matter cross-power spec-

<footnote n="13">¹³ https://cosmo.obspm.fr/public-datasets/</footnote>
trum becomes a sensitive probe of the density-convergence spectrum. The agreement with CLASS predictions (which include both MB and RSDs) is excellent at these redshifts, but not at lower ones.

- **Weak lensing:** The dominant effect is MB, whose contribution rises from ~5–10% to ~30% between $\ell = 200$ and $\ell = 5000$. Magnification bias is not included in CLASS. Since the MB contribution to the shear power spectrum is small, it indicates that the gravitational convergence and shear power spectrum can differ by ~20–30% at large $\ell$. The reason for the lensing bias is that highly magnified regions are less sampled than de-magnified ones. However, the magnification is directly linked at first order to the convergence. The convergence is more correlated to itself than to the shear, and thus the convergence lensing bias is much larger than the shear lensing bias. This bias can have observational implications for the interpretation of surveys that probe sources of known size or known luminosity.

- **Galaxy-galaxy lensing:** The dominant non-trivial contribution to overdensity-convergence cross-spectra is MB, which plays a role in both the overdensity and the convergence. The contribution is of order -20% in $C_{\delta_1\kappa_2}(\ell, z_1 = 0.7, z_2 = 1.8)$, ~ -50% in $C_{\delta_1\kappa_2}(\ell, z_1 = 0.7, z_2 = 0.7)$, and ~100% in $C_{\delta_1\kappa_2}(\ell, z_1 = 1.8, z_2 = 1.8)$ and $C_{\delta_1\kappa_2}(\ell, z_1 = 1.8, z_2 = 0.7)$. The three last cross-spectra are usually neglected in cosmological analysis (and in particular the last one), but their amplitude is not negligible compared to $C_{\delta_1\kappa_2}(\ell, z_1 = 0.7, z_2 = 1.8)$. For thin shells or when the density shell is behind the convergence shell, the density-convergence cross-spectrum becomes a sensitive probe of the gravitational convergence power spectrum. CLASS captures only the effect of MB on the overdensity. Redshift-space distortions are usually negligible but may play a role at the ~10% level at low $\ell$ for some specific configurations.

- **Redshift evolution:** At low redshift, $z < 0.5$, we observe the same trends, although MB plays a smaller role compared to RSDs. Even though we consider angular spectra, fingers-of-God effects are visible at large $\ell$ as a damping of the matter power spectra. These effects are not included in CLASS.

- **Cosmological dependence:** By changing the cosmology to wCDM, we observe a relative variation in the ten cross-spectra. However, the amplitudes of the variations are different for each spectrum. For the specific model we considered, the largest variations are observed for convergence spectra and for the density-convergence spectra dominated by relativistic effects, followed by the relativistic matter-density cross-spectrum and the standard density-convergence spectrum. This indicates that relativistic 3×2 points spectra are powerful cosmological probes.

- **Correlation function:** There is a good agreement between the convergence angular correlation function estimated from convergence maps weighted by the observed density (i.e. the inverse transform of the spectra presented in this article) and the direct pair-based estimate from catalogues. The agreement is not observed without density weight nor when removing the monopole and dipole. However, MB plays an important role in all cases (including with a pixel-based estimate without a density weight).

The publicly available catalogues are very generic (they have not been built for a specific community, but instead from first principle). They may help when investigating any N-point cross-correlation between several observables (each being sensitive to a particular combination of relativistic effects). We identify some examples of applications of RAYGAL data, but the reader can

probably find many others. We also point to some illustrative articles on these topics that have already made use of RAYGAL data:

- **Galaxies (weak-lensing studies):** This work can be extended to shear and ellipticities, more realistic selection functions, and higher-order statistics.

- **Galaxies (clustering and RSD studies):** The impact of relativistic effects on higher-order statistics, the impact on baryon acoustic oscillations, and two-points correlation functions (e.g., Breton et al. 2019; Di Dio & Seljak 2019; Taruya et al. 2020; Di Dio & Beutler 2020; Saga et al. 2020) can be studied.

- **Galaxy cluster, group, and void studies:** The profile (e.g., Corasaniti et al. 2018), mass function, gravitational redshift, weak lensing, and relativistic cross-correlation.

- **Supernovae, standard candles, and standard rulers:** The fluctuation of the cosmic distances (e.g., Breton & Fleury 2021).

- **CMB studies:** The ISW-RS effect, weak lensing, and CMB cross-correlation.

- **Other sources:** Relativistic effects in the Lyman-$\alpha$ forest, in intensity mapping, and so on.

In the future, we plan to expand our source catalogue (to include galaxies and active galactic nuclei), to explore alternative cosmological models (modified gravity, etc.), to extend our work to strong lensing (i.e. by implementing a multiple root finder since our current version finds one image in the case of strong lensing), to consider alternative messengers (gravitational waves), and to investigate the role of baryonic physics. Relativistic ray tracing is a promising tool of modern cosmology: it provides a natural framework for combining multiple cosmological probes and may pave the way to shedding light on the nature of the dark sector.

**Acknowledgements.** We acknowledge financial support from the DIM ACAV of the Région Île-de-France, the Action fédératrice Cosmologie et structuration de l'univers as well as the JSPS Grant L16519. This work was granted access to HPC resources of TGCC/CINES through allocations made by GENCI (Grand Équipement National de Calcul Intensif) under the allocations 2016-042287, 2017-A0010402287, 2018-A0030402287, 2019-A0050402287 and 2020-A0070402287. A.T. acknowledges the support from MEXT/JSPS KAKENHI Grant Nos. JP17H06359, JP20H05861, JP21H01081, and JST AIP Acceleration Research Grant No. JP20317829, Japan. We thank E. Di Dio for pointing out the “hidden” CLASS options to tune, S. Prunet for help on MPGRAFIC and POLSPICE, R. Teyssier for help on RAMSES, B. Li for sharing his TSC routine, C. Murray for useful comments on the catalogues, J. Adamek for tips about HEALPIX. We also thank J.-M. Alimi, S. Colombi, I. Achitouv and A. Le Brun for fruitful discussions. We deeply thank S. Mene for tremendous efforts in making the local servers working at LUTH. S.A. was supported in part by the project “Combining Cosmic Microwave Background and Large Scale Structure data: An Integrated Approach for Addressing Fundamental Questions in Cosmology”, funded by the MIUR Progetti di Rilevante Interesse Nazionale (PRIN) Bando 2017, Grant No. 2017YJYZAH.

## References

- Adamek, J., Daverio, D., Durrer, R., & Kunz, M. 2016a, Nat. Phys., 12, 346

- Adamek, J., Daverio, D., Durrer, R., & Kunz, M. 2016b, JCAP, 7, 053

- Adamek, J., Clarkson, C., Coates, L., Durrer, R., & Kunz, M. 2019, Phys. Rev. D, 100, 021301

- Alestas, G., Kazantzidis, L., & Perivolaropoulos, L. 2020, Phys. Rev. D, 101, 123516

- Alimi, J. M., Bouillot, V., Rasera, Y., et al. 2012, in Proc. Int. Conf. on High Performance Computing, Networking, Storage and Analysis, 73, 11

- Benitez-Llambay, A. 2015, py-sphviewer: Py-SPHViewer v1.0.0.

- Bernardeau, F., Bonvin, C., & Vernizzi, F. 2010, Phys. Rev. D, 81, 083002

- Beutler, F., & Di Dio, E. 2020, JCAP, 2020, 048

- Bonvin, C. 2008, Phys. Rev. D, 78, 123530

- Bonvin, C., & Durrer, R. 2011, Phys. Rev. D, 84, 063505
Bonvin, C., Clarkson, C., Durrer, R., Maartens, R., & Umeh, O. 2015, JCAP, 2015, 040

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

Bouillot, V. R., Alimi, J.-M., Corasaniti, P.-S., & Rasera, Y. 2015, MNRAS, 450, 145

Breton, M.-A., & Fleury, P. 2021, A&A, 655, A54

Breton, M. A., & Reverdy, V. 2022, A&A, in press https://doi.org/10.1051/0004-6361/202142661

Breton, M.-A., Rasera, Y., Taruya, A., Lacombe, O., & Saga, S. 2019, MNRAS, 483, 2671

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

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

Chisari, N. E., Richardson, M. L. A., Devriendt, J., et al. 2018, MNRAS, 480, 3962

Chon, G., Challinor, A., Prunet, S., Hivon, E., & Szapudi, I. 2004, MNRAS, 350, 914

Coates, L., Adamek, J., Bull, P., Guandalin, C., & Clarkson, C. 2021, MNRAS, 504, 3534

Corasaniti, P. S., Ettori, S., Rasera, Y., et al. 2018, ApJ, 862, 40

Crocce, M., Pueblas, S., & Scoccimarro, R. 2006, MNRAS, 373, 369

DES Collaboration (Abbott, T. M. C., et al.) 2022, Phys. Rev. D, 105, 023520

Deshpande, A. C., Kitching, T. D., Cardone, V. F., et al. 2020, A&A, 636, A95

Di Dio, E., & Beutler, F. 2020, JCAP, 2020, 058

Di Dio, E., & Seljak, U. 2019, JCAP, 2019, 050

Di Dio, E., Montanari, F., Lesgourgues, J., & Durrer, R. 2013, JCAP, 2013, 044

Fidler, C., Rampf, C., Tram, T., et al. 2015, Phys. Rev. D, 92, 123517

Fidler, C., Tram, T., Rampf, C., et al. 2016, JCAP, 9, 031

Fleury, P., Larena, J., & Uzan, J.-P. 2017, Phys. Rev. Lett., 119, 191101

Fleury, P., Larena, J., & Uzan, J.-P. 2019, Phys. Rev. D, 99, 023525

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

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

Francfort, J., Cusin, G., & Durrer, R. 2021, CQG, 38, 245008

Ghosh, B., Durrer, R., & Sellentin, E. 2018, JCAP, 2018, 008

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

Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759

Grasshorn Gebhardt, H. S., & Jeong, D. 2020, Phys. Rev. D, 102, 083521

Guandalin, C., Adamek, J., Bull, P., et al. 2021, MNRAS, 501, 2547

Guillet, T., & Teyssier, R. 2011, J. Comput. Phys., 230, 4756

Heitmann, K., Lawrence, E., Kwan, J., Habib, S., & Higdon, D. 2014, ApJ, 780, 111

Heitmann, K., Bingham, D., Lawrence, E., et al. 2016, ApJ, 820, 108

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

Kilbinger, M. 2015, Rep. Prog. Phys., 78, 086901

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

Komatsu, E., Smith, K. M., Dunkley, J., et al. 2011, ApJS, 192, 18

Lepori, F., Adamek, J., Durrer, R., Clarkson, C., & Coates, L. 2020, MNRAS, 497, 2078

Lesgourgues, J. 2011, ArXiv e-prints [arXiv:1104.2932]

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

Michaux, M., Hahn, O., Rampf, C., & Angulo, R. E. 2021, MNRAS, 500, 663

Petri, A., Haiman, Z., & May, M. 2017, Phys. Rev. D, 95, 123503

Planck Collaboration VI. 2020, A&A, 641, A6

Prunet, S., Pichon, C., Aubert, D., et al. 2008, ApJS, 178, 179

Rasera, Y., Corasaniti, P.-S., Alimi, J.-M., et al. 2014, MNRAS, 440, 1420

Reverdy, V. 2014, PhD Thesis, Laboratoire Univers et Théories, https://hal.archives-ouvertes.fr/tel-02095297/document, France

Riess, A. G., Casertano, S., Yuan, W., Macri, L. M., & Scolnic, D. 2019, ApJ, 876, 85

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

Saga, S., Taruya, A., Breton, M.-A., & Rasera, Y. 2020, MNRAS, 498, 981

Schmidt, F., Rozo, E., Dodelson, S., Hui, L., & Sheldon, E. 2009, ApJ, 702, 593

Schneider, P., van Waerbeke, L., Kilbinger, M., & Mellier, Y. 2002, A&A, 396, 1

Sgier, R., Fluri, J., Herbel, J., et al. 2021, JCAP, 2021, 047

Szapudi, I., Prunet, S., Pogosyan, D., Szalay, A. S., & Bond, J. R. 2001, ApJ, 548, L115

Takahashi, R., Oguri, M., Sato, M., & Hamana, T. 2011, ApJ, 742, 15

Takahashi, R., Sato, M., Nishimichi, T., Taruya, A., & Oguri, M. 2012, ApJ, 761, 152

Taruya, A., Saga, S., Breton, M.-A., Rasera, Y., & Fujita, T. 2020, MNRAS, 491, 4162

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

Teyssier, R., Pires, S., Prunet, S., et al. 2009, A&A, 497, 335

Turner, E. L., Ostriker, J. P., & Gott, J. R. I. 1984, ApJ, 284, 1

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

Yoo, J., Fitzpatrick, A. L., & Zaldarriaga, M. 2009, Phys. Rev. D, 80, 083514
## Appendix A: Power-spectrum expressions under the linear mapping assumption

Under the linear mapping assumption (which is only used for analytical predictions), the matter overdensity $\delta_i$ of a shell $i$ can be decomposed as

$$ \delta_i \approx \delta_i^{\text{com}} + \delta_i^{\text{MB}} + \delta_i^{\text{RSD}}, \quad (A.1) $$

where $\delta_i^{\text{com}}$ is the comoving overdensity, $\delta_i^{\text{MB}} = \frac{1+\delta_i^{\text{com}}}{\mu_i} - (1 + \delta_i^{\text{com}}) \approx -2\kappa_i$ is the MB contribution dominated by the dilution term ($\mu_i$ is the magnification and we assumed a slope $s=0$), and $\delta_i^{\text{RSD}}$ is the contribution from redshift perturbations (see Breton et al. 2019 and reference therein). The matter-density auto- and cross-power spectra are given by

$$ \begin{aligned} C_{\delta_i \delta_j}(l) &\approx C_{\delta_i^{\text{com}} \delta_j^{\text{com}}}(l) \\ &\quad + C_{\delta_i^{\text{Perm}} \delta_j^{\text{MB}}}(l) + C_{\delta_i^{\text{MB}} \delta_j^{\text{MB}}}(l) \\ &\quad + C_{\delta_i^{\text{Perm}} \delta_j^{\text{RSD}}}(l) + C_{\delta_i^{\text{MB}} \delta_j^{\text{RSD}}}(l) + C_{\delta_i^{\text{RSD}} \delta_j^{\text{RSD}}}(l), \end{aligned} \quad (A.2) $$

where $C_{\delta_i^{\text{Perm}} \delta_j^{\text{MB}}}(l) = C_{\delta_i^{\text{MB}} \delta_j^{\text{MB}}}(l) + C_{\delta_i^{\text{MB}} \delta_j^{\text{MB}}}(l)$. The first line is the comoving power spectrum. The second line adds the weak-lensing contribution (i.e. MB with $s=0$). The third line adds the RSD contributions.

The convergence can be approximated at first order as

$$ \kappa_i \approx \kappa_i^{\text{Born}}, \quad (A.3) $$

where $\kappa_i^{\text{Born}}$ is the convergence according to Born approximation. However, in cosmological analysis, one needs to perform some averaging procedure to extract some statistics. While the above expression is recovered with direction averaging, it is not with source averaging (which is used in the analysis of observational data, in the calculation of power spectra and correlation functions). The average convergence depends on the apparent overdensity and is therefore sensitive to both MB and RSDs. The average becomes $\langle \kappa_i \rangle = (\Sigma_{\text{sources}}(1+\delta_i)\kappa_i)(\Sigma_{\text{sources}}(1+\delta_i))^{-1}$ where $\Sigma_{\text{sources}}$ is the sum over the sources. As a consequence, the (source) average convergence can be written as

$$ \langle \kappa_i \rangle_s \approx \kappa_i^{\text{Born}} + \kappa_i^{\text{MB}} + \kappa_i^{\text{RSD}}, \quad (A.4) $$

where $\kappa_i^{\text{MB}} \approx \frac{\kappa_i^{\text{Born}}}{\mu_i} - \kappa_i^{\text{Born}} + \kappa_i^{\text{PB}} \approx -2(\kappa_i^{\text{Born}})^2 + \kappa_i^{\text{PB}}$ is the MB contribution (assuming $s=0$; see Schmidt et al. 2009) and $\kappa_i^{\text{RSD}} \approx (\Sigma_{\text{sources}}(1+\delta_i^{\text{com}} + \delta_i^{\text{RSD}})\kappa_i^{\text{Born}})(\Sigma_{\text{sources}}(1+\delta_i^{\text{com}} + \delta_i^{\text{RSD}}))^{-1} - \kappa_i^{\text{Born}}$ is the contribution from redshift perturbation. The MB contribution includes the sub-dominant post-Born correction $\kappa_i^{\text{PB}}$, which plays a percent level role at large $\ell > 2000$. It is dominated by a geometric contribution (i.e. the integral of second derivative of the potential along the true light-ray path) and a lens-lens coupling contribution (see Petri et al. 2017).

In the following we omit the average to make the equations more readable. The resulting power spectrum for the gravitational convergence can be expressed as

$$ \begin{aligned} C_{\kappa_i \kappa_j}(l) &\approx C_{\kappa_i^{\text{Born}} \kappa_j^{\text{Born}}}(l) \\ &\quad + C_{\kappa_i^{\text{Perm}} \kappa_j^{\text{MB}}}(l) + C_{\kappa_i^{\text{MB}} \kappa_j^{\text{MB}}}(l) \\ &\quad + C_{\kappa_i^{\text{Perm}} \kappa_j^{\text{RSD}}}(l) + C_{\kappa_i^{\text{MB}} \kappa_j^{\text{RSD}}}(l) + C_{\kappa_i^{\text{RSD}} \kappa_j^{\text{RSD}}}(l), \end{aligned} \quad (A.5) $$

where the first line is the convergence spectrum under the Born approximation, the second line adds MB (assuming $s=0$) and the third line adds the RSD contribution.

Given our decomposition of the overdensity and convergence, we can proceed similarly to evaluate the galaxy-galaxy lensing power spectrum,

$$ \begin{aligned} C_{\delta_i \kappa_j}(l) &\approx C_{\delta_i^{\text{com}} \kappa_j^{\text{Born}}}(l) \\ &\quad + C_{\delta_i^{\text{com}} \kappa_j^{\text{MB}}}(l) + C_{\delta_i^{\text{MB}} \kappa_j^{\text{Born}}}(l) + C_{\delta_i^{\text{MB}} \kappa_j^{\text{MB}}}(l) \\ &\quad + [C_{\delta_i^{\text{RSD}} \kappa_j^{\text{Born}}}(l) + C_{\delta_i^{\text{RSD}} \kappa_j^{\text{MB}}}(l) + C_{\delta_i^{\text{RSD}} \kappa_j^{\text{RSD}}}(l) + \\ &\quad C_{\delta_i^{\text{com}} \kappa_j^{\text{RSD}}}(l) + C_{\delta_i^{\text{MB}} \kappa_j^{\text{RSD}}}(l)], \end{aligned} \quad (A.6) $$

where the first line is the density-convergence spectrum under the comoving and Born approximation, the second line adds MB (assuming $s=0$) and the third line and fourth line (within the square bracket) adds the RSD contribution.

In CLASS, all the weak-field relativistic effects and wide angle effects are included. All the terms are evaluated from the knowledge of the 3D matter power spectrum. The approximations we identified are the following: the code makes use of the above linear mapping assumption (which is an approximation in particular for $\delta_i^{\text{RSD}}$), the matter power spectrum is evaluated from HALOFIT, the Born approximation is used, the Doppler contribution to $\delta_i^{\text{RSD}}$ is evaluated approximately and the terms $\kappa_i^{\text{MB}}$, $\kappa_i^{\text{RSD}}$ are not included.

## Appendix B: Convergence two-point correlation function

In this section we perform configuration-space measurements for the convergence two-point correlation function (the harmonic-space measurements can be found in Sect. 3.3). This allows us to test the methodology used in this paper, based on the HEALPIX maps and angular power spectrum estimators.

To do so, we used the code ATHENA (Schneider et al. 2002). First, we verified that the configuration-space measurement on the HEALPIX map containing the convergence using the Born approximation gave the same results as its counterpart using POLSPICE.

In Fig. B.1, we show the relative difference between the measurements made on the catalogue containing the dilution bias with respect to Born at $z = 1.8$. We see that the measurements from ATHENA and POLSPICE agree extremely well, and yield a ~20% difference with respect to Born at the arcminute scale and ~50% at 100 arcmin. This confirms that the dilution bias has a strong impact on the convergence two-point statistics as in Fig. 10, where we used the same POLSPICE measurement. We note that had we used POLSPICE and subtract the average and dipole to the data (which is common practice), then the results would vastly differ, especially at larger angular scales. This most likely comes from the fact that when averaging over sources (which differs from averaging over directions; Kibble & Lieu 2005; Bonvin et al. 2015; Breton & Fleury 2021), $\langle \kappa \rangle = -2\langle \kappa^2 \rangle \neq 0$. This effect, in the form of $\langle \kappa \rangle \langle \kappa \rangle$ gains a relative importance as the signal decreases with angular scale, hence our results. Moreover we expect some leakage from small $\ell$ to large and intermediate angles. It is interesting to see that there is already a ~10% discrepancy at 10 arcmin.

Finally, we show the relative difference when computing the cosmic shear. We see that in this case, the results are very close to the Born approximation. This shows that the very large effect of lensing bias mostly impacts the convergence, and only to a lesser the extent cosmic shear. This comes from the fact that the source-averaged convergence is non-zero (while the averaged shear is, neglecting intrinsic alignment) and the convergence three-point
<figure><img src="image_21.png" /></figure>

**Fig. B.1.** Relative difference between lensing angular two point correlation function on the source catalogue accounting for the dilution bias and the Born convergence angular two point correlation function. In red and green diamonds we show the measurements of cosmic shear and convergence correlation function using ATHENA, and in blue and light blue we show the results using the same methodology as in Sect. 3, keeping and removing the monopole and dipole, respectively.

correlation function contribution is large (i.e. $\langle\kappa\kappa\rangle/\langle\kappa\kappa\rangle$ is larger than $\langle\gamma\gamma\kappa\rangle/\langle\gamma\gamma\rangle$).