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

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. 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. RayGalGroupSims (or in short RayGal) consists of two N-body simulations of $(2625\,h^{-1}\,{\rm Mpc})^3$ volume with $4096^3$ particles of a standard flat $\Lambda$CDM model and a non-standard $w$CDM phantom dark energy model. Light-cone data from the simulations have been generated using a parallel ray-tracing algorithm that has accurately solved billion geodesic equations. Catalogues and maps with relativistic weak-lensing which include post-Born effects, magnification bias (MB) and redshift space distortions (RSD) due to gravitational redshift, Doppler, transverse Doppler, 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 10 (cross-)power spectra (3$\times$2 points) of the matter density fluctuation field and the gravitational convergence at $z=0.7$ and $z=1.8$. We find $1-30\%$ level contribution from both MB and RSD to the matter power spectrum, while the Fingers-of-God effect is visible at lower redshift in the non-linear regime. MB 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. MB also plays a significant role in the galaxy-galaxy lensing by decreasing the density-convergence spectra by $20\%$, while coupling non-trivial configurations (such as the one with the convergence at the same or even lower redshift than the density field).


Introduction
Observations of the distribution of galaxies in the Universe carry a wealth of cosmological information pertaining to the cosmic Article number, page 1 of 25 arXiv:2111.08745v2 [astro-ph.CO] 28 Jul 2023 matter content, the state of expansion, the amplitude of matterdensity 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 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 twopoints 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;. 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(Fidler et al. , 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 , 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 Raytracing Galaxy Group Simulations, hereafter RayGal) simulation suite 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 nonnegligible 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 Model Ω m σ 8 w ΛCDM 0.25733 0.80101 -1.0 wCDM 0.27508 0.85205 -1.2 Notes. The columns from left to right give the values of the cosmic matter density, Ω m , the linear root-mean-square fluctuations smoothed on the 8 h −1 Mpc scale, σ 8 , and the redshift-independent equation of state, w. In all models the baryon density is set to Ω b = 0.04356, the density of the relativistic species to Ω r ≃ 0.00008, the curvature to Ω k = 0, the reduced Hubble constant to h = 0.72, and the slope of the primordial spectrum to n S = 0.963.
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; we present a case study of the imprint of relativistic effects on weaklensing observables using the RayGal data in Sect. 3; and finally we discuss the conclusions in Sect. 4.

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 Ω m and σ 8 of the wCDM model have been chosen so as to be within the 1σ confidence region of the WMAP inferred constraints in the Ω m − σ 8 plane and along the degeneracy line of σ 8 − w contours (see beige filled contours in Fig. 1). Henceforth, this model is statistically indistinguishable at 1σ 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σ (TT, TE, and EE, the temperature and polarisation spectra Planck Collaboration et al. 2020) (see red-yellow filled contours in Fig. 1). We notice that values of the equation of state varying in the range −1.3 ≲ w ≲ −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.

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 h −1 Mpc) 3 sampled by 4096 3 dark-matter particles. The optimal compromise between mass resolution (of order 0.7 0.8 0.9 1.0 1. ∼ 2 × 10 10 h −1 M ⊙ depending on the simulated cosmology) and large statistics makes them ideal to investigate the cosmologydependent distribution and properties of halos from Milky Way 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. ∼ 0.64 h −1 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).
Simulations were performed with (an optimised version of) the parallel AMR N-body code RAMSES 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 max = 0.5, 2, and 10, respectively) and apertures (full sky, 2500, and 400 deg 2 , respectively) using the onion-shells techniques of Fosalba et al. (2008);Teyssier et al. (2009). Unlike many other studies, the maximum redshift z 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φ, 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 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 ∆ m = 50 to ∆ 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.

Relativistic ray tracing and catalogue creation
To produce realistic observables, we use the ray-tracing framework Magrathea-pathfinder (Breton & Reverdy 2021) 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 with ϕ the gravitational potential, η 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 ν k ν = 0, where k α = (dη/dλ, dx i /dλ) and λ 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 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: (2) with ϕ the gravitational potential and a ′ = da/dη the time derivative of the scale factor. In particular, we use the precise value of ∂ϕ/∂η 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.
6.87 · 10 10 4096 7 8 4.41(4.46) · 10 11 1.88(2.01) · 10 10 5 46(51) 0.5 CAMB Notes. L box is the box length, n part the number of particles, n x the grid size, n ref the number of refinements, m ref the refinement threshold, n cell the final number of AMR cells, m p the particle mass, ∆ x the spatial resolution, z i the starting redshift, C dt the Courant-like factor that determines the size of the integration time step (with respect to the largest stable time step), and P lin (k) specifies the source of the linear power spectrum used to generate initial conditions. All calculations have been performed in double precision. 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 methodology described in Breton et al. (2019), which consists of the following 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 milliarcsecond). 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-A&A proofs: manuscript no. aanda cludes 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 position of the source, θ, is known, as is the comoving angular position, β, we can compute the deflection angle and evaluate the lensing distortion matrix, A, along the true photon trajectory. This reads as where κ and γ = (γ 1 , γ 2 ) are the weak-lensing quantities related, respectively, to the convergence and shear, while ω 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 ∼ 10 h −1 kpc located at a comoving distance of ∼ 1 h −1 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 ob-servable 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 straightforward. 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. However, 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 observable 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 spectrum 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 proposals using the polarisation of radio-galaxies (Francfort et al. 2021). Here we provide both the fundamental lensing quantities 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 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 ≈ 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 sourceaveraging, 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(Fleury et al. , 2019Breton & Fleury 2021). Finally, even if strong lensing events occur, our geodesic finder would most likely find the principal image of the source.

Maps
Maps were built by projecting the matter distribution and its properties onto a sphere. We discretised the maps by consider-ing 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 mean and (half-)width ∆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.

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 Ray-Gal data are available on the website of the COS team at LUTH 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 (κ, γ 1 , γ 2 , and 1/µ, respectively the convergence, shear components, and inverse magnification).
Data | z z i 3.00 2.33 2.00 1.50 1.00 0.67 0.50 0.43 0.25 0.11 0.00 other z Part. x Notes. We note that halo properties are available for about 30 to 40 'other z' per cosmology, but the redshifts may differ from one cosmology to another. z is the redshift, z i the starting redshift, 'Part.' stands for particles, and 'prop.' stands for properties. Notes. 'Charac.' stands for characteristics, z max 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.
We provide Dirac radial selection maps for several redshifts within the range of our light cones. Top-hat radial selection maps (for κ, γ 1 , γ 2 and δ) 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 Table 3 and Table 4. An illustration of the apparent matter distribution and the various relativistic effects within the light cone 'Narrow 2' are shown in Fig. 2 and Fig. 3 6 . In this cone, the source angular density (∼ 40 arcmin −2 ) 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 −3 h 3 Mpc −3 ) 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, δ) and that of lensing (i.e. the gravitational convergence, κ) 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. 6 We made use of Py-SPHViewer (Benitez-Llambay 2015) for the realisation of these images.

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 × N 2 side 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.
We considered two shells (1 and 2) extracted from the ΛCDM 2500 deg 2 light cone. We isolated two possible redshifts for each shell: z = 0.7 ± 0.2 and z = 1.8 ± 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).

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 (κ Born ). The mean Born convergence for 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. a population of source galaxies (see Kilbinger 2015) is given by with χ the comoving distance (in a flat universe), χ S the comoving distance of the source, χ lim the limiting comoving distance of the galaxy sample, n(χ S ) the source probability distribution, and ∆ ⊥ ϕ 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 shell and constant density, n(χ S ) ∝ 4πχ 2 S /V shell . It is worth noting that a common additional approximation is to assume with H 0 the Hubble constant. We do not make this assumption in our Born calculation. Moreover in this expression δ is often estimated within shells of typical size between 10 − 100 h −1 Mpc.
Here we estimate ∆ ⊥ ϕ 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 δ com . The overdensity is estimated as the local excess of particles (in real space) with respect to the mean matter density, where N com (θ) 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 galaxygalaxy lensing studies (see Ghosh et al. 2018).
We therefore call 'relativistic effects' all the non-trivial GR effects beyond these two references (δ com and κ 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 + δ = (1 + δ com )µ −1 µ 2.5s with µ the magnification and s the logarithmic slope of the luminosity function. The first term µ −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 A&A proofs: manuscript no. aanda  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. 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 Fig. 4 and Fig. 5.
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) 8 . The density auto-spectra are polluted by shot noise at large ℓ. 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). 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.

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 where δ com i is the comoving overdensity, δ MB i ≈ −2κ i is the MB contribution dominated by the dilution term (κ 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 δ RSD i is the contribution from redshift perturbations. In the context of source averaging (as in galaxy surveys) one can define an average convergence 9 as where κ Born i is the convergence according to Born approximation, κ MB i ≈ −2(κ Born i ) 2 is the MB contribution (assuming s = 0) and κ RSD i is the contribution from redshift perturbation. More details as well as the expression of the angular spectra under the linear mapping assumption are provided in Appendix A.
By including relativistic effects in our simulations, we follow a strategy similar to that of studies of the CMB. The two most widely used CMB codes in the literature are Camb (Lewis et al. 2000) and Class (Lesgourgues 2011), which provide predictions for these effects (Challinor & Lewis 2011;Di Dio et al. 2013). Here, we confront our results to analytical predictions from the latter code, since by adjusting the name list 10 it automatically computes the ten cross-spectra we are interested in.
Class 11 is a linearised Einstein-Boltzmann solver, very accurate in the linear regime. It predicts lensing, density, and densitylensing angular power spectra without relying on Limber approximation. Lensing predictions are performed at first order (Born without MB) and without RSDs. However, density and density-lensing predictions include MB and RSDs (Doppler effect, gravitational redshift, ISW-RS). Class can also make predictions in the non-linear regime and we use them as our reference predictions. In the non-linear regime it relies on several assumptions, such as: (i) weak lensing under the Born approximation (MB cannot be accounted for, and post-Born effects are neglected), (ii) RSDs with the linear mapping assumption from real space to redshift space and approximate treatment of the Doppler effect in the non-linear regime (i.e. no finger-of-God effect), and (iii) the Halofit prescription (Takahashi et al. 2012) for the non-linear matter power spectrum P(k).
Using the Halofit prescription induces up to 5% errors in the predicted matter power spectrum around k = 1 h Mpc −1 (Heitmann et al. 2014). The question of the error of P(k) due to the Halofit prescription is beyond the scope of this work, which is focused on relativistic effects. We therefore prefer to feed Class with a matter power spectrum interpolated from the RayGal suite. The precision of the fit becomes ∼ 1% between k = 0.002 h Mpc −1 and k = 20 h Mpc −1 (beyond which we still use Halofit prescription) for redshift between z = 0 and z = 2.33. Below k = 0.2 h Mpc −1 , the effect of cosmic variance was reduced by subtracting the difference between the linear spectrum in the simulation (i.e. a specific realisation of initial conditions) and the linear spectrum from the Camb code (i.e. ensemble average). We notice that at large wavenumbers we did not remove the shot noise because its effect was found negligible in the relevant range of scales for the current study (ℓ < 5000).
As we will see in the following, one of the main effects that cannot be captured by Class (beyond non-linear RSDs at low redshift) is the effect of MB on convergence-convergence spectra and density-convergence spectra since it involves the bispectrum ⟨κκκ⟩. To cross-check the validity of our results (which is obtained from an accurate but complicated ray-tracing procedure) we computed from our Born maps, the absolute value of the inverse magnification |µ −1 Born |. From this we computed a (rough) overdensity map including MB and a (rough) convergence map including MB by weighting the comoving density and the Born convergence by |µ −1 Born |. More precisely the approximate density and convergence maps are given by δ weight MB = |µ −1 Born |(1 + δ com ) − 1 and κ weight MB = |µ −1 Born |κ Born . We then computed the '|µ −1 Born | weight MB spectra' by cross-correlating these maps. They are not analytical predictions but rather an estimate of the MB effect directly from the Born maps allowing for further cross-check of our ray-tracing results.
To conclude, we computed the ten auto-and cross-spectra for two redshift shells both analytically and using the RayGal data. For each case we computed a reference spectrum (comoving and/or Born), a lensed spectrum accounting for MB correction and a relativistic spectrum accounting for both MB corrections and relativistic redshift perturbations (MB+RSDs). Additionally, we computed a rough estimate of MB effect from Born inverse magnification maps: this |µ −1 Born | weight MB spectrum serves as a cross-validation of our MB results when analytical predictions are not available.

Results: Analysis of density-lensing angular power spectra and cross-spectra (3×2 points)
Here, we present the results of the computation of the angular matter-density and gravitational convergence auto-and crosspower spectra, C δ,δ ; C δ,κ ; C κ,κ , from the RayGal data.  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 Λ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 ΛCDM (continuous lines) and wCDM (dashed lines) cosmologies. Percent level agreement is shown up to k = 1 − 2 h Mpc −1 , where resolution effects damp the power spectrum. Baryonic effects are also expected to damp the power spectrum near this wavenumber.

Validation of the 3D Cartesian real-space matter power spectrum
Although the real-space distribution is not the target of this study, we perform an analysis of the 3D matter power spectrum (computed from a snapshot of periodic-boundary-conditions simulations at a fixed redshift) as a test of robustness of the numerical analysis. First, we determine the range of validity of the matter power spectrum. Second, the evaluated matter power spectrum is used as input for our analytical predictions of relativistic effects.
We compute the 3D matter power spectrum in real space at z = 0 and z = 1 (for both the ΛCDM and wCDM runs) with the fast Fourier transform based parallel estimator Powergrid (Prunet et al. 2008). The density is computed using a Cloud-In-Cell assignment scheme within a cartesian grid with 8192 3 elements (the power spectrum is then deconvolved from the assignment window). The comparison with spectra from the emulator Cosmicemu (Heitmann et al. 2016) shows a 1% agreement up to k = 1-2 h Mpc −1 for both cosmologies (see Fig. 6). We notice that the location of the drop in power, due to finite mass resolution effect, is similar to that induced by baryonic physics and in particular active galactic nuclei feedback (Chisari et al. 2018) in cosmological simulations with hydrodynamics. It indicates that a higher-resolution simulation would not necessarily improve matter power spectra predictions unless (costly) baryonic physics is properly and accurately taken into account (which is a nontrivial task). The increase in power near k = 7 − 10 h Mpc −1 is due to both shot noise and aliasing, which have not been subtracted. The oscillations at small wavenumbers are mostly related to sample variance: they are limited to the 2 − 3% level down to k = 0.02 h Mpc −1 where linear theory starts to reach percent level precision. The percent level agreement between Cosmicemu and RayGal spectra over 2 decades of wavenumbers (k = 0.02 to 2 h Mpc −1 ) is a good cross-check of the precision of our real-space data.

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 ≤ ℓ ≤ 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 autospectra C δ 1 δ 2 (ℓ, z 1 = 0.7, z 2 = 0.7) and C δ 1 δ 2 (ℓ, 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 (ℓ < 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 ℓ = 50 and ℓ = 500 validates our simulation measurement. The non-linear regime begins near ℓ ∼ 500, corresponding approximately to k ∼ 0.2 h Mpc −1 at z ∼ 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 ℓ = 500 and ℓ = 5000. We halt at ℓ = 5000, which corresponds approximately to k ∼ 2 h 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 ℓ we consider is of order of N side , we also expect the aliasing effect to be small. The green curve shows the cross-spectrum C δ 1 δ 2 (ℓ, 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 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,  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 crossspectrum 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 ℓ = 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 ℓ) due to the finite area of the light cone.
this is a good cross-validation of Class and RayGal's random subset matter angular spectra.

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 autospectra, 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 ℓ = 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 ℓ = 50 to ℓ = 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) δ i = δ com i + (5s − 2)κ 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 The RSD contributions to the multipoles of the 3D galaxygalaxy 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 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 ℓ ∼ 100. For the shell thickness we use, it can reach 5% at ℓ = 20 (for the z = 0.7 ± 0.2 spectra) and 20% at ℓ = 90 (for the z = 1.8 ± 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 ℓ ≈ 1000 using our MB estimate based on Born inverse magnification weight (we discuss the origin of the small discrepancy at large ℓ and low redshift in the lensing section). For both shells, the RSD effect is too small at large ℓ 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 fingersof-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 ℓ = 50 and ℓ = 5000). We note that we normalised the spectrum variation by the Class spectrum (even for RayGal spectrum). 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 (|µ Born | −1 ) 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 |µ 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, MB+RSDs in purple, and the |µ Born | −1 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.
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 δκ (l) becomes the dominant one compared to C δ com δ com (l) (and C κκ (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 κκ (l) and C δκ (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.

Lensing: Convergence-convergence angular spectra
In this section we focus on the weak-lensing spectra, C κ 1 κ 2 (ℓ, z 1 , z 2 ), for shells 1 and 2 at redshifts z 1 ∈ {0.7, 1.8} and z 2 ∈ {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).

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 autospectra one since κ 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 z1=0.700−z2=0.700 z1=1.800−z2=1.800 z1=0.700−z2=1.800 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 2 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 ℓ) due to the finite area of the light cone. We also note a shift at the ∼5-10% level at large ℓ and low redshifts, which vanishes at larger redshifts.  Top left: Gravitational convergence auto-spectrum C 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 (|µ 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 |µ 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 |µ 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).
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 ℓ = 50 and ℓ = 2000. The deviation below ℓ = 50 is related to the limited angular size of the cone. Between ℓ = 2000 and ℓ = 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 ℓ that seems to vanish at large redshift, this provides a cross-validation of RayGal Born convergence maps and the Class predictions.

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 ℓ = 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 ℓ = 100 and ℓ = 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)), κ i = κ Born i (1 + (5s − 2)κ Born i ) (with s = 0 in our mass-limited sample), and as a consequence, the convergence two-point spectrum (see Eq. (A.5)) is approximated by The convergence spectrum depends on some ⟨κ i κ 2 j ⟩ 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 ℓ > 1000. These most likely come from the fact that the value of the magnification (computed with he 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 opposite of our s = 0 correction). They found the correction to increase with ℓ and with redshift. For a source at z = 1 the (norm of the) correction increase from 2% at ℓ = 500 to 7% at ℓ = 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 κ (and not to γ) and the effect on the convergence power spectrum is therefore larger. The bi-spectrum ⟨γγκ⟩ is indeed much smaller than the bispectrum ⟨κκκ⟩. An important consequence is therefore that the convergence power spectrum (or correlation function) can differ from the shear power spectrum (or correlation 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.

Galaxy-galaxy lensing: Density-convergence angular cross-spectra
In this section we focus on the galaxy-galaxy lensing spectra, C δ 1 κ 2 (ℓ, z 1 , z 2 ), for shells 1 and 2 at redshift (z 1 , z 2 ) ∈ {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 nontrivial 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.

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 δ 1 κ 2 (ℓ, z 1 = 0.7, z 2 = 1.8) is not equal to C δ 1 κ 2 (ℓ, 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 δ 1 κ 2 (ℓ, 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 δ 1 κ 2 (ℓ, 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 δ 1 κ 2 (ℓ, 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 prediction and the numerical data points from ℓ = 50 to ℓ = 5000 for the main spectrum C δ 1 κ 2 (ℓ, z 1 = 0.7, z 2 = 1.8). The deviations at lower ℓ are related to sample z1=0.700−z2=1.800 z1=1.800−z2=0.700 (OUTSIDE) 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 ℓ) due to the finite area of the light cone.
variance. The spectrum C δ 1 κ 2 (ℓ, z 1 = 0.7, z 2 = 0.7) is noisy but compatible with Class within the error bars. The spectra C δ 1 κ 2 (ℓ, z 1 = 1.8, z 2 = 1.8) seems roughly in agreement with the expectation between ℓ = 100 and ℓ = 5000 but the oscillations are extremely large and can reach almost 100%. Finally, the spectrum C δ 1 κ 2 (ℓ, 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.

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 δ 1 κ 2 (ℓ, z 1 = 1.8, z 2 = 1.8), which has been seriously boosted by relativistic effects. The most striking spectrum is C δ 1 κ 2 (ℓ, z 1 = 1.8, z 2 = 0.7), which is now of the same order as C δ 1 κ 2 (ℓ, 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 σ 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 δ 1 κ 2 (ℓ, z 1 = 0.45, z 2 = 0.45) of order 10% for ℓ < 80). The dominant contribution is again MB, where the apparent density becomes (at first order) δ i = δ com i + (5s − 2)κ i and the convergence becomes κ i = κ Born i (1 + (5s − 2)κ Born i ) (with s = 0 here). As a consequence, there is an additional lensing contribution in the cross-spectrum (see Eq. (A.6)): 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 δ 1 κ 2 (ℓ, z 1 = 0.7, z 2 = 1.8) is presented in the bottom-right panel. The dominant relativistic contribution is MB with −15% between ℓ = 50 and ℓ = 500 and −20% between ℓ = 500 and ℓ = 5000. Class predicts a negative effect of order 10% beyond ℓ = 50. Given the size of the error bars, the deviation is of order 3σ below ℓ = 500 and 5σ 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 κ − κ 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 δ 1 κ 2 (ℓ, z 1 = 0.7, z 2 = 0.7). This negative effect ranges from −50% to −100%. The shape is highly non-trivial with a bump near ℓ = 100, a minimum near ℓ = 1000 and another bump near ℓ = 5000. All these effects are well reproduced by Class (with deviations at the 2σ level ) and our |µ −1 |weight MB spectra below ℓ = 1000. The bottom-left panel shows the contributions to C δ 1 κ 2 (ℓ, z 1 = 1.8, z 2 = 1.8), which are about +100% since C δ com κ Born (ℓ, 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 δ 1 κ 2 (ℓ, 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 κ 1 κ 2 (ℓ, z 1 = 1.8, z 2 = 0.7). The corresponding κ − κ 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, Top left: Matter densitygravitational 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 (|µ 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 |µ 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 |µ 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. it is therefore possible to probe lensing spectra from galaxygalaxy 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.

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 ± 0.045 and z 2 = 0.45 ± 0.025 extracted from the full-sky catalogue. The resolution of the Healpix maps built from the catalogue is N side = 1024 while the resolution of the Born map is N side = 2048. These maps are then degraded to N 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 highand low-redshift cross-spectra -we illustrate some of them in Fig. 14, Fig. 15, and Fig. 16. It shows the impact of relativistic effects on three spectra. Fig. 14 shows the impact of relativistic effects on C δ 1 δ 2 (ℓ, 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 ℓ ∼ 100 (corresponding roughly to k ∼ 0.15 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  . 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 |µ 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.  Relativistic effects at low redshift on matter power spectrum C δ 1 δ 2 (ℓ, 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. thus reinforcing our claim of the damping of the angular matter power spectrum at large ℓ 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 κ 1 κ 2 (ℓ, z 1 = 0.45, z 2 = 0.45) is illustrated in Fig. 15. The RSD effects are still negligible. The effect of MB is also much smaller, being of order −2% over a wide range of multipoles from ℓ ∼ 10 to ℓ ∼ 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 δ 1 κ 2 (ℓ, z 1 = 0.45, z 2 = 0.45)). The impact of relativistic effects on the overdensity-convergence power spectrum C δ 1 κ 2 (ℓ, z 1 = 0.225, z 2 = 0.45) is illustrated in Fig. 16. The MB effect is negative of order −2% to −5% between ℓ = 100 and ℓ = 1000.   . 16. Relativistic effects at low redshift on density-convergence spectrum C δ 1 κ 2 (ℓ, 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.
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 δ 1 κ 2 (ℓ, z 1 = 0.45, z 2 = 0.45) for ℓ < 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.

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 l 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 Ray-Gal 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.   . 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 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 40%. The agreement with Class is good except for ℓ > 500, which shows a discrepancy of ∼ 10% . 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σ 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 z1=1.800−z2=0.700 z1=0.700−z2=1.800 z1=1.800−z2=1.800 z1=0.700−z2=0.700 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 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 40%. The agreement with Class is good. 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. In Fig. 17, Fig. 18 and Fig. 19 we investigate the relative variation in the ten densityconvergence 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 ℓ = 5000 between our results and Class calculation, which is an encouraging validation of our work. In Fig. 17, the matter-density autospectrum in wCDM cosmology differs from the one of ΛCDM cosmology by ∼ −10% near ℓ ∼ 30, ∼ 0% near ℓ ∼ 250 and ∼ +10% near ℓ ∼ 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 ℓ ∼ 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 ℓ ∼ 30, 30 − 40% at ℓ ∼ 1000 and 30% at ℓ ∼ 5000. There is a good agreement with Class expectation up to ℓ ∼ 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 δ 1 κ 2 (ℓ, z 1 = 0.7, z 2 = 1.8) (with smaller error bars in purple) rises from 0% near ℓ ∼ 40 up to ∼ 20% at ℓ ∼ 1000 and beyond. The agreement with Class is good at the 2-σ level. The spectrum C δ 1 κ 2 (ℓ, z 1 = 0.7, z 2 = 0.7) shows a similar trend rising from 0% near ℓ ∼ 100 up to ∼ 20% at ℓ ∼ 2000 and beyond. Finally, the spectra C δ 1 κ 2 (ℓ, z 1 = 1.8, z 2 = 1.8) and C δ 1 κ 2 (ℓ, z 1 = 1.8, z 2 = 0.7) dominated by relativistic effects show a larger variation rising from 20% at ℓ ∼ 30 to 40% near ℓ ∼ 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.
-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 crosscorrelation 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. 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. 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 funtion 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. ⟨κκκ⟩/⟨κκ⟩ is larger than ⟨γγκ⟩/⟨γγ⟩).