Observing relativistic features in large-scale structure surveys – II. Doppler magnification in an ensemble of relativistic simulations

: The standard cosmological model is inherently relativistic, and yet a wide range of cosmological observations can be predicted accurately from essentially Newtonian theory. This is not the case on ‘ultralarge’ distance scales, around the cosmic horizon size, however, where relativistic effects can no longer be neglected. In this paper, we present a novel suite of 53 fully relativistic simulations generated using the gevolution code, each covering the full sky out to z ฀ 0.85, and approximately 1930 deg2 out to z ฀ 3.55. These include a relativistic treatment of massive neutrinos, as well as the gravitational potential that can be used to exactly calculate observables on the past light cone. The simulations are divided into two sets, the first being a set of 39 simulations of the same fiducial cosmology (based on the Euclid Flagship 2 cosmology) with different realizations of the initial conditions, and the second that fixes the initial conditions, but varies each of seven cosmological parameters in turn. Taken together, these simulations allow us to perform statistical studies and calculate derivatives of any relativistic observable with respect to cosmological parameters. As an example application, we compute the cross-correlation between the Doppler magnification term in the convergence, ฀v, and the CDM + baryon density contrast, ฀cb, which arises only in a (special) relativistic treatment. We are able to accurately recover this term as predicted by relativistic perturbation theory, and study its sample variance and derivatives with respect to cosmological parameters. ABSTRACT The standard cosmological model is inherently relativistic, and yet a wide range of cosmological observations can be predicted accurately from essentially Newtonian theory. This is not the case on ‘ultra-large’ distance scales, around the cosmic horizon size, however, where relativistic eﬀects can no longer be neglected. In this paper, we present a novel suite of 53 fully relativistic simulations generated using the gevolution code, each covering the full sky out to 𝑧 ≈ 0 . 85 , and approximately 1930 square degrees out to 𝑧 ≈ 3 . 55 . These include a relativistic treatment of massive neutrinos, as well as the gravitational potential that can be used to exactly calculate observables on the past light cone. The simulations are divided into two sets, the ﬁrst being a set of 39 simulations of the same ﬁducial cosmology (based on the Euclid Flagship 2 cosmology) with diﬀerent realisations of the initial conditions, and the second which ﬁxes the initial conditions, but varies each of seven cosmological parameters in turn. Taken together, these simulations allow us to perform statistical studies and calculate derivatives of any relativistic observable with respect to cosmological parameters. As an example application, we compute the cross-correlation between the Doppler magniﬁcation term in the convergence, 𝜅 𝑣 , and the CDM+baryon density contrast, 𝛿 cb , which arises only in a (special) relativistic treatment. We are able to accurately recover this term as predicted by relativistic perturbation theory, and study its sample variance and derivatives with respect to cosmological parameters.


INTRODUCTION
The large-scale distribution of matter in our Universe will be studied at ever greater detail with upcoming astronomical surveys such as LSST (Abell et al. 2009), Euclid (Laureĳs et al. 2011), DESI (Levi et al. 2013), and the Roman Telescope (Spergel et al. 2015). The complicated astrophysical feedback processes that plague the dynamics of galaxies and clusters become subdominant at cosmological distance scales, making large-scale structure an ideal laboratory for probing gravity, which is effectively the only important force at those scales. This is usually done by extracting summary statistics such as power spectra, bispectra etc. (Leclercq et al. 2014;Bonvin & Durrer 2011;Yoo et al. 2009). As the volume of surveys increases, these will be measured more accurately and for an increasing number of modes. As a result, previously unconstrained small effects can be detected with high significance.
Cosmological -body simulations provide a powerful and versatile means to predict these summary statistics given a cosmological model. These simulations commonly use Newtonian theory (Springel 2005;Teyssier 2002) which is sufficient for many purposes, in particular in the context of the ΛCDM concordance model. At extremely large distance scales, the interpretation of such simulations becomes subtle, however (Rigopoulos & Valkenburg 2015;Green & Wald 2012;Chisari & Zaldarriaga 2011). One way of maintaining consistency with general relativity at leading order is ★ l.j.c.coates@qmul.ac.uk by using so-called Newtonian motion gauges (Fidler et al. 2017), but the full machinery for analysing simulations in this context still needs to be developed.
Alternatively, general relativity can be implemented in the simulations explicitly. Employing techniques from numerical relativity, this has been explored e.g. in Giblin et al. (2016) and Macpherson et al. (2017). The main drawback of this formulation is the requirement to keep track of the wave-like solutions of the gravitational field, which needs extremely fine time resolution and thus leads to practical limitations. In the weak-field regime relevant for cosmology, however, one can easily perform a scalar-vector-tensor decomposition of the gravitational field in order to isolate these wave-like components. They can then be treated with fast approximate methods that completely remove the limitation on the time stepping. Such an approach is implemented in the weak-field relativistic -body code gevolution (Adamek et al. 2016a,b) that we employ in this work. Relativistic effects that appear at extremely large distance scales are naturally included in our numerical simulations. Physical quantities, whenever they are gauge-dependent, are computed in the Poisson gauge, which is widely employed in practical calculations. Furthermore, using ray tracing we can compute physical observables exactly and directly.
While the final observables (i.e. the observed summary statistics) in actual surveys contain the fully aggregated information of how matter is distributed and observed on our past light cone, in the weak-field regime it can often be useful to study different effects separately. In the literature, one often finds that any effects that are not explained by weak lensing or specifically the Kaisertype redshift-space distortions (Kaiser 1987) are called relativistic effects, even though confusingly some of them would still be captured in Newtonian simulations, and even though weak lensing and redshift-space distortions are both arguably relativistic in nature as well. Some examples include corrections in the two-point correlation function, which can be as large as 10% (Lorenz et al. 2018;Tansella et al. 2018;Bertacca et al. 2012;Bonvin 2014;Beutler & Di Dio 2020;Raccanelli et al. 2016;Yoo & Desjacques 2013), or in bispectra where the signal-to-noise ratio for the relativistic part is ∼10 for a survey like Euclid De Weerd et al. 2020;Clarkson et al. 2019;Jolicoeur et al. 2018;Umeh et al. 2017;Bertacca et al. 2018). Such corrections are also relevant for the study of non-Gaussianity and bias (Bruni et al. 2012;Umeh et al. 2019;Wang et al. 2020;Alonso et al. 2015;Camera et al. 2015;Fonseca et al. 2015). The benefit of inherently relativistic simulations is that all such effects are transparently included, and so there is no ambiguity in the predictions for observable quantities.
In this paper, we present the UNITY simulations, a set of 53 fully relativistic -body simulations for which we have retained~35 TB of data in the form of HEALP maps of different fields as a function of comoving distance from a chosen observation point, and many power spectra etc. derived from these fields. A large portion of our simulations are run using the same fiducial cosmology, but using a varying random seed to generate the initial conditions, so as to give us multiple random realisations of the same underlying cosmology. The rest of the simulation suite contains pairs of simulations with the same random initial conditions, but each of the cosmological parameters varied by some percentage around the fiducial value. This allows us to compute numerical derivatives to determine the effect of each parameter on various summary statistics and observables. A full ray-tracing procedure can be applied to the stored data in order to produce light cones, or more selective treatments can be used to study particular fields and observables in isolation. The data are available on request 1 .
The layout of the paper is as follows. In Sect. 2, we describe the simulations in more detail, including the data products that are available. We then show some examples where the simulation data are used to extract a relativistic signal in Sect. 3, and then finally we conclude in Sect. 4.

SIMULATIONS
For our simulations, we use the -body code gevolution. The code is described in full in Adamek et al. (2016a), but we will give a brief overview of its workings here. gevolution employs the Friedmann-Lemaître-Robertson-Walker (FLRW) metric with perturbations in Poisson gauge and a weak-field setting where all gravitational fields ( , , , and ℎ ) are small. These metric quantities are stored on a Cartesian grid and are evolved together with the -body ensemble that describes the cold dark matter (CDM) in phase space. The joint evolution is therefore computed using a particle-mesh approach, keeping a fixed and uniform resolution on the mesh. The metric quantities can be saved on spatial hypersurfaces at any redshift alongside a full particle snapshot. The newest version of gevolution 2 also allows the data to be saved on a series of light cones given a specific (or several) observers. In this case, the metric quantities are saved on a series of HEALP 3 (Górski et al. 2005; Zonca et al. Table 1. Details of the UNITY simulations. The first row is based on the fiducial cosmology, with multiple different realisations. This gives the ability to study the covariance. The second row is also run using the fiducial cosmology, but gevolution is run in the -body gauge. These can be used to compare and determine gauge effects in our simulation. Finally we list a series of pairs of simulations where we vary one of the parameters around the fiducial value. This allows us to calculate finite differences on different statistics.

The UNITY suite of simulations
Our simulations have a box volume of(4032 Mpc/h) 3 , where all of the metric quantities are calculated on a Cartesian grid with 2304 3 grid points. This means that we have spatial and mass resolutions of 1.75 Mpc/h and ≈ 4.6 × 10 11 M ⊙ /ℎ, respectively. We use 2304 3 particles to sample cold dark matter and baryons, while the massive neutrino species are treated with a grid-based approach similar to the one described in Brandbyge & Hannestad (2009) that uses the linear transfer functions. For the case where we explore a nonstandard equation of state parameter of dark energy, 0 = −0.9 instead of 0 = −1, we use a very similar approach to account for the perturbations of the dark energy fluid, see Dakin et al. (2019) and, in particular, Hassani et al. (2019) for details. As these perturbations are often neglected in Newtonian simulations of evolving dark energy, we also run a simulation for comparison where the dark energy is perfectly homogeneous (in Poisson gauge).
Each simulation has two light cones with the same observer at the corner of the box. The first light cone extends out to a distance of 2015 Mpc/h and covers the full sky 4 occupying a corner of the box., whereas the second one, illustrated in Fig. 1, extends further out to a distance of 4690 Mpc/h with a disk-shaped survey area of approximately 1932 square degrees.The pointing at the centre of the disk is towards the opposite corner of the box to allow us the maximum distance without repeating data, which could add spurious correlations (especially on large scales).
The UNITY simulations consist of a total of 53 simulations, the parameters of which are summarised in Table 1. Of these, 39 use the same fiducial cosmology of = 0.96, = 2.1 × 10 −9 , ℎ = 0.67, = 0.021996, = 0.121203, = 0.06 eV and 0 = −1. These parameters match the fiducial cosmology of the Euclid Flagship 2 simulations . For each of these simulations we vary the random seed, which gives us a different realisation of initial conditions for the same cosmology.
For the remaining simulations we vary several of the parameters individually, in turn, adding or subtracting a small amount from the fiducial value in such a way that we can form finite difference derivatives with respect to the parameters at a later stage, where is an observable quantity of interest and is the respective cosmological parameter. We varied each parameter by 5%, as if Δ is too small, the results for the finite difference will be dominated by numerical noise, since the difference between the observables will in many cases be minimal. For and 0 we instead use a larger variation of 10%. For the simulations where is varied we take a slightly different approach, since the neutrino mass scale is less well constrained and we would therefore like to allow for larger excursions. We considered cases with a total neutrino mass of 0.1 eV and 0.2 eV (while ensuring the same squared-mass differences between the three mass eigenstates). Note that we did not keep Ω fixed in these cases.

Data products
We retain several different types of data product from the simulations, giving us enough flexibility to calculate a wide range of observable quantities, but without needing to retain the full particle catalogue and metric perturbations on a grid for many snapshots, which would result in a very large data volume.
HEALP maps -Any calculations that we wish to do on the light cone are simplified by saving the data in spherical shells about a pre-specified observing location, instead of on a Cartesian grid for many snapshots in time. We use the HEALPix pixelisation to store the data for each field and for each shell. These maps are saved at a spatial resolution in radial slices of comoving width 1.75 Mpc/h and cover both the full-sky and the pencil-beam light cone. The side of the maps varies as a function of the distance from the observer, always maintaining a sampling of the fields close to the resolution of the simulation and reaching a maximum of side = 2048.
The quantities saved on these maps are the scalar potential , the peculiar velocity field (coarse-grained at the scale of the pixel) projected onto the line of sight, the CDM+baryon (cb) density, and the neutrino ( ) density. These quantities are calculated, as usual in gevolution, on a Cartesian grid, but are interpolated onto an appropriate set of HEALP maps around the observer at each time step using trilinear interpolation, and then saved to disk.
In Fig. 2, we show examples of these HEALP maps for both the cb and the fields. On the left panel we have these saved for the full-sky light cone at a redshift of ≈ 0.1, and on the right we show the pencil-beam light cone at a redshift of ≈ 1.
Power spectra -For each simulation, we also store a set of power spectra for 14 different redshift slices: = 50, 30, 10, 4, 3, 2.5, 2, 1.5, 1, 0.75, 0.5, 0.25, 0.1 and 0. Thanks to gevolution's ability to calculate all metric data for the entire simulation, we are able to store power spectra for a multitude of variables. In this case, we extract power spectra for both the scalar and vector gravitational potentials, as well as for the CDM+baryon density and the total matter density (i.e. including the massive neutrinos). This results in 14 × 4 = 56 power spectra per simulation, which can be used directly without needing to reanalyse the simulation data.
In Fig. 3 we plot the mean of all of the CDM+baryon power spectra over the 34 random realisations of the fiducial cosmology at ≃ 0, plus confidence intervals showing the expected sample variance (estimated from the simulations) and the error on the mean. We also plot theoretical predictions for comparison, calculated using a linear power spectrum from CLASS (Blas, Lesgourgues & Tram 2011) and a non-linear power spectrum model from HMcode (Mead et al. 2016). The deviations at large scales are due to the different gauges (Poisson gauge for gevolution and the linear prediction of CLASS, while HMcode uses synchronous gauge), while there is good agreement on intermediate and small scales (up to the expected non-linear corrections). Fig. 4 we present the angular power spectra computed using full-sky light cones from the 34 random realisations of the fiducial cosmology. We employ the simple quadratic estimator

Angular power spectra -For completeness, in
where ℓ are the spherical harmonic expansion coefficients of some scalar field projected onto the unit sphere (Leistedt et al. 2013). For this example, we took = cb and performed the analysis using the HEALP anafast routine, which is suitable for the full-sky case. The shaded region in Fig. 4 corresponds to the standard deviation ofˆℓ over the 34 samples, with the theoretical angular power spectrum ℓ computed as where cb ( , ) is the power spectrum accounting for baryons and cold dark matter only, smoothed by the cloud-in-cell (CIC) kernel, as discussed in Sec. 3.2. We projected each map at ≈ 0.58, which is equivalent to a comoving distance of ( ) = 1500 Mpc/ℎ. In this case, the theoretical calculations are broadly consistent with the measured angular power spectrum, but a slight amplitude offset can be seen, particularly for the non-linear theoretical curve at ℓ 100. From several consistency checks, we have found a contribution coming from the mildly-and fully-nonlinear scales that leaks to large angular scales from the projection integral, Eq. (3). By comparing the output angular power spectrum from CLASS, for CDM+baryons, the same offset was observed. However, as we can see, both linear and non-linear predictions match the mean value of the simulations at a level of 10% for ℓ 100. On the other hand, the discrepancies at higher ℓ are primarily caused by resolution effects.
Ray-traced quantities -Given a set of HEALP maps of the gravitational potential on the light cone, it is straightforward to construct other interesting quantities that are linearly related to    ≈ 0.58. The mean and standard deviation over the 34 random realisations of the fiducial cosmology are shown as the solid red line and blue shaded region respectively, while the orange region shows the corresponding standard error of the mean (as in Fig. 3). We also show theoretical predictions based on the linear power spectrum for CDM and baryons from CLASS, and the non-linear power spectrum from HMcode, convolved with the CIC kernel (see Sec. 3.2). Bottom panel: Fractional difference between the theoretical predictions and the mean. The light and dark grey regions show 10% and 5% differences between the theory and the mean respectively. or its time derivative. Examples include the integrated Sachs-Wolfe effect and the weak-lensing potential. Using the Born approximation, HEALP maps of these quantities can be computed directly in pixel space by adding together maps of with appropriate weights. For instance, the weak-lensing potential, which is defined in Lewis & Challinor (2006) as can be constructed by such a procedure, explained in more detail in Lepori et al. (2020). Here we make the assumption that ≈ , which is an excellent approximation here, and neglect the effect of frame dragging. The HEALPix maps can then be easily converted into linear weak-lensing convergence and shear maps by using where Δ is the Laplacian, and the derivatives are taken on the map. The convergence (we use the subscript " " to distinguish this term from Doppler magnification, as discussed below) and the shear parametrise the amplification matrix, that relates lensed images to unlensed ones if lensing is treated linearly (Bartelmann & Schneider 2001). It is worth pointing out that non-linear ray tracing is also possible with the data, e.g. using the full machinery developed in Lepori et al. (2020), although we do not pursue this here.

RESULTS
In this section, we show one of the many ways these simulations can be used to construct relativistic observables by using the Doppler magnification effect as an example. This was first highlighted and investigated as an observable in its own right by Bonvin et al. (2006); Bonvin (2008). Later works have shown that this signal should be detectable with modern day optical and radio surveys, (e.g. Bonvin et al. 2017;Andrianomena et al. 2019). A perturbative expression was derived for the cross-correlation between the Doppler magnification signal and the matter density in Bacon et al. (2014). Doppler magnification is a relativistic effect that is caused by the relative motion of the source and the observer, which correlates with matter density as sources will tend to fall towards areas of high density. It is an inherently (special) relativistic effect, and while it is also possible to derive it from Newtonian simulations, this requires a more careful handling of gauge issues etc. than is needed here. To illustrate their utility, we will go through the steps in calculating this signal within our suite of simulations and then compare with the perturbative results from Bacon et al. (2014).

Doppler magnification
In the linear weak-lensing regime, the true shape of the source is related to the observed image through the Jacobi map, where¯is the angular diameter distance to the source in the background metric, and is the amplification matrix given in Eq. (7).
Doppler magnification appears if one uses the observed redshift as a distance indicator, i.e. the Jacobi map is written at fixed observed redshift as where one defines Here is the peculiar velocity of the source, and and are the redshift and the comoving radial distance to the source, respectively. We define the direction of to be pointed from the observer to the source. It can be seen that galaxies that have a velocity vector directed towards the observer will create a negative at low redshift, which means that they will appear smaller in angular size and dimmer when compared to a typical source at the same observed redshift. In contrast, if the source is moving away from the observer, the will be positive, which means it will appear brighter with a larger angular size. This effect comes about because a surface of fixed observed redshift does not coincide with a surface of fixed comoving distance. Here and in the following, we assume that the Doppler contribution from the source peculiar motion is the only relevant redshift perturbation, i.e. we neglect gravitational redshift and other subdominant corrections. Note that all such corrections can be included by combining the appropriate fields when raytracing, however.
Another important point to note from Eq. (10) is that at high redshifts the term in brackets will decrease, causing a lower amplitude of Doppler magnification. Using this fact and the estimate | | ∼ 0 / , one expects that the Doppler magnification is only important on large scales and at low redshift (Bolejko et al. 2013;Bacon et al. 2014). Fig. 2 of Bacon et al. (2014) does indeed show that the Doppler magnification term is dominant at medium-to-low redshifts and wavenumbers (ℓ 1000 at = 0.2, and 100 at = 0.4). The relative importance of the different contributions to the observed convergence can also be judged from Fig. 5, where we show the Doppler magnification signal and the weak lensing signal within our simulation. To calculate the Doppler magnification signal we use Eq. (10) together with maps of the redshift space distortion field that are included in our data products. For the full-sky maps, we show the signal at ≈ 0.1, where it can be seen that the Doppler magnification is much stronger than the weak lensing signal. When we look at the pencil beam at ≈ 2, we see that the weak-lensing effect dominates instead, as this integrated effect becomes stronger with increasing distance.

Density-convergence cross-correlation
Since matter tends to collapse onto massive structures, there is an obvious correlation between the Doppler magnification and the density field. At low redshift, the Doppler magnification of sources on the far side of a large concentration of matter tends to be negative, while the opposite is true for sources on the near side. In order to measure this correlation, we take the average density over a small interval in distance from the observer, [ , + Δ ], and compute the angular cross-power with the Doppler convergence, , evaluated at the far end of the interval. We then compute the average of the resulting cross-power over a larger distance range in order to accumulate a larger total signal. A full derivation of the angular cross-correlation from perturbation theory is presented in Bacon et al. (2014). We present only the final result in Poisson gauge here. We obtain where cb ( ) is the CDM + baryon power spectrum at redshift zero and ( ) is the linear growth factor, 5 defined as ( , ) = ( ) ( , 0). Note that, in a slight variation to Bacon et al. (2014), we employ weights that are uniform in comoving distance instead of weights that are uniform in redshift.
Matter overdensities only have an appreciable gravitational influence over short distances, expected to be somewhere in the region of tens of megaparsecs. Without trying to make an optimal choice, we set Δ = 52.5 Mpc/ℎ for the distance window in which the density is computed, ignoring longer-range correlations. The crosscorrelation signal ℓ ( ′ ) can then be averaged over a broader bin to get the average angular cross-correlation within the bin, which we will denote as ℓ .
To compute this value in our simulations we use the HEALP maps of the line-of-sight peculiar velocity field · and the CDM+baryon density as described in Sec. 2.2. Specifically, for each radial shell we use maps from the two consecutive simulation time steps that together enclose the light cone at the given distance and %For each of these radial shells, we save two additional shells at time steps on either side. This allows us to calculate the data on the null hypersurface by linear interpolation in conformal time. From these data, we create a thin bin of Δ = 52.5 Mpc/ℎ where we sum up the density maps, and then cross-correlate with the Doppler magnification map directly on the far side of the bin, which we compute from Eq. (10). The cross-correlation is calculated using the HEALP anafast function. We then repeat for all thin bins within the thick bin and take the average to obtain ℓ .
To accurately compare the simulations to the perturbative prediction, it is necessary to convolve the power spectrum in the perturbative calculation with the cloud-in-cell (CIC) kernel (Hockney & Eastwood 1981). This effectively gives a smoothed power spectrum which accounts for the fact that the density field in the simulations is coarse-grained at a finite resolution. The expression for the monopole of this smoothed power spectrum is where the CIC kernel CIC ( ) is defined as where is the -th component of and is the Nyquist wavenumber. The smoothed monopole power spectrum of Eq. (12) is then substituted for cb in Eq. (11).
In Fig. 6 we show the distribution of ℓ measured from all realisations of the fiducial cosmology, with the mean shown as a  solid red line and the standard deviation shown as a blue shaded region. We also plot the perturbative prediction calculated from Eq. (11) as a black dashed line. The results are shown for two redshift bins, the first at ≈ 0.1−0.3, and the second at ≈ 0.6−0.8.
In the lower redshift bin, perturbation theory overestimates the signal beyond ℓ 30 or so. This difference is mostly due to nonlinear effects on small and intermediate scales that are not included in the linear power spectrum model, caused by orbit crossings that generate both vorticity and velocity dispersion and at the same time reduce the power in the velocity divergence (Pueblas & Scoccimarro 2009;Hahn et al. 2015;Jelic-Cizmek et al. 2018). As expected, this effect is also present in the density-velocity cross-correlation. Since this bin is at relatively low redshift, non-linear effects are important even at quite low values of ℓ. In the higher redshift bin, on the other hand, the value measured from our simulations fits more closely to the perturbative prediction, although the strength of the signal has decreased by around an order of magnitude by this point. Also note the shift of the peak in the cross-correlation to correspondingly smaller angular scales.
In Fig. 7 we show the numerical derivatives of ℓ by using finite differences (Eq. (1)) of this quantity in the two redshift bins, ≈ 0.1 − 0.3 (upper panels) and ≈ 0.6 − 0.8 (lower panels), for seven pairs of simulations with the same initial conditions but different values of the cosmological parameters , , ℎ, , , and 0 as described in Table 1. In the case of varying 0 , the pair of simulations used are the baseline cosmology and the one where we include perturbations in the dark energy field (see Sec. 2.1). Derivatives of this kind are useful for Fisher forecasting studies, and give a direct measure of the sensitivity of an observable to a given parameter.
Theoretical predictions of the derivatives from perturbation theory are also shown as dashed lines in Fig. 7. As in previous figures, we see differences of a few percent between the perturbation theory and simulated quantities in the lowest redshift bin (upper panels), which is consistent with the growing importance of non-linear effects at these redshifts, as discussed above. At higher redshift (lower panels), these effects are subdominant however, and the agreement between perturbation theory and the simulated quantities is good to within a couple of percent even on smaller angular scales.
A notable feature of Fig. 7 is the wiggle-like features in some of the curves for the higher redshift bin. These are due to shifts in the location of the baryon acoustic oscillations, which could also be seen (albeit at quite a low level) in Fig. 6. While the signal-to-noise ratio of any practical measurement of the Doppler magnification signal is unlikely to be sufficient to detect the BAO feature for the foreseeable future, it is interesting that they can in principle be picked up by this observable.

CONCLUSIONS
In this paper, we have presented a suite of novel general relativistic cosmological simulations. These were run using the gevolutionbody code (Adamek et al. 2016a), which takes into account all relevant relativistic effects, including frame dragging and relativistic neutrino effects. While most of those effects are small in comparison with the conventional 'Newtonian' terms in most large-scale observables, the rapidly increasing precision of upcoming surveys will soon make it impractical to ignore them without risking biases in cosmological parameter estimates.
The suite of simulations contains a total of 53 runs, divided into subsets that are envisioned to have two main uses. The first subset contains 39 simulations that all use the same fiducial cosmology, but vary the random seed used to generate the initial conditions, which essentially gives us different realisations of the same underlying cosmology. Possible applications of this subset include statistical studies of observables, estimators, and data extraction methods, plus some rudimentary kinds of simulation-based covariance estimation.
The second part of this suite consists of 7 pairs of simulations with cosmological parameters that are systematically varied around the fiducial cosmology (which matches the Euclid Flagship 2 cosmology), while maintaining the same (random) initial conditions. These allow us to study the derivatives of any observable that we can calculate with respect to a set of cosmological parameters. Possible applications of this suite include Fisher forecasting, where derivatives of observables are used to estimate the uncertainties on measurements that can be achieved by future experiments.
We have stored a range of data products for each simulation, including all of the metric degrees of freedom and other fields needed to reconstruct any cosmological observable on large scales. These fields have been determined in a spherical coordinate system about a fiducial observer, and can be fed into a ray-tracing algorithm to produce precise predictions of observables on the past light cone. The geometry of the simulations has been chosen to maximise the sky area and depth of the light cones that can be simulated, with the full sky accessible out to = 0.85 and a large area (1930 sq. deg.) available out to = 3.55. These specifications are well-matched to a variety of current and near-future large-scale structure surveys, including the ESA Euclid mission, the Roman Space Telescope, the VRO Legacy Survey of Space and Time, and the Square Kilometre Array. While the simulations do not have sufficient resolution to produce suitable dark matter halo catalogues for these surveys, biased tracers can be painted onto the simulations using other means (e.g. Borzyszkowski et al. 2017;Bull 2017;Witzemann et al. 2019;Yip et al. 2019;Ramanah et al. 2020;Farr et al. 2020) To showcase the potential use-cases of our suite of simulations, we calculate the cross-correlation of the inherently special relativistic Doppler magnification signal, , and the matter density contrast, cb , from our 34 random realisation simulations and compare it to the perturbation theory result from Bacon et al. (2014). We find good agreement with the perturbation theory calculation in a relatively high redshift bin of ≈ 0.6 − 0.8, but find non-negligible corrections in a lower redshift bin of ≈ 0.1 − 0.3, where linear theory overestimates the signal. This is due to non-linear effects that appear on small scales which are not included in the linear prediction, and which (due to projection effects) affect most of the relevant angular scales at the low redshifts where the Doppler magnification signal is largest. While an improvement, replacing the linear matter power spectrum in the perturbation theory calculation with a non-linear power spectrum model does not fully capture these effects (c.f. Figs. 3 and 4). This highlights the value of having fully relativistic cosmological simulations on hand to make predictions of such observables.
We also calculated the Doppler magnification cross-correlation signal for matched pairs of simulations with the same initial conditions but different values for the cosmological parameters. This allowed us to approximate the derivatives of the observable ℓ and therefore see which parameters it is most sensitive to. While we again saw generally good agreement with predictions from perturbation theory, especially in the higher redshift bin, non-negligible corrections remained. Since reasonably any large-scale structure observable can be constructed from our suite of simulations, including many different combinations of cross-correlations and even high-order statistics, it should be possible to make Fisher matrixtype forecasts for a very wide range of surveys and observables using these data.
In conclusion, in this paper we have described a suite of fully relativistic -body simulations, and shown a particular example (the Doppler magnification term in the density-convergence cross-correlation) in which such simulations are needed in order to make accurate predictions for next-generation surveys. While perturbation theory calculations were able to capture most features of the target signal, non-linear effects made few-percent differences at low redshifts. To accurately model these in perturbation theory, one would likely have to include higher-order corrections in redshiftspace which are difficult to compute.