Observing relativistic features in large-scale structure surveys -- I: Multipoles of the power spectrum

Planned efforts to probe the largest observable distance scales in future cosmological surveys are motivated by a desire to detect relic correlations left over from inflation, and the possibility of constraining novel gravitational phenomena beyond General Relativity (GR). On such large scales, the usual Newtonian approaches to modelling summary statistics like the power spectrum and bispectrum are insufficient, and we must consider a fully relativistic and gauge-independent treatment of observables such as galaxy number counts in order to avoid subtle biases, e.g. in the determination of the $f_{\rm NL}$ parameter. In this work, we present an initial application of an analysis pipeline capable of accurately modelling and recovering relativistic spectra and correlation functions. As a proof of concept, we focus on the non-zero dipole of the redshift-space power spectrum that arises in the cross-correlation of different mass bins of dark matter halos, using strictly gauge-independent observable quantities evaluated on the past light cone of a fully relativistic N-body simulation in a redshift bin $1.7 \le z \le 2.9$. We pay particular attention to the correct estimation of power spectrum multipoles, comparing different methods of accounting for complications such as the survey geometry (window function) and evolution/bias effects on the past light cone, and discuss how our results compare with previous attempts at extracting novel GR signatures from relativistic simulations.


INTRODUCTION
The next generation of galaxy surveys -such Euclid, VRO/LSST, and SKA -will be both wide and deep, covering a broad range of redshifts as well as large areas of the sky, therefore mapping out an unprecedentedly large volume of space and time. On the one hand, this will significantly increase the amount of information available for existing types of cosmological analyses, reducing the same variance uncertainties on observables such as the BAO scale, redshift-space distortions, and the lensing shear power spectrum. On the other hand, the sheer size of these surveys will also allow qualitatively different cosmological observations to be made. In particular, they will be large enough to access modes on the order of the matter-radiation equality scale k eq (e.g. Philcox et al. 2020), and possibly even up to the comoving horizon scale k H ∼ (aH). These represent the very largest observable scales in the Universe, where novel observational features of inflationary and gravitational physics arise that cannot be constrained on the smaller scales probed by existing surveys (e.g. Liguori et al. 2010;Camera et al. 2015;Baker & Bull 2015;Alonso et al. 2015; caroline.guandalin@usp.br (2009) has shown that relativistic effects induce odd multipoles to appear in the cross power spectrum of dark matter tracers, a characteristic with no Newtonian counterpart (Bonvin, Hui & Gaztañaga 2016;Gaztañaga, Bonvin & Hui 2017;De Weerd et al. 2020). This is, by itself, a new cosmological observable allowing us to probe the equivalence principle at cosmological scales via the Euler equation (Bonvin & Fleury 2018), as well as the gravitational redshift effect (McDonald 2009;Bonvin, Hui & Gaztañaga 2014). Moreover, by bearing a strong dependence on the Weyl potential, this provides an alternative test for theories of gravity, while the dependence on astrophysical parameters like the magnification and evolution biases opens a new window to a better understanding of the LSS. Other approaches to constraining deviations from GR via the behaviour of the relativistic effects have also been considered, e.g. Lombriser, Yoo & Koyama (2013); Baker & Bull (2015).
In this paper, we develop the basic building blocks of an analysis pipeline that is capable of extracting the relativistic effect signatures from large-scale structure data. As previously mentioned, standard LSS analysis techniques often rely on Newtonian assumptions or the distant-observer approximation, and so it is necessary to adapt them in order to account for the relativistic effects. Relativistic effects also introduce additional dependencies on the astrophysical properties of the source galaxy population(s) that must be accounted for, such as the magnification bias and evolution bias. Using a mock dark matter halo catalogue extracted from the past light cone of a fully-relativistic N-body simulation generated by the gevolution 1 N-body code, we show how these complications can be overcome in the case of relatively idealised catalogue data, with a view to later extending our pipeline to more realistic scenarios.
For the sake of simplicity, we focus only on the detection of odd multipoles caused by relativistic corrections to the redshiftspace power spectrum. The relativistic effects that arise in the odd multipoles have the advantage of having a leading-order scaling that goes like H /k, making them easier to detect on scales k H as compared with the O(H 2 /k 2 ) corrections that affect even multipoles. The dipole is the most straightforward to model and detect, and has the advantage of having previously been detected in the two-point correlation function and power spectrum of halos in the RayGal simulation 2 by Breton et al. (2019) and Beutler & Di Dio (2020), respectively, at low redshift. This makes it a suitable target for comparison, although we choose to study higher redshifts of around z ∼ 2 − 3 in order to differentiate our paper from these previous works.
This paper is organized as follows. In Section 2 we review the theory of relativistic effects in the two-point statistics of biased tracers. In Section 3 we describe the gevolution light-cone simulation used in this analysis. In Section 4 we review the Fast Fourier Transform (FFT) estimator for the power spectrum multipoles, and present our results in Section 5. Finally, we conclude in Section 6. For the sake of completeness, we also include Appendix A, which explains the details of the halo catalogues derived from the simulated light cone, and Appendix B, where we review the standard method to account for the window function and present some additional results from our measurements.

RELATIVISTIC EFFECTS IN THE POWER SPECTRUM
Contrary to the simplistic view of N-body simulations, which give us the three-dimensional positions of objects at a fixed time slice, the true observed quantity in a galaxy survey is the number of dark matter tracers (e.g. galaxies or halos) N(z,n) in a pixel given by a solid angle dΩ around a directionn = (θ, ϕ), defined with respect to the observer's line of sight (LOS), and at a redshift bin [z, z+dz] (Bonvin & Durrer 2011;Bonvin 2014). The number overdensity of some tracer α can thus be defined as where the equality is obtained by relating the number counts with the number density as n(z,n) ≡ N(z,n)/V(z,n). In the above equation,N α (z) is the selection function of the tracer α, obtained by angular averaging over the tracer number count. The quantities defined in equation (1) are in redshift space, meaning that they are characterized by the observed (comoving) coordinates s = (s, θ, ϕ), with the radial comoving coordinate s being connected to the observed redshift by some cosmological model 3 . The standard treatment (Kaiser 1984), relating the number of sources in a perfect Friedmann-Lemaître-Robertson-Walker (FLRW) universe with the truly observed density field via the conservation of number counts, gives rise to the so-called redshiftspace distortions. This allows us to relate the theoretical predictions in a homogeneous universe with the observed quantities with the addition of departures from the perfect FLRW metric.
In Kaiser (1984), corrections to the angular pair of coordinates (θ, ϕ) are not considered, and perturbations to the radial coordinate s come solely from the peculiar velocities of the sources. Even though it describes satisfactorily observations limited to subhorizon scales, where the Newtonian treatment is well suited, this is not a truly observed quantity, as it is gauge-dependent. Furthermore, future galaxy surveys and cosmological observations that rely on the largest (near-horizon) scales demand a proper treatment of the LSS clustering. At smaller scales, the improved sensitivity will also hold the potential for a detection of subleading corrections (for an example see Saga et al. 2020).
Relativistic corrections that appear by considering the covariant definition of redshift have been widely developed in the past decade, and became a paradigm to study large cosmological scales. In addition to solving well-known gauge issues manifested at these scales, it accounts for a number of effects with no Newtonian counterpart. For instance, gravitational redshift and lensing effects are concisely included in equation (1), and we refer the reader to equation (3.23) of Yoo (2014) and equation (16) of Bonvin (2014) for its full expression.
By collecting the terms proportional to v · n we end up with (Bonvin 2014;Clarkson et al. 2019): where is called Doppler term, H −1 ∂ r (v · n) is the standard Kaiser term, H = aH is the comoving Hubble factor, s α ∝ ∂ ln r ln(r 2 φ α ) is called magnification bias, is the evolution bias and b α is the linear bias. With the exception of the true density perturbation δ α , all other terms appear due to departures from a perfect FLRW universe. Within the linear theory, we can relate quantities in configuration space with their Fourier counterpart to arrive at the main equation with µ k ≡ (k ·r) to keep the explicit dependence with the LOS.
Assuming that all objects in the survey posses the same LOS, i.e.k ·r = µ is a constant (flat-sky approximation), the crossspectrum P (s) αβ (k) = δ α (k)δ * β (k) of two tracers α and β is given by: In this equation, α and β refers to distinct tracers, which could be different types of galaxies or dark matter halos of different masses, f is the growth rate, parametrised by f (z) ∼ Ω m (z) γ , with γ being the growth index, and P (r) (k) is the matter power spectrum in real space.
In this case, isotropy is broken by the choice of LOS and we can expand P where Neglecting the quadratic terms O(H /k) 2 , the coefficients of the expansion are given by 4 In the absence of these quadratic corrections, the monopole, quadrupole and hexadecapole are the same as in the Newtonian case. Still, the imaginary term appearing from the relativistic corrections in equation (2) gives rise to the dipole term manifested in the cross-spectrum of LSS tracers: While it scales as H /k for the cross-correlation of LSS tracers, a fact that makes this signal a smoking gun for relativistic effects in the galaxy clustering, it is identically zero for the auto-correlation. We also call the reader's attention to the fact that this dipole term is anti-symmetric, meaning that δ α (k)δ * β (k) = − δ β (k)δ * α (k) . In Figures 1 and 2 we illustrate the dipole term in both the Fourier and configuration spaces, respectively, for three linear and evolution bias differences (different colours) at a fixed redshift z = 2.
In what follows we explore the detection of (14) in a relativistic simulation of a light cone, described in Section 3. Since we will be dealing with dark matter halos, the magnification bias s α in the Doppler term vanishes. Therefore, in addition to the linear bias of the halos, the remaining parameter entering the theoretical predictions is the evolution bias (4). The procedure for fitting b e from the halo samples is described in Appendix A3, with the results discussed in Section 4.

SIMULATION
In this work we make use of a large N-body simulation performed with the relativistic code gevolution (Adamek et al. 2016a,b). The simulation has a comoving volume of (2.4 Gpc/h) 3 with dark matter particles of mass 2.64 × 10 9 (M /h), and represents a typical ΛCDM cosmology: h = 0.67556, ω b = 0.022032, ω cdm = 0.12038, T CMB = 2.7255 K, A s = 2.215 × 10 −9 , n s = 0.9619, N ur = 3.046, and N ncdm = 0. In order to avoid replications in the light cone, the pencil beam was carefully oriented in the periodic domain. The initial conditions for the simulation were set at a redshift of z = 127.
Unlike the standard approach to building light cones ( Figure 1, but for the cross-correlation function dipole of different tracers at redshift z = 1.9. Differences in the linear and evolution bias are shown in the legend. Solid lines represent the case where there is no magnification bias s α = 0, whereas shaded regions represent the effect of different magnification biases among the tracers. Dotted lines shows the limiting case where s α is smaller than s β by 40%, whilst the dashed ones shows the opposite case, with s α larger than s β by a factor of 40%. of generating many simulation snapshots with a sufficient small redshift step between them to avoid time discretization effects in the final light cone, the light cone output from gevolution records particle positions and velocities on the fly. During the simulation, particles are identified that are within a proper comoving distance interval from a predefined observer that would cause them to be placed in the final catalogue. These particles are then shifted by a fractional time step and recorded on the null FLRW hypersurface given by the past light cone of the observer. Hence, there are no time discretization artefacts and no need to generate an enormous amount of snapshots to build the light cone. In our case, no replications whatsoever were performed in order to cover the whole light cone volume, which has the advantage of removing any concerns about spurious correlations on large scales due to periodicity for example. The gevolution code does not employ the adaptive mesh refinement (AMR) method and thus has a low accuracy at small scales. However, while AMR can improve the 1-halo term by better resolving halo substructures, it does not significantly impact the large scales dominated by the 2-halo term, which is the focus of this work. As will be pointed out in Section 3.2, all sub-halos are discarded in our analysis in any case.

Ray tracing
We apply a ray tracing algorithm to our simulation as a postprocessing tool. The algorithm was previously described in Lepori et al. (2020), but we give a brief review of it here.
The purpose of the ray tracer is to add extra information on source objects within the simulation to the catalogue, such as their angular diameter distance (D A ) relative to a specific observer, the respective observed redshift (z), or the ellipticity ( ) which is closely related to the weak-lensing shear (γ). In contrast to the more common case where ray tracing is applied to Newtonian Nbody simulations, in gevolution the metric perturbations and the source positions are both provided in Poisson gauge, which makes the treatment of gauge issues transparent. Our algorithm also does not rely on the Born approximation to model the light path. Importantly, incorrectly modelling the lensing probability distribution function can lead to errors in estimating cosmological parameters, as shown in (e.g.) Adamek et al. (2019).
The algorithm is similar to the one presented in Breton et al. (2019), and works by integrating the geodesic equations backwards in time from the observer to the source of interest on the observer's past light cone. A physical definition of source, such as a halo or a dark matter particle, is required, as a 4-velocity vector is needed to define the source's rest frame. This allows us to get the observed redshift of the source in a gauge-independent way. For each of these sources, we use the background FLRW model to give us the initial direction vector (n) for each light ray towards a source. We then integrate backwards in time with the fully perturbed metric until the light ray reaches its closest approach to the event on the light cone. At this point, we can now calculate a "deflection angle" by which the initial n must be corrected to achieve a closer approach to the source. We repeat this process several times until suitable convergence is achieved.
This process works well in the weak-lensing regime, as only a single null ray exists between the observer and each source. In the strong lensing regime, multiple images can be formed, which complicates matters. The number of sources where this phenomenon is observed is negligible however, and so we concentrate only on weak lensing. Since strong lensing will only affect our results on very small scales where an image could be duplicated, this choice has a negligible impact on our analysis.
Ray tracing is the key step in properly incorporating relativistic corrections in our analysis. For example, instead of using the redshift output directly from the halo finder which would only include the background expansion and the Doppler correction, we are able to use the 'observed' redshift, which includes all relativistic effects. We can also calculate the perturbed position of sources on the sky, which is important for any n-point correlation calculations done using the catalogue. The algorithm also output D A and both the real and the imaginary parts of the shear component separately (γ 1 + iγ 2 − 4 ), although these are not needed in the current analysis.

Halo catalogue
From the real space particles, the halo catalogue was created with the Rockstar halo finder (Behroozi, Wechsler & Wu 2012), using a friends-of-friends (FOF) algorithm with linking length b = 0.28 in order to detect over 10 7 halos in the light cone.
After going through the ray-tracer algorithm, which is crucial to connect the halos and the observer, the perturbed threedimensional positions of halos were obtained and the resulting file consists of three mock surveys contained within the range 0.0 z 7.1, with different survey areas. The survey that will be used in this work spans the range of comoving look-back distance from 275 up to 4560 Mpc/h.
We limit ourselves to the high-redshift region between z min = 1.7 and z max = 2.9, with redshift bins of size ∆z = 0.4 which kept the variation of the growth function within the 5% limit 5 . Each redshift bin has an effective volume of ∼ 0.7 (Gpc/h) 3 given the chosen cosmology and the sky fraction f sky ∼ 0.01. After this redshift selection, we were left with 8.5 × 10 6 dark matter halos. The high-redshift binning was chosen to deliver a reasonable volume necessary for the observation of the relativistic features at large scales, giving an effective fundamental mode of k F = 2π/V 1/3 ∼ 7 × 10 −3 (h/Mpc). In a future work we will present the results of the same analysis, but in the full-sky case. The current survey area of ∼ 400 deg 2 is compatible with the current survey areas available for a cross-correlation analysis (Zhao et al. 2020).
The final halo catalogue was then separated into three halo samples per redshift bin, of different masses such that, at each redshift bin, the number of halos was the same in each sample. The main properties of these samples are detailed in Table 1. Because more massive halos are expected at lower redshifts, the effective redshiftz of each halo sample varies slightly, but only by less than 0.5%; therefore, we considered the values shown in the Table as the respective central redshift. The biases were computed by fitting the ratio between the real space power spectrum of the halos and dark matter (see Appendix A for a throughout discussion and comparison with the Tinker et al. (2010) fitting function). The halo population incorporating all halos is referred to as H all in what follows.

POWER SPECTRUM MULTIPOLE ESTIMATOR
To compute the power spectrum multipoles we make use of the standard approach proposed by Yamamoto et al. (2006);Bianchi et al. (2015) and Scoccimarro (2015) (for pioneering work see also Yamamoto, Nishioka & Taruya 2000), which we dub YBS estimator. It is built upon the practical algorithm developed by Feldman, Kaiser & Peacock (1994) to optimally estimate the power spectrum of galaxy surveys with a varying selection function. As mentioned in section 2, the selection functionN encodes the spatial modulations of the mean number density of objects. For both spectroscopic and photometric surveys, the selection function accounts for all non-cosmological effects, being sensitive, for example, to the different intrinsic brightness of galaxies.
The selection function gives an estimate of the probability that a galaxy brighter than a certain threshold, at a distance s, is included in the sample. Hence, it is intrinsically related to the notion of luminosity function Φ(L) (Martínez & Saar 2002). In Wang et al. (2020) a clear example of such fact is given, with the luminosity function of eBOSS quasars (QSO) being used to fit the QSO number density and derive the evolution and magnification biases.
To resume the construction of the estimator, N X (x i jk ) denotes either the count-in-cells of the the random catalogue, X = r, or the data (halo) catalogue, X = h, where x i jk is the position of each cell in a three-dimensional grid obtained by a mass assignment scheme, e.g. Nearest Grid Point (NGP), Cloud In Cell (CIC), or Triangular Shaped Cloud (TSC). In this analysis we consider the simplest NGP assignment.
We begin by defining the weighted galaxy fluctuation, or the overdensity field 6 , as where is the number density which will be written as a grid, after mass assignment scheme is chosen. Therefore, in practice n h (x) = N h (x i jk ) is the count-in-cells grid and n r (x) is the corresponding quantity for the random catalogue, which is obtained by randomly sampling α −1 times more objects within the survey volume, with the same selection function as the real data.
The results presented here do not employ a weighting scheme, i.e. w(x) = 1, and we follow Jeong (2010) for the implementation of the quadratic estimator. The normalization factor will be given by and the shot noise, only relevant for the monopole term, will be P shot ≈ x y z 1 + α α where i ≡ L i /n i is the size of each cell dimension in units of Mpc/h. In this work we choose i = 10 Mpc/h. After the construction of these quantities, the power spectrum multipoles can be obtained by the YBS estimator, which we now briefly discuss. The whole idea of this method relies in generalizing the power spectrum to local regions in space, where statistical homogeneity may be assumed. These regions are defined by a single middle line-of-sight d = (s 1 + s 2 )/2, as shown on the left of 6 It is more instructive to write F(x) = w(x)[n(x) −n(x)] = w(x)n(x)δ(x) = W (x)δ(x), where we call W (x) the window function. Then, in Fourier space F(k) is the convolution of the window with the density contrast: F(k) = (2π) −3 ∫ d 3 qW (k − q)δ(q), and one can show that F(k)F(−k) = (2π) −3 ∫ d 3 q |W (k − q) | 2 P(q) + ∫ d 3 xw 2 (x)n(x). Therefore, in Feldman, Kaiser & Peacock (1994) it is considered the overdensity field divided by the magnitude of the window function, W = [ ∫ d 3 xW 2 (x)] 1/2 , which we called N (Jeong 2010). Figure 3. Then, the corresponding power spectrum at this region is Notice that e ik 1 ·s 1 e −ik 2 ·s 2 = e id ·q e −ik ·s , where the Fourier transform of the local configuration is shown on the right of Figure 3. With this change of coordinates, one can see that the local power spectrum can be obtained by taking the Fourier transform of the first component of the local correlation function ξ(s 1 , s 2 ) = ξ(s, d). Finally, the multipoles of the local power spectrum can be obtained from the Legendre expansion Lastly, the inversion of this relation yields the power spectrum multipoles: where the brackets correspond to an average over k-shells and S is the shot-noise term, only relevant for the monopole.
In order to speed up the computation of the multipoles by means of FFTs, the YBS estimator takes the end-point LOS s 1 . It is worth mentioning that this LOS intrinsically generates odd multipoles which may impact the signal we are trying to measure. Hence, this must be accounted for in the window function, as discussed in Appendix B.
With the adoption of this LOS, the monopole can estimated aŝ is the Fourier transform of the overdensity field with no weight, and the dipole is obtained bŷ where and For higher order multipoles, we refer the interested reader to Bianchi et al. (2015) and Beutler, Castorina & Zhang (2019) for the even and odd ones, respectively. Finally, to compute the cross-dipole, F 0 (k) and F 1 (k) are built from the first and second tracers, respectively, and the normalization becomes (Beutler & Di Dio 2020): Still for the specific case of dark matter halos, the magnification bias is identically zero, s = 0, and the only term accounting for the mass function variations is the evolution bias, which explores its dependency with time: this is related to the fact that halos can merge to form more massive structures and, thereby, their number counts are not conserved. The evolution biases b i e (z) of each halo sample i = {H 0 , H 1 , H 2 , H all } considered in this work, within each redshift bin, is shown in Figure 4. The procedure to compute b i e (z) is described in Appendix A and is based on the work of Beutler & Di Dio (2020).
We work with 29 bandpowers (Fourier bins) linearly spaced between k min ≈ 0.006 and k max ≈ 0.157 with ∆k ≈ 0.005. As already mentioned in Section 4, we work with three-dimensional grids containing n x = n y , n z < n x cells of side i = 10 Mpc/h, with the number of cells n i varying between the different redshift bins. This binning was chosen to deliver a less noisy measurement at large scales. Figure 5 shows the monopole estimated from the four halo samples in the first redshift bin (z ≈ 1.89). Apart from the amplitude of the monopole, not much change occurs between different redshifts, hence we only show the first z-bin here. Our error  Table 1, at each redshift bin considered in the analysis. Notice that, even though the number density is approximately the same for all the halo populations, the intrinsic evolution of each halo population with the redshift gives rise to very different evolution biases. The computation is described in Appendix A3, following Beutler & Di Dio (2020). bars are estimated from the standard deviation of 100 log-normal mocks generated with the same characteristics of the original data: box dimensions, evolution and linear biases, selection function and survey mask. For this last step, we first generated the log-normal mocks for the whole box encompassing the different redshift bins of the light cone, and then applied the proper mask to select the specific angular region. However, because these error bars only quantify the variance of the estimator, the particular survey features do not matter much, and thus it should be possible to compute the variance of mocks inside the box, with the only caveat of following the mean number density of tracers to properly incorporate the shot noise. The impact of the window function in the standard deviation of the log-normal samples generated minor changes at very large scales, and we do not explore this further.
To what concerns the variance of our measurements, the window function had a negligible impact. Nonetheless, to avoid any sort of complications, we explore the asymmetry of the relativistic signal as suggested in Beutler & Di Dio (2020): by computing ∆P 1 = P αβ 1 − P βα 1 it is possible to isolate the relativistic contribution and get rid of the impact of the window function, which is symmetric. As pointed in Section 2, the Doppler term is anti-symmetric, δ α (k)δ * β (k) = − δ β (k)δ * α (k) , and thus ∆P 1 ∼ 2 δ α (k)δ * β (k) . Consistently, the theoretical prediction for the last crosscorrelation H 2 × H all is positive, as ∆b = b 2 − b all > 0. However, as shown in Figure 6, we conclude that no detection can be claimed with this pencil-beam light cone. Very similar results were obtained for the other two redshift samples z 1 ≈ 2.29 and z 2 ≈ 2.69.
Still, we point out the possibility of exploring optimal weighting schemes to enhance the signal in the light of the work carried out by Castorina et al. (2019). Lastly, the advantages of employing different tracers, coupled with the low densities of the catalogues, as well as the need for robust statistics, suggests the use of optimal weights for a more efficient combination of tracers (Abramo, Secco & Loureiro 2016;Montero-Dorta et al. 2020).

CONCLUSIONS
We have performed a power spectrum multipole analysis on data from a light cone generated from a fully-relativistic N-body simulation. We focused on the dipole signal in the cross-correlation between different dark matter halo sub-populations, which is a purely relativistic (non-Newtonian) effect. The simulation was generated by the gevolution code, which employs a novel ray-tracing Figure 5. Monopole estimated from the four halo samples in the first redshift bin:z ≈ 1.89. Error bars are computed from the standard deviation of 100 log-normal mocks generated with the same box dimensions, evolution and linear biases, selection function and survey mask. Solid lines represent the convolution of the theory with the window function (see Appendix B). The bottom panel shows k P 0 (k) for the less and more massive halo samples (lower and upper curve, respectively) for better visibility of the errors and the larger scales; it also shows the theoretical monopoles without considering the window function (dash-dotted lines). method to connect the halos with the observer, and which is capable of incorporating all relevant general relativistic effects on cosmologically-relevant distance scales. We showed in detail how the survey window function and quantities such as the evolution bias can be estimated on the past light cone, allowing a rigorous comparison with gauge-invariant theoretical calculations at linear order.
Similar studies of relativistic observables in simulations have been made in the past. For example, Breton et al. (2019) and Beutler & Di Dio (2020) used the full-sky RayGal simulation, which is limited to the redshift range of 0.05 < z < 0.465, with an effectivez ∼ 0.341. While the simulation that we based our study on in principle covers the redshift range 0 ≤ z ≤ 7.1, our analysis focused on a particular high redshift bin in the range 1.7 ≤ z ≤ 2.9, covering a sky fraction of only f sky = 0.01. This is a similar sky area to the overlap region between different LSS tracers (luminous red galaxies and emission line galaxies) in the multi-tracer analysis of the final eBOSS data (see Table 2 of Zhao et al. 2020), although these data are from lower redshift, z ∼ 1.
While we were able to robustly test our analysis methods using these simulated data, no conclusive detection of the dipole signature was possible due to the limited volume of the redshift bin, a challenge that is of paramount importance for current surveys too. Beutler & Di Dio (2020) studied the possibility of subtracting various contributions to the total signal in order to iso-late the Doppler contribution and remove sample variance. Since the Doppler term is expected to increase in amplitude with redshift, one could also consider developing an optimal weighting scheme to enhance the signal and improve the prospects of detection . We leave this, and other schemes (Abramo, Secco & Loureiro 2016;Abramo & Bertacca 2017;Montero-Dorta et al. 2020) to enhance detectability of the signal, to be explored in future work however.
We did not incorporate wide-angle effects in our modelling, as they are not relevant for the solid angle and redshift range of our analysis. A careful account of these effects should also be explored in the context of wider survey areas however, particularly in the case of future surveys such as Euclid, LSST, and SKA, which are expected to cover an appreciable fraction of the sky.
Similarly, integrated effects (e.g. lensing), while fully included in our mock data, were neglected in our analytical model, but are known to impact large angular scales. Despite the Doppler term being the largest contribution to the relativistic effects for our particular setup, non-local terms should also be modelled and properly included for analyses that go to larger scales.
In this paper, we have limited our analysis to a single highredshift bin with a relatively narrow survey area, and have pursued only a limited set of observables, i.e. the multipoles of the relativistic power spectrum. In future work, we will relax these limitations by moving to larger survey volumes more representative of the next generation of large-scale structure surveys, while also including wide-angle and integrated effects, and extending our analysis to two-point correlation functions and multipoles of the relativistic bispectrum.
withρ m,0 the comoving background matter density today, f (σ) the multiplicity function, and σ the overdensity variance smoothed in a sphere of radius R. The multiplicity function can be computed analytically from the spherical collapse model (Press & Schechter 1974) or the ellipsoidal collapse (Sheth & Tormen 1999), or from numerical fits (Tinker et al. 2008). The RayGal simulation consists of a set of high-resolution Newtonian N-body simulations, whose halos have been ray-traced to the redshift space position, rendering them with almost all properties of our halos. The RayGal light cone was build from 300 snapshots to avoid time discretization effects.
Our analysis was based on halo masses M 200b ≡ M defined within the density thresholds of ∆ = 200, whose correspondence with the parameters fit of the Tinker mass function (Tinker et al. 2008) (solid lines in Figure A1) is straightforward. In Corasaniti et al. (2018), RayGal halo mass functions were computed from the snapshots and were based on the Sheth-Tormen (Sheth & Tormen 1999) fit, with the halo identified with the spherical overdensity (SO) method, and thus the halos are more closely connected to the ellipsoidal collapse employed in the Sheth-Tormen fit (Desjacques, Jeong & Schmidt 2018).
However, for the RayGal light cone, halos were identified via a friends-of-friends (FOF) algorithm, just like in our catalogue. In Smith et al. (2017), differences with numerical fits seem at the lowmass end are also present, and they conclude that such discrepancies are associated with the comparison between different halo finder methods (SO and FOF). They computed the mass function using the SO correspondent, and just as in our case, found the same behaviour at low masses. In the catalogue employed in our analysis, ∼ 1.1 × 10 6 halos with M vir ∈ [0.518, 4.862] × 10 12 (M /h) were discarded for having their respective M 200b null. As pointed out in Smith et al. (2017), small overdensities in large FOF groups might be identified as part of the larger group, leading to a lack of such structures.
Such discrepancies are important if one wishes to paint galaxies to the halos via, e.g., a halo occupation distribution. For our current purposes, the lack of a proper function to describe the light cone halo mass function impacted only our ability to predict the linear halo bias (see next section), and thus did not pose an issue for the analysis.

A2 Halo bias
The halo mass function describes the fraction of matter inside dark matter halos. So in order to obtain the correct halo statistics, we must account for their position in space. The halo bias, which is defined by the ratio of the halo power spectrum, P hh (k), to the linear dark matter power spectrum, P lin (k) (Tinker et al. 2010 is best understood within the context of the Peak-Background Split (PBS), where the long-wavelength modes enhances the probability Figure A1. Comparison between the halo mass function of the catalogue employed in this analysis (stars), with the one from the full-sky RayGal simulation (dots) from Breton et al. (2019), with both catalogues in real space. Solid curves correspond to the Tinker mass function fit (Tinker et al. 2008), while the bottom panel shows the relative difference between the fit prediction and the mass function computed from the simulations. Vertical dashed line corresponds to the limit of gevolution halos with at least 600 particles. The cosmological parameters in the RayGal simulation that differ from ours are h = 0.72, A s = 2.431 × 10 −9 , Ω m = 0.257, T cmb = 2.726. We stress that the RayGal mass function was multiplied by a 0.1 factor for a cleaner visualization, as the values were very similar. Differences between the Tinker and Sheth-Tormen (Sheth & Tormen 1999) mass functions were minor, so we only present the former.
of forming halos by decreasing the threshold δ c (z = 0) = 1.686 for overdensities which are located at the peak of large-scale (background) fluctuations. It can be either derived from analytical mass functions, giving the Press-Schechter and Sheth-Tormen halo biases, or from equation A2 via numerical simulations. From Tinker et al. (2010), the bias is given by where ν = δ c /σ and A, B, C, a, b, and c are parameters fitted from simulations, depending on the matter perturbations at virialization which is chosen to be ∆ = 200. This phenomenological fit proved to be unsatisfactory for our halo samples, for the reasons described in the previous section.
We proceeded then to the definition of equation (A2), with the polynomial fit b 2 (k) = b 2 1 + b 2 2 k, and considered the linear term as the fit for the linear halo biases, neglecting the scale dependence emerging from non-linear effects in the power spectrum. Notice that the estimated spectra P hh employed in this fit are for the real space halos. The results are shown in Table A1 for each halo sample. This approach worked reasonably well for the first redshift slice, but proved to be inaccurate for the higher redshift samples, where the impact of the window function was larger. The most likely reason is the mode coupling with the window function, present in our real space light cone.
Hence, instead of convolving the theory with the window Figure A2. The comoving number density of halos, normalized by the total amount N 0 of objects inside each mass bin. Solid lines represent the linear fit of equation (A7), whereas dashed lines represent the truen(z) of the simulation, computed by dividing each redshift bin into 50 points. The smooth curves are the result of a cubic spline interpolation, for visual reasons. survey window function (random catalogue) contains 10 8 particles to completely fill the survey region, one possibility is to compute the power spectrum multipoles of the random catalogues. We adopt this approach, obtaining Q (k) by means of FFTs, which are then taken to configuration space. From the straightforward product of equation (B2), the Legendre expansion of (B2) results in (Wilson et al. 2017;Beutler et al. 2017;Beutler, Castorina & Zhang 2019;Beutler & Di Dio 2020): ξ 0 (s) = ξ 0 (s)Q 0 (s) + 1 5 ξ 2 (s)Q 2 + 1 9 ξ 4 (s)Q 4 (s) + . . .
From an inverse Hankel (1D Fourier) transform of equations (B9) and (B10) we finally obtain the convolved power spectrum mul-tipoles of equation (B4) which we use for comparison with the estimated quantities. This paper has been typeset from a T E X/L A T E X file prepared by the author.