LIGER: mock relativistic light-cones from Newtonian simulations

TL;DR

LIGER is a fast method that transforms Newtonian simulation data into relativistic mock galaxy catalogues, including lensing and Doppler effects, aiding survey forecasts and covariance estimations.

Contribution

The paper introduces LIGER, a novel technique to incorporate relativistic effects into Newtonian simulation-based galaxy catalogues efficiently.

Findings

Magnification bias can be measured with high confidence in Euclid-like surveys.

Relativistic Doppler effects can be detected at 5.5σ significance with SKA.

LIGER enables rapid generation of realistic mock catalogues for survey analysis.

Abstract

We introduce a method to create mock galaxy catalogues in redshift space including general relativistic effects to linear order in the cosmological perturbations. We dub our method LIGER, short for `light cones with general relativity'. LIGER takes a (N-body or hydrodynamic) Newtonian simulation as an input and outputs the distribution of galaxies in comoving redshift space. This result is achieved making use of a coordinate transformation and simultaneously accounting for lensing magnification. The calculation includes both local corrections and terms that have been integrated along the line of sight. Our fast implementation allows the production of many realizations that can be used to forecast the performance of forthcoming wide-angle surveys and to estimate the covariance matrix of the observables. To facilitate this use, we also present a variant of LIGER designed for large-volume simulations with low mass resolution. In this case, the galaxy distribution on large scales is obtained by biasing the matter-density field. Finally, we present two sample applications of LIGER. First, we discuss the impact of weak gravitational lensing onto the angular clustering of galaxies in a Euclid-like survey. In agreement with previous analytical studies, we find that magnification bias can be measured with high confidence. Second, we focus on two generally neglected Doppler-induced effects: magnification and the change of number counts with redshift. We show that the corresponding redshift-space distortions can be detected at 5.5$\sigma$ significance with the completed Square Kilometre Array.

Figures

Equations

$$x_{\mathrm{s}}^{\mu }=({\eta }_{\mathrm{s}},{\mathbf{x}}_{\mathrm{s}})=({\eta }_0-{\chi }_{\mathrm{s}},{\chi }_{\mathrm{s}}{\mathbf{n}}_{\mathrm{s}}),$$

$$x_{\mathrm{r}}^{\mu }[{\chi }_{\mathrm{r}}({\chi }_{\mathrm{s}})]=x_{\mathrm{s}}^{\mu }({\chi }_{\mathrm{s}})+\Delta x^{\mu }({\chi }_{\mathrm{s}}),$$

$$x_{\mathrm{r}}^{\mu }({\chi }_{\mathrm{r}})$$

$$\Delta x^0({\chi }_{\mathrm{s}})$$

$$\Delta x^i({\chi }_{\mathrm{s}})$$

$$\mathrm{d}s^{2}=a^{2}(\eta)\Big{[}-\big{(}1+2\Psi\big{)}\mathrm{d}\eta^{2}+\big{(}1-2\Phi\big{)}\delta_{ij}\mathrm{d}x_{\rm r}^{i}\mathrm{d}x_{\rm r}^{j}\Big{]}\;,$$

$$1+z=\frac{(u_{\mu }{p}^{\mu }){|}_{\mathrm{e}}}{(u_{\mu }{p}^{\mu }){|}_{\mathrm{o}}},$$

$$\delta \chi$$

$$\delta x^0$$

$$\delta x^i$$

$$\mathcal{M}={\left(\frac{D_{\mathrm{L}}}{{\bar{D}}_{\mathrm{L}}}\right)}^{-2},$$

$$\mathcal{M}$$

$$\kappa =\frac{1}{2}{\int }_0^{{\chi }_s}\left({\chi }_{\mathrm{s}}-\chi \right)\frac{\chi }{{\chi }_{\mathrm{s}}}{\nabla }_{\perp }^2\left(\Phi +\Psi \right)\mathrm{d}\chi$$

$${\nabla }_{\perp }^2={\nabla }^2-{\left(n_{\mathrm{s}}^i{\partial }_i\right)}^2-\frac{2}{\chi }{n}_{\mathrm{s}}^i{\partial }_i.$$

$$\Phi =\Psi =\varphi ,{\nabla }^2\varphi =4\pi G{a}^2{\bar{\rho }}_{\mathrm{m}}{\delta }_{\mathrm{sim}},\mathrm{and}{v}^i=v_{\mathrm{sim}}^i,$$

$${\delta }_{\mathrm{s}}={\delta }_{\mathrm{sim}}+{\delta }_{\mathrm{RSD}},$$

$${\delta }_{\mathrm{RSD}}$$

$$\delta \ln a={\Psi }_o-(n_{\mathrm{s}}^i{v}_i{)}_o-{\Psi }_{\mathrm{e}}+(n_{\mathrm{s}}^i{v}_i{)}_{\mathrm{e}}-{\int }_0^{{\chi }_s}{\partial }_0\left(\Phi +\Psi \right)\mathrm{d}\chi .$$

$${\delta }_{\mathrm{g},\mathrm{s}}=b{\delta }_{\mathrm{sim}}+\mathcal{Q}(\mathcal{M}-1)+\mathcal{E}(\delta \ln a-\mathcal{H}{\varphi }_v)+{\delta }_{\mathrm{RSD}}.$$

$$\displaystyle{\cal Q}=-\frac{\partial\ln\bar{n}_{\rm g}}{\partial\ln L}\Big{|}_{L=L_{\rm lim}},$$

$$\mathcal{E}=-\frac{\partial \ln {\bar{n}}_{\mathrm{g}}}{\partial \ln (1+z)}.$$

$${\delta }_{\mathrm{g},\mathrm{s}}=(b-1){\delta }_{\mathrm{sim}}+{\delta }_{\mathrm{s}}+\mathcal{E}\delta \ln a+\mathcal{Q}(\mathcal{M}-1),$$

$$n_{\mathrm{g},\mathrm{s}}$$

$$n_{\mathrm{g},\mathrm{s}}({\hat{\mathbf{n}}}_{\mathrm{s}},z)$$

$${\delta }_{\mathrm{g},\mathrm{v}}=-\frac{1}{\mathcal{H}}{\left(\frac{\partial n_{\mathrm{s}}^i{v}_i}{\partial {\chi }_{\mathrm{s}}}\right)}_{\mathrm{e}}-\frac{\alpha ({\chi }_{\mathrm{s}})}{\mathcal{H}{\chi }_{\mathrm{s}}}\left[{\left(n_{\mathrm{s}}^i{v}_i\right)}_{\mathrm{e}}-{\left(n_{\mathrm{s}}^i{v}_i\right)}_{\mathrm{o}}\right],$$

$$\alpha ({\chi }_{\mathrm{s}})={\gamma }_0+{\gamma }_1\mathcal{H}{\chi }_{\mathrm{s}}$$

$${\gamma }_0=2(1-\mathcal{Q})\mathrm{and}{\gamma }_1=1-\frac{3}{2}{\Omega }_{\mathrm{m}}(z)-\mathcal{E}+2\mathcal{Q}.$$

$$a_{lm}^{(i)}=\int d^2\text{mn@boldsymbol}\theta {\sigma }_g^{(i)}(\text{mn@boldsymbol}\theta )Y_{lm}^{\ast }(\text{mn@boldsymbol}\theta ),$$

$$C_l^{(ij)}=\frac{1}{2l+1}\sum _{m=-l}^l{a}_{lm}^{(i)}{a}_{lm}^{(j)\ast }.$$

$$\mu =C_l^{(\mathrm{A})}{\Sigma }_{lm}^{-1}{C}_m^{(\mathrm{A})}$$

Full text

††pagerange: LIGER: mock relativistic light-cones from Newtonian simulations–A††pubyear: 2015

1 Introduction

The advent of galaxy redshift surveys has revolutionised our understanding of the large-scale structure of the Universe and provided us with multiple ways to constrain the cosmological model. Mock catalogues of synthetic galaxies play a threefold role in the analysis of these datasets (Cole et al., 1998; Blaizot et al., 2005; Kitzbichler & White, 2007; Sousbie et al., 2008; Carlson & White, 2010; Merson et al., 2013). (i) They shape theoretical predictions into structures that closely match observations. (ii) They form a straightforward tool to derive biases and covariance matrices of estimators for statistical descriptions of the large-scale structure (e.g. correlation functions or their Fourier analogues). (iii) Related to that, as forecasting tools, they provide key information to designing new surveys by minimising the impact of statistical errors and systematic effects on selected observables.

Since the 1970s, the size of galaxy catalogues has constantly increased in terms of solid-angle and redshift coverage as well as in sampling rate. The next generation of surveys will provide us with the possibility to measure galaxy clustering on scales comparable with the Hubble radius (e.g. Laureijs et al., 2011; Levi et al., 2013; Maartens et al., 2015). Theoretical studies suggest that a number of general relativistic effects might be detectable on these scales. In order to fully exploit the potential of the new datasets, it is therefore imperative to develop analysis tools (and thus mock catalogues) that include these effects. In this paper, we present a method to create mock galaxy catalogues that incorporate relativistic corrections and are built upon the output of either common Newtonian simulations of galaxy formation or semi-analytic models based on standard N-body simulations.

Relativistic effects arise from the fact that we observe galaxies on our past lightcone. The presence of perturbations superimposed to a Friedmann-Robertson-Walker (FRW) background alters the null geodesics of the photons emitted by distant galaxies. In consequence, all the direct observables for a galaxy are different than in a smooth universe: its redshift, angular position on the sky and the flux in any given waveband. Galaxy peculiar velocities, for instance, distort the radial pattern of the galaxy distribution (Kaiser, 1987; Hamilton, 1997). Similarly, magnification due to gravitational lensing modifies the observed number counts in flux-limited samples (Turner, 1980; Turner, Ostriker & Gott, 1984; Sasaki, 1987; Matsubara, 2000). Many recent studies have demonstrated the existence of several additional corrections that, although suppressed on smaller scales, might generate observable signals on distances comparable with the Hubble radius (e.g. McDonald, 2009; Yoo, 2009; Yoo, Fitzpatrick & Zaldarriaga, 2009; Bonvin & Durrer, 2011; Challinor & Lewis, 2011; Bertacca et al., 2012; Jeong, Schmidt & Hirata, 2012; Yoo et al., 2012; Di Dio et al., 2013, 2014; Bonvin, Hui & Gaztañaga, 2014; Montanari & Durrer, 2015; Yoo & Seljak, 2015; Bonvin, Hui & Gaztanaga, 2016; Cardona et al., 2016; Di Dio et al., 2016; Raccanelli et al., 2016a, b, c; Gaztanaga, Bonvin & Hui, 2017). At linear order in the perturbations, these additional corrections include Doppler terms plus Sachs-Wolfe (standard and integrated, see Sachs & Wolfe, 1967; Rees & Sciama, 1968) and (Shapiro) time-delay contributions. Robust models of galaxy clustering on large scales should thus include these modifications that, most likely, will be key to extracting unbiased information on the dark sector of the Universe (i.e. on the nature of dark energy and dark matter) and to improve constraints on primordial non-Gaussianity. This can be done following different approaches. One possibility is to study structure formation using N-body simulations that include dynamical space-time variables in the weak-field approximation (Adamek et al., 2013; Adamek, Durrer & Kunz, 2014; Adamek et al., 2016a, b), within a post-Friedmann framework (Milillo et al., 2015; Bruni, Thomas & Wands, 2014; Thomas, Bruni & Wands, 2015), or by employing full numerical relativity (Bentivegna & Bruni, 2016; Giblin, Mertens & Starkman, 2016). Alternatively, one can correct a posteriori the results of Newtonian simulations to account for lightcone effects (Chisari & Zaldarriaga, 2011). This is feasible because, at linear order in the perturbations, the mathematical description of a pressureless fluid can be formulated so that there is agreement between general relativity and its Newtonian approximation (Haugg, Hofmann & Kopp, 2012; Rigopoulos & Valkenburg, 2015; Fidler et al., 2015, 2016)

In this paper, we follow the latter approach to develop the LIGER (light cones using general relativity) method. As illustrated in Fig. 1, schematically, LIGER takes a Newtonian simulation as an input and, after selecting an observer, outputs the distribution of galaxies in ‘comoving redshift space’ (i.e. as it would be inferred by applying the background metric to the observed galaxy properties). The algorithm combines the original snapshots of the simulation at constant background time to produce the galaxy distribution on the perturbed light cone. This is achieved by using a coordinate transformation that includes local terms and contributions that are integrated along the line of sight.

Multiple efforts have been made in the literature to investigate the detectability of subtle relativistic effects from forthcoming survey data. Generally, these studies are based on the Fisher-information matrix, use idealised survey characteristics and neglect systematics. The ultimate test to discern which relativistic effects will be observable is to apply the very same estimators that are used for the data to the LIGER mocks. This exciting perspective provides the main motivation for our work. The numerical implementation of the LIGER method will be made publicly available in due course.

The paper is structured as follows. In Section 2, we introduce the LIGER method and describe its numerical implementation. In Section 3, we present two straightforward applications of our code. As an illustration of LIGER’s functionality, we first re-analyse a result which has already been discussed in the literature, namely, the impact of magnification bias in the observed cross-correlation of galaxy samples at substantially different redshifts. Subsequently, we discuss the more challenging detection of Doppler terms in the galaxy angular power spectrum at low redshift. Finally, in Section 4, we conclude. Throughout, we adopt units in which the speed of light is one and define the space-time metric tensor to have signature $(-,+,+,+)$. Greek indices indicate space-time components (i.e. run from 0 to 3) while Latin indices label spatial components (i.e. run from 1 to 3). The Einstein summation convention is adopted.

2 The LIGER method

2.1 Theory

2.1.1 Redshift-space distortions

We observe galaxies as they are at the time in which their worldline intersects our past light cone. The comoving location of a galaxy can be inferred from two basic observables: its position on the sky, ${\bf n}_{\rm s}$ (a unit vector defined in terms of two angles), and its redshift, $z$. In fact, these data are sufficient to build three-dimensional maps of the galaxy distribution provided that we assume to live in an unperturbed FRW universe with a fixed set of cosmological parameters. In reality, such a ‘redshift-space’ map gives a distorted portrayal of the cosmic web due to the presence of inhomogeneities (Sargent & Turner, 1977). Galaxies are artificially shifted both in the radial and tangential directions due to their peculiar motions and the bending of the light they emit. These effects are collectively known under the name of redshift-space distortions.

In mathematical terms, redshift-space is characterized by a set of coordinates that ‘flatten’ our past light cone (e.g. Bertacca, 2015). For instance, the null geodesic from an observed galaxy to us can be described in terms of the following conformal space-time coordinates:

[TABLE]

where $\eta_{0}$ is the present-day value of conformal time (i.e. at observation), $\chi_{\rm s}$ denotes the comoving distance (from the observer) of events located along the geodesic in the unperturbed model universe and $n^{i}_{\rm s}=x_{\rm s}^{i}/\chi_{\rm s}$. The full distance to the galaxy corresponds to the observed redshift $z$, in compact notation $\chi_{\rm s}(z)$.

For a given photon path (see Fig. 2), we want to define a mapping from real to redshift space,

[TABLE]

where $x_{\rm r}^{i}$ denotes the actual comoving position located at distance $\chi_{\rm r}$ along the direction $n^{i}_{\rm r}=x_{\rm r}^{i}/\chi_{\rm r}$. [In general, we use the subscripts ‘s’ and ‘r’ to distinguish redshift-space quantities from their real-space counterparts.] Perturbing $x^{\mu}_{\rm r}$ around $x^{\mu}_{\rm s}$ and writing $\chi_{\rm r}=\chi_{\rm s}+\delta\chi$, we obtain, at linear order,

[TABLE]

By using $\chi_{\rm s}$ as the affine parameter for the null geodesic, we write the total derivative along the past light cone as $\mathrm{d}/\mathrm{d}\chi_{\rm s}=-\partial/\partial\eta_{\rm s}+n^{i}_{\rm s}\partial/\partial x_{\rm s}^{i}$. Since ${\mathrm{d}x_{\rm s}^{i}}/{\mathrm{d}\chi_{\rm s}}=n^{i}_{\rm s}$ [to zero order], linear redshift-space distortions can be written as

[TABLE]

The first term on the right-hand side of equation (5) corresponds to the change in the affine parameter while the second one derives from the perturbation of the photon path and has both tangential and radial components. Note that the real-space distance to the galaxy does not coincide with $\chi_{\rm s}(z)$.

2.1.2 Perturbations of the past light cone

In order to compute explicitly all the terms in equations (4) and (5), we need to work out how metric perturbations alter null geodesics. In what follows, we model the matter content of the universe as a collisionless fluid in the single-stream regime. This common assumption provides a suitable approximation on the large scales we are interested in. We use the subscripts ‘e’ and ‘o’ to denote the fluid properties evaluated at the position of the light source (the galaxy) when the photons were emitted and at the location of the observer when the photons were received, respectively. For linear scalar perturbations in the Poisson gauge111The restricted Poisson gauge containing only scalar perturbations is also known as the longitudinal or conformal Newtonian gauge., the space-time metric can be expressed as

[TABLE]

where $a$ denotes the scale factor of the expanding universe while $\Psi$ and $\Phi$ indicate the Bardeen potentials of the inhomogeneities. In this framework, the redshift of a galaxy measured by an observer is

[TABLE]

where $u^{\mu}$ denotes the four-velocity of the matter fluid (we assume there is no velocity bias) and $p^{\mu}$ is the photon four-momentum. By perturbing the photon geodesic around the FRW solution, we derive expressions for $\delta x^{\mu}$ and $\delta\chi$. The final result for the galaxy shift is (see also Yoo, 2009; Yoo, Fitzpatrick & Zaldarriaga, 2009; Bonvin & Durrer, 2011; Challinor & Lewis, 2011; Jeong, Schmidt & Hirata, 2012)

[TABLE]

[TABLE]

where $\mathcal{H}=\partial_{0}\ln a$ and $\chi_{\rm s}$ are evaluated at the observed redshift of the galaxy while $v^{i}$ is the peculiar velocity. Here, local corrections express the Sachs-Wolfe and the Doppler effects. Those integrated along the line of sight derive from gravitational lensing, the Shapiro time-delay and the integrated Sachs-Wolfe effect.

2.1.3 Magnification

Metric perturbations also alter the solid angle under which galaxies are seen by distant observers thereby enhancing or decreasing their apparent flux. In terms of the luminosity distance, $D_{\rm L}$, the magnification of a galaxy is defined as

[TABLE]

where $\bar{D}_{\rm L}$ denotes the luminosity distance in the background model universe evaluated at the observed redshift of the galaxy.

At linear order in the perturbations, we can write (e.g. Challinor & Lewis, 2011; Bertacca, 2015)

[TABLE]

where

[TABLE]

corresponds to the classical convergence and the differential operator $\nabla_{\perp}^{2}$ is defined as

[TABLE]

Note that the magnification includes contributions from different physical effects. We refer to the term proportional to $n^{i}_{\rm s}v_{i}$ as ‘Doppler lensing’. In order to model statistical observables (e.g. correlation functions) on small-scales, it is acceptable to replace the two-dimensional Laplacian $\nabla_{\perp}^{2}$ with the three-dimensional one $\nabla^{2}$, so that $\kappa$ can be expressed in terms of the matter overdensity using the Poisson equation. In fact, the effective lensing weight, $(\chi_{\rm s}-\chi)\,\chi/\chi_{\rm s}$, varies on scales comparable to the Hubble radius and the line-of-sight integral that defines $\kappa$ heavily suppresses the contribution of radial Fourier modes with smaller wavelengths (e.g. Kaiser, 1992). Although this approximation has been implemented to produce full-sky mock catalogues (e.g. Fosalba et al., 2008, 2015), we do not use it since we want to study galaxy clustering at wide angular separations.

2.1.4 Link with N-body simulations

To evaluate $\delta x^{\mu}$, $\delta\chi$ and ${\cal M}$, we need to compute the gravitational potentials appearing in equation (6) as a function of space and time. Since we want to apply our results to simulations, we need to derive the potentials starting from the particle distribution in the computer models. This corresponds to using the matter density contrast in the synchronous comoving gauge, i.e. $\delta_{\rm sim}\equiv\delta_{\rm syn}$. Fortunately, to linear order in the perturbations and for a pressureless fluid in a universe with $\Lambda$CDM background, the source equation for $\Psi$ in the Poisson gauge can be re-written in terms of $\delta_{\rm syn}$ as the standard Poisson equation (e.g. Chisari & Zaldarriaga, 2011; Green & Wald, 2012). Therefore, the complete dictionary we use to translate from the simulations to the Poisson gauge is:

[TABLE]

where $G$ denotes Newton’s gravitational constant and $\bar{\rho}_{\rm m}$ is the matter density in the FRW background.

2.2 Light cones from simulations

In this section we explain the numerical methods we use to implement the theory discussed above and build mock light cones starting from the output of a simulation. We begin with the calculation of the gravitational potential. Following a standard procedure, we use the particle distribution in each snapshot to compute the matter density contrast on a regular Cartesian grid with the cloud-in-cell method (e.g. Hockney & Eastwood, 1988). We then solve the Poisson equation using a fast Fourier transform and obtain $\phi({\bf x}_{\rm r},t)$ as well as its spatial derivatives (by spectral differentiation). Partial time derivatives of the potential are computed with a finite-difference method that combines several consecutive snapshots at fixed comoving position ${\bf x}_{\rm r}$.

There are a few subtleties at play in the calculation of the galaxy shift and the magnification given in equations (8), (9) and (12). All these quantities include local terms evaluated at a specific position and non-local parts that are expressed as integrals along the line of sight to the observer. The integrals should be taken in redshift-space where the photon path is a straight line. Since $|\phi|\ll 1$ and deflections are generally small, we take the integrals in real space which is correct to linear order in the perturbations as in the Born approximation in quantum mechanics. We use the fast voxel traversal algorithm by Amanatides & Woo (1987) to perform the integrals within the grid over which the gravitational potential is evaluated. All functions appearing in the integrands are interpolated in time (here converted into the line-of-sight distance to make sure that everything is computed on the backward light cone of the observer) such that their first derivatives are continuous. Note that the integration path starts at the observer and ends at a fixed redshift-space position which is unknown for all the simulated galaxies. Although to linear accuracy we could use the real-space position of the galaxies, we implement the following procedure which is slightly more accurate (see also Fig. 3). (i) We evaluate all local terms in equations (8) and (9) and shift the galaxies accordingly. (ii) The non-local terms are estimated with integrals that run from the observer to the position of the galaxy shifted by the local terms. (iii) An additional shift due to the non-local terms is imposed to obtain the final redshift-space position. In principle, steps (ii) and (iii) could be iterated until numerical convergence is achieved. However, this is not necessary in practice since the local terms generate much larger shifts than the non-local ones. This is fortunate because the integration along the line of sight is by far the slowest element of the LIGER code. Magnification is computed along the same lines (similar line-of-sight integrations have been used by White & Hu, 2000; Kiessling et al., 2011; Fosalba et al., 2008, 2015, for weak-lensing studies).

Of course we do not shift all the galaxies at all times as this would significantly slow down the code and also be useless. We first identify the snapshots within which a given galaxy would cross the backward light cone of the observer in the absence of metric perturbations. We then calculate and apply the redshift-space displacements considering a few outputs surrounding this time. Finally, we compute the intersection of the world line of the galaxy with the straight light cone of the observer in redshift space and we save this position and the corresponding magnification.

Each light cone identifies a sub-region in space-time corresponding to a three-dimensional ball in comoving redshift space. Its radius is limited by the box size of the underlying simulation. To avoid replications or spurious correlations due to the periodic boundary conditions applied in cosmological simulations, we limit the radius of the balls to one third of the box size. This way each light cone covers nearly 15 per cent of the simulation volume and we can place five different observers from a single run avoiding intersections.

2.3 Light cones based on dark-matter-only simulations

The LIGER method is general and can be used with all kinds of cosmological simulations. Whenever galaxy positions and luminosities are available (from either a hydrodynamic simulation or a semi-analytic model based on an N-body run), it is straightforward to apply shifts and magnifications at their locations. However, the transverse size of the light cones rapidly increases with redshift so that very large simulation boxes are required to cover wide opening angles. In this case, running simulations with sufficient spatial and mass resolution to follow galaxy formation is computationally challenging. For this reason, LIGER has the option to shift the dark-matter particles themselves and create the galaxy density field a posteriori. The problematic step is to account for galaxy biasing. We illustrate how our implementation works by reasoning in terms of continuous densities.

To linear order in the perturbations, we can write the matter density contrast in redshift space as

[TABLE]

where we have conveniently collected the corrections due to the metric distortions into the term (Yoo, 2009; Yoo, Fitzpatrick & Zaldarriaga, 2009; Bonvin & Durrer, 2011; Challinor & Lewis, 2011; Jeong, Schmidt & Hirata, 2012)

[TABLE]

in which $\phi_{v}$ is the linear velocity potential222This term originates because $\delta_{\rm sim}$ is defined in the synchronous comoving gauge while all the rest is set in the Poisson gauge. at the galaxy position (i.e. $v_{i}=\partial_{i}\phi_{v}$) and the apparent redshift change $\delta\ln a=\delta z/(1+z)$ due to the perturbations is

[TABLE]

To the same accuracy, galaxy clustering in redshift space can be modelled in terms of three redshift-dependent bias parameters, $b$, ${\cal Q}$ and ${\cal E}$, encoding information about different properties of the galaxy population under study, namely (Challinor & Lewis, 2011; Jeong, Schmidt & Hirata, 2012)

[TABLE]

The expression above assumes that the intrinsic perturbation in the galaxy number density is $\delta_{\rm g,r}=b\,\delta_{\rm sim}$ with $b$ the linear bias parameter. It also considers that lensing magnification alters the observed number density of galaxies. This effect is quantified by the magnification-bias parameter

[TABLE]

where $\bar{n}_{\rm g}(>L)$ denotes the comoving number density of galaxies with luminosity larger than $L$ and the derivative is evaluated at the (redshift-dependent) limiting luminosity of the survey.333For simplicity, we assume that the list of targets for spectroscopic observations is flux limited. In case also a size cut is applied, another redshift-dependent function should be added to ${\cal Q}$ since gravitational lensing also alters the size of galaxy images (Schmidt et al., 2009). Finally, equation (19) takes into account that the comoving number density of galaxies in the sample might change with redshift. This phenomenon is described by the ‘evolutionary bias’ parameter

[TABLE]

Our goal is to connect $\delta_{\rm g,s}$ with the particle density in the N-body simulations. For simplicity, we assume that $|{\cal H}\phi_{v}|\ll|\delta\ln a|$ and neglect the velocity potential which could influence galaxy clustering only on scales comparable with the Hubble radius. Then equations (16) and (19) give

[TABLE]

which can be used to derive the local galaxy number density $n_{\rm g,s}=\bar{n}_{\rm g}\,(1+\delta_{\rm g,s})$. By expressing the matter fields in terms of the density of N-body particles, i.e. $\delta_{\rm sim}+1={n_{\rm sim,r}}/{\bar{n}_{\rm sim}}$ and $\delta_{\rm s}+1={n_{\rm sim,s}}/{\bar{n}_{\rm sim}}$, we can write

[TABLE]

where $w_{\rm g}=\bar{n}_{\rm g}/\bar{n}_{\rm sim}$ denotes the mean number of galaxies per simulation particle at a given redshift. The products $w_{\rm g}n_{\rm sim,r}$ and $w_{\rm g}n_{\rm sim,s}$ rescale the unbiased density fluctuations in the simulations to the galaxy mean density. The magnification term in equation (23) reflects the relative change of the galaxy counts per particle which is proportional to $w_{\cal Q}={\cal M}^{\cal Q}$. For $|{\cal M}-1|\ll 1$, we can thus write $w_{\rm g}w_{\cal Q}n_{\rm sim,s}=w_{\rm g}[{\cal Q}({\cal M}-1)+1]n_{\rm sim,s}$, so that

[TABLE]

By using the definition of ${\cal E}$ and linearising, $w_{\rm g}(z)w_{\cal Q}(\hat{\bf n}_{\rm s},z)n_{\rm sim,s}(\hat{\bf n}_{\rm s},z)+\bar{n}_{\rm g}(z){\cal E}\,\delta\ln a$ coincides with $w_{\rm g}(\bar{z})w_{\cal Q}(\hat{\bf n}_{\rm s},z)n_{\rm sim,s}(\hat{\bf n}_{\rm s},\bar{z})$ where $\bar{z}=z-\delta z$ is the redshift in absence of perturbations (note that to first order it is equivalent to evaluate $w_{\cal Q}$ at $z$ or $\bar{z}$). Eventually, making explicit the arguments of all functions, we obtain

[TABLE]

We use this expression to compute $n_{\rm g,s}$ from the simulations (see Fig. 3 for a schematic representation). In practice, we weigh the shifted and unshifted dark-matter particles according to equation (25). Once the light cone for the matter has been constructed, it is very fast to build the galaxy mocks for many different galaxy populations. This corresponds to changing the functions $w_{\rm g}(z)$, $b(z)$ and ${\cal Q}(z)$.

2.3.1 Doppler terms

Isolating the terms proportional to the velocity field in the right-hand side of equation (19) we obtain

[TABLE]

where $\partial/\partial\chi_{\rm s}=n^{i}\partial_{i}$,

[TABLE]

and, assuming a flat $\Lambda$CDM universe,

[TABLE]

Equation (26) coincides with the seminal result for the linear redshift-space distortions derived by Kaiser (1987). In the classical literature, the function $\alpha$ is often written as $\alpha=2+{\partial\ln\bar{n}_{\rm g}}/{\partial\ln\chi_{\rm s}}$ (Kaiser, 1987; Zaroubi & Hoffman, 1996; Hamilton, 1998). Taking into account that we observe galaxies on our past light-cone reveals that several physical effects influence $\alpha$ (see also McDonald, 2009; Yoo, 2009; Bertacca et al., 2012; Raccanelli et al., 2016a). Equations (28), get contributions from geometric distortions, redshift evolution (or redshift-dependent selection effects), Doppler magnification and cosmic acceleration.

Following a standard practice in cosmology, we label the expression proportional to $\alpha$ in equation (26) with the collective name of Doppler terms. Their contribution is usually neglected in clustering studies. In fact, for an ideal galaxy sample with $\alpha\simeq 2$ and if the depth of a galaxy redshift survey is much larger than the comoving wavelength of interest, the Doppler induced $\delta_{\rm g,v}$ is heavily suppressed (due to the $\chi_{\rm s}^{-1}$ scaling) with respect to the signal generated by the radial velocity gradient which is always comparable to density perturbations (Kaiser, 1987). This reasoning relies upon the distant-observer approximation. However, it has been shown that the Doppler corrections can alter the galaxy autocorrelation function at large angular separations in a significant way (Papai & Szapudi, 2008; Raccanelli, Samushia & Percival, 2010). In this case, the two terms on the right-hand side of equation (26) can be of comparable sizes. Moreover, $Q$ and ${\cal E}$ can drive $\alpha$ sensibly away from 2 and thus enhance the chance of detecting the Doppler terms from observational data. We will return to this issue in Section 3.6.

3 Examples

We present two sample applications of LIGER: first, we estimate the importance of magnification bias in an Euclid-like survey and then we investigate the detectability of Doppler terms in a low-redshift galaxy catalogue based on the concept of the Square Kilometre Array (SKA). To begin with, we introduce the numerical simulations and the statistical methods we use. We then describe the specifications for the surveys and discuss our results in Sections 3.4 and 3.6.

3.1 N-body simulations

We run a large number of cosmological N-body simulations using the L-PICOLA code (Howlett, Manera & Percival, 2015) and subsequently apply LIGER to their outputs. L-PICOLA is an implementation of the COLA method (Tassev, Zaldarriaga & Eisenstein, 2013) in which the large-scale dynamics is solved using second-order Lagrangian perturbation theory while a particle-mesh algorithm is used for the small scales. This technique is orders of magnitude faster than standard N-body codes and accurately simulates the clustering of matter on large scales. This makes it an ideal tool to build large mock catalogues for studying galaxy clustering although it does not resolve the internal dynamics of dark-matter haloes.

Our simulations include $1024^{3}$ particles in a periodic cube with side length $L$. In order to cover the relevant volumes, we use very large values of $L$, namely $12h^{-1}$ Gpc for our first application and $5h^{-1}$ Gpc for the second one. This way we obtain 165 light cones extending to redshift 2.3 for the Euclid-like mocks and 125 light cones extending to redshift 0.6 for the more local mocks (an example is shown in Fig. 4). In all cases, the gravitational potential in LIGER is evaluated on a grid with $512^{3}$ cells.

It is worth stressing that LIGER is completely general and can be applied to the output of any N-body code. Here we use L-PICOLA because it is ideal for our purposes. Note that we do not make use of the built-in feature to build light cones on the fly implemented in L-PICOLA. However LIGER could be merged with it in the future.

3.2 Angular power spectra

Our examples focus on large-scale galaxy clustering that we quantify in terms of the angular power spectrum. We first divide our mock light cones into multiple redshift bins and measure the projected galaxy number density contrast on the sky, $\sigma_{\rm g}^{(i)}(\mn@boldsymbol{\theta})$, for each of them (labelled by the index $i$). We then decompose $\sigma_{\rm g}^{(i)}(\mn@boldsymbol{\theta})$ in spherical harmonics, $\sigma_{\rm g}^{(i)}(\mn@boldsymbol{\theta})=a_{lm}^{(i)}\,Y_{lm}(\mn@boldsymbol{\theta})$, with

[TABLE]

and measure the angular auto- and cross-spectra between all redshift bins using

[TABLE]

In practice, we use the Healpix algorithm (Górski et al., 2005) to build digitized maps of $\sigma_{\rm g}^{(i)}(\mn@boldsymbol{\theta})$ from which we calculate the power spectra. In all cases, we make sure that the pixel size does not affect our results in the range of scales of interest. The galaxy density in a pixel is computed from the distribution of the N-body particles as described in Section 2.3. Since our particles are rather massive ($1.2\times 10^{14}\,h^{-1}M_{\odot}$ for the $12\,h^{-1}$ Gpc boxes and $8.8\times 10^{13}\,h^{-1}M_{\odot}$ for the $5\,h^{-1}$ Gpc ones), each of them ‘contains’ multiple galaxies (i.e. $w_{\rm g}\sim 10$). This is unavoidable given the extremely large volumes covered by our simulations and the obvious limitations in computing time and memory usage. Although the resulting overdensity field has increased shot noise with respect to the actual galaxy distribution, our statistical analysis is not influenced by it. In fact, we never attempt to subtract shot noise from the power spectra as our study is based on the comparison of different sets of mock catalogues. Even more importantly, in all cases, shot noise is by far subdominant with respect to the sample and cosmic variance of the clustering signal we are interested in (note that our main results are based on the analysis of cross statistics between galaxy samples at different redshift).

We also take into account that most survey geometries do not cover the full sky and exclude wide regions surrounding the galactic plane. In order to simulate a realistic setting, we mask an appropriate amount of the sky around the equator of the observer and consider two distinct regions around the poles covering a fraction $f_{\rm sky}$ of the celestial sphere. Spherical harmonics are not orthogonal over finite solid angles and the pseudo power spectrum $\tilde{C}_{l}^{(ij)}$ measured from the cut sky does not coincide with $\hat{C}_{l}^{(ij)}$ (Peebles, 1973; Wandelt, Hivon & Górski, 2001). We use the standard method by Hivon et al. (2002) to construct an unbiased estimate of the full-sky spectra $\hat{C}_{l}^{(ij)}$ which is obtained multiplying $\tilde{C}_{l}^{(ij)}$ by the inverse of a mode-mode coupling matrix that depends on the survey mask.

3.3 Statistical analysis

Our sample applications aim at quantifying the detectability of some specific redshift-space effects from measurements of two-point statistics. We assume that we can perfectly model the different contributions to the clustering signal and check whether including or excluding some of them improves or worsen the fit to the mock data including the full physics. In particular, we proceed as follows. We isolate a particular effect (say, e.g., magnification bias) and denote its expected partial contribution to the model galaxy power spectrum as $C_{l}^{\rm(A)}$ so that its complementary part is $C_{l}^{\rm(B)}=C_{l}-C_{l}^{\rm(A)}$ (here the superscripts $(ij)$ indicating the redshift bins are understood to simplify notation). We then fit the power spectra extracted from our mock catalogues with the model $M_{l}=\epsilon C_{l}^{\rm(A)}+C_{l}^{\rm(B)}=C_{l}+(\epsilon-1)C_{l}^{\rm(A)}$ where the coefficient $\epsilon$ can only assume the values zero or one. The question we want to address is to what statistical significance the data favour $\epsilon=1$, i.e. how necessary it is to add $C_{l}^{\rm(A)}$ to the model in order to fit the data $\hat{C}_{l}$ . All this boils down to comparing the quality of the fit obtained using $\epsilon=0$ and $\epsilon=1$. This exercise can be performed following different statistical procedures which give very similar results.

3.3.1 Frequentist approach: simple hypotheses (SH)

We want to test the null hypothesis $H_{0}:\epsilon=0$ against the alternative hypothesis $H_{1}:\epsilon=1$. Let ${\cal L}_{0}$ and ${\cal L}_{1}$ denote the likelihood of the data under $H_{0}$ and $H_{1}$, respectively. Based on the Neyman-Pearson lemma, the likelihood-ratio statistic $\lambda={\cal L}_{0}/{\cal L}_{1}$ provides the most powerful test for two simple hypotheses. If the data do not support $H_{0}$, then the likelihood ratio should be small. Thus, we reject the null hypothesis with confidence level $\tau$, if $\lambda\leq\omega$, where $\omega$ is a constant such that the probability $P(\lambda\leq\omega|H_{0})=\tau$. In order to determine $\tau$ it is thus necessary to determine the probability distribution of the test statistic under $H_{0}$. Working with the log-likelihood $\chi^{2}=-2\ln{\cal L}$, the rejection condition becomes $\Delta\chi^{2}=\chi_{1}^{2}-\chi_{0}^{2}>\ln\omega$. Assuming Gaussian errors for the angular power spectrum with covariance matrix $\Sigma_{lm}=\langle\hat{C}_{l}\hat{C}_{m}\rangle-\langle\hat{C}_{l}\rangle\langle\hat{C}_{m}\rangle$, we find that, under $H_{0}$, $\Delta\chi^{2}$ follows a Gaussian distribution with mean

[TABLE]

and variance $4\mu$ (see Appendix A). Therefore, we reject $H_{0}$ at the 95 per cent confidence level if $\Delta\chi^{2}>\mu+3.29\sqrt{\mu}$. The coefficient 3.29 should be replaced with 4.652 to get a 99 per cent confidence level. A formal $5\sigma$ rejection is obtained for $\Delta\chi^{2}>\mu+10\sqrt{\mu}$.

It is worth noticing that, if the covariance matrix of the measurements does not depend on $\epsilon$, then the expected value of $\lambda$ under $H_{1}$ is $E(\lambda|H_{1})=-E(\lambda|H_{0})=-\mu$. Therefore, the mean values $E(\lambda|H_{1})$ and $E(\lambda|H_{0})$ will be separated by more than ${\cal N}$ standard deviations of the $\lambda$ distribution only if $\sqrt{\mu}>{\cal N}$. This is why $\sqrt{\mu}$ is often denoted as the signal-to-noise ratio, $S/N$, of $C_{l}^{\rm(A)}$. On the other hand, if the covariance depends on $\epsilon$, then $E(\lambda|H_{1})\neq-E(\lambda|H_{0})$ and the dispersions around the mean of $\lambda$ under $H_{0}$ and $H_{1}$ will be different (see Appendix A).

3.3.2 Frequentist approach: composite hypothesis (CH)

We also consider a generalized likelihood-ratio test with a compound alternative hypothesis. In this case we contrast the null hypothesis $H_{0}:\epsilon=0$ with $H_{1}:\epsilon\neq 0$. Let ${\cal L}_{\rm max}$ be the maximum value of the likelihood of the data when $\epsilon$ is varied (between 0 and 1) and ${\cal L}_{0}$ the corresponding likelihood under the null hypothesis. We form the ratio $\lambda={\cal L}_{0}/{\cal L}_{\rm max}$ which is always between 0 and 1. We assess the statistical significance of the test by comparing $\Delta\chi^{2}=-2\ln\lambda$ to the chi-square distribution with one degree of freedom (as we only tuned one parameter to determine ${\cal L}_{\rm max}$). Basically we convert $\Delta\chi^{2}$ into the corresponding percentile of the chi-square distribution. This is the confidence level at which the null hypothesis can be rejected. In this case, 95 (99) per cent confidence corresponds to a critical value of $\Delta\chi^{2}=3.84$ (6.64).

3.3.3 Fisher information

An alternative approach consists of quantifying the Fisher information that the angular power spectrum carries about $\epsilon$ (which, in this case, is assumed to be a real number). Assuming Gaussian errors, we obtain that the Fisher ‘matrix’ for $\epsilon$ is

[TABLE]

(we adopt the ‘field’ perspective as in Carron, 2013). It follows from the Cramèr-Rao inequality that $\mu^{-1/2}$ gives a lower bound for the expected uncertainty on $\epsilon$ (i.e. the mean curvature of the likelihood function ${\cal L}(\epsilon)$ at its peak). Therefore, $\sqrt{\mu}$ corresponds to the signal-to-noise ratio with which $\epsilon$ can be measured.

3.3.4 Bayesian approach (BA)

Assuming a flat prior for $\epsilon$, the posterior probabilities for $\epsilon=0$ and $\epsilon=1$ are $P_{0}=\frac{{\cal L}_{0}}{{\cal L}_{0}+{\cal L}_{1}}$ and $P_{1}=\frac{{\cal L}_{1}}{{\cal L}_{0}+{\cal L}_{1}}$, respectively. The Bayes factor, $K=P_{0}/P_{1}$, thus corresponds to the likelihood ratio ${\cal L}_{0}/{\cal L}_{1}$. According to the Jeffrey scale, there is strong (decisive) evidence against $H_{0}$ if $K>10$ (100) which gives $\Delta\chi^{2}>4.605$ (9.21). Using the alternative scale by Kass and Raftery gives very strong evidence for $\Delta\chi^{2}>10$.

3.3.5 Likelihood estimation

All the statistical methods listed above require the calculation of the likelihood function for $\epsilon$ in each mock realisation. To do this, we need a model and the covariance matrix for the power spectra (we assume Gaussian measurement errors). We build an ‘exact’ and unbiased model by averaging the angular power spectra obtained for all our mock light cones. Since there are no particularly deviant realizations, the average spectra are smooth. In parallel, we use the maximum-likelihood estimator to get a first approximation for the covariance matrix, $\hat{\Sigma}_{lm}$. It is well known that the precision matrix obtained by inverting $\hat{\Sigma}_{lm}$ is not very accurate. Although we make sure to consider enough mock skies so that the covariance matrix of our data vector is invertible, it still contains considerable noise. We thus use the shrinkage method (Schäfer & Strimmer, 2005) to reduce the noise in $\hat{\Sigma}_{lm}$. As a target we use a diagonal matrix which is always compatible with our estimates.

3.4 Magnification bias in a Euclid-like survey

As a simple application of the LIGER method, we discuss the detectability of magnification bias in a Euclid-like survey. Related work has been presented by Di Dio et al. (2014) and Montanari & Durrer (2015). These authors focused on the Fisher information matrix as a forecasting tool, while we base our study on the statistical analysis of a large number of mock catalogues.

3.4.1 Euclid spectroscopic sample

Euclid is a medium-class mission of the European Space Agency planned for launch in 2020. It will map the distribution of star-forming galaxies through their redshifted H$\alpha$ emission in the regions with galactic latitude larger than 30 degrees ($f_{\rm sky}=0.36$). Low-resolution spectroscopy in the near infrared will be used to measure galaxy redshifts in the range $0.7<z<2.0$.

The specifics of the Euclid redshift survey depend on the poorly known properties of emission-line galaxies at moderate redshifts. In order to calculate the redshift distribution of the galaxies as well as $Q(z)$ and ${\cal E}(z)$ we use the redshift-dependent luminosity function by Pozzetti et al. (2016, model two) and assume a limiting line flux of $3.0\times 10^{-16}$ erg cm*-2* s*-1*. Further, we use the linear fit for the galaxy bias $b(z)$ given in Pozzetti et al. (2016). Our results are summarised in Fig. 5.

3.4.2 Angular power spectra

We split each of the 165 Euclid mocks into four redshift bins (with boundaries $0.7,0.86,1.06,1.35,2.0$) that on average contain the same number of galaxies ($\sim 1.6\times 10^{7}$). We then compute the auto- and cross-power spectra among all bins. To reduce noise we do not resolve individual multipoles and use ten bandpowers with $\Delta l=20$. Fig. 6 shows $C_{l}^{(44)}$ (left-hand panel) and $C_{l}^{(14)}$ (right-hand panel), where the indices 1 and 4 indicate the lowest and highest redshift bin, respectively. The blue curves correspond to the mean power spectra extracted from the mock catalogues that include all redshift-space effects (hereafter GR). Note that the auto-spectrum is more than ten times larger than the cross-spectrum. The shaded regions indicate the standard deviation of the spectra over the 165 realisations in the Euclid-like (cross pattern) and in the full-sky mocks (light blue).

It is interesting to analyse the different contributions to the spectra. The most commonly considered source of redshift-space distortions is the so-called Kaiser effect due to the gradient of the galaxy radial peculiar velocities, i.e. the first term on the right-hand side in equation (26). In our approach, this correction derives from the $(n_{\rm s}^{i}v_{i})_{\rm e}$ terms in the particle shift and is always mixed with the Doppler contribution. In order to evaluate the relative importance of the velocity-induced shift, we build a new set of Euclid mock catalogues (based on the same N-body simulations as the GR ones) in which we replace equations (8), (9) and (12) with $\delta\chi=-(n_{\rm s}^{i}v_{i})_{\rm e}/\mathcal{H}$, $\delta x^{i}=0$ and ${\cal M}=1$, respectively (this is the standard way to implement redshift-space distortions in simulations and omits the terms proportional to ${\cal Q}$ in $\alpha$). We dub these light cones KD, a short for ‘Kaiser and Doppler’. The red curves in Fig. 6 show the mean clustering signal extracted from the KD mocks (which from now on we refer to as the KD model for the auto- and cross-spectra). The lack of power with respect to the GR results is evident: $C_{l}^{(44)}$ is underestimated by $\sim 3$ per cent and $C_{l}^{(14)}$ oscillates around zero. In the left-hand panel of Fig.  7, we demonstrate that these differences are highly significant. The histogram on the left-hand side displays the distribution of $\Delta\chi^{2}=\chi^{2}_{\rm GR}-\chi^{2}_{\rm KD}$ obtained fitting all the $C_{l}^{(ij)}$ from the 165 GR mocks with the GR and the KD models, respectively. Similarly, the histogram on the right-hand side shows the corresponding distribution for the KD mocks. Based on the fact that the histograms are well separated, we conclude that an Euclid-like survey should be able to detect the signature of redshift-space distortions that are not included in the KD model. To better quantify how inaccurately the KD model fits the mock GR data, we apply the statistical tests we have introduced in Section 3.3. The cumulative distribution over the 165 GR mocks of the significance with which the KD model is rejected is shown in the right-hand panel of Fig. 7. Typically the data disfavour the simpler model with 8$\sigma$ confidence or a Bayes factor of $10^{13}$. The precise statistical significance of this result is very sensitive to the assumed ${\cal Q}(z)$. For instance, it increases to 14$\sigma$ if we use model 3 from Pozzetti et al. (2016).

This result is not surprising. It is well known that weak gravitational lensing alters the observed clustering signal in deep magnitude-limited surveys (Turner, 1980). The influence of lensing is twofold: (i) the actual magnitude limit of the survey fluctuates on the sky and with redshift; (ii) the surface density of galaxies on the sky (and thus their volumetric density in redshift space) is changed. The last two terms in equation (17) and the $Q$-dependent term in equation (19) summarise the net effect on the galaxy overdensity field. Following some early detections (Bartelmann & Schneider, 1994; Norman & Williams, 2000), the weak-lensing effect on clustering has been measured with high statistical significance ($8\sigma$) by cross correlating samples of distant quasars and background galaxies in the Sloan Digital Sky Survey (Scranton et al., 2005). Given this premise, we build a third set of mock light cones (labelled $\kappa$KD) in which we account for the redshift-space distortions due to both the peculiar velocities and weak lensing assuming that the convergence is the only source of magnification, i.e. ${\cal M}=1+2\kappa$. To account for both magnification bias and the volume corrections due to lensing, we simply weigh the N-body particles proportionally to $[{\cal M}(\hat{\bf n}_{\rm s},z)]^{{\cal Q}(z)-1}$ (instead of the standard $[{\cal M}(\hat{\bf n}_{\rm s},z)]^{{\cal Q}(z)}$) and use the same shifts as in the KD mocks.444 The additional ${\cal M}^{-1}$ results from a volume distortion due to lensing. To test the consistency of our code, we evaluate this effect in two ways. In the GR mocks the change of volume is realized by the particle shift. Alternatively, we weigh particles proportionally to ${\cal M}^{-1}$ in the KD mocks. We find the same result. The resulting spectra (yellow lines in Fig. 6 which are barely distinguishable by eye from the blue ones) provide an excellent fit to the $C_{l}^{(ij)}$ derived from the GR mocks. This suggests that the measurable differences between the full signal and the KD model are due to gravitational-lensing convergence.

In order to make a quantitative analysis and investigate whether other (more subtle and interesting) light-cone effects (e.g. Doppler lensing and potential terms) might be detectable with a Euclid-like survey, we once again resort to statistics. Fig. 8 shows that, in almost all mock GR realizations, the null hypothesis that the data are generated under the $\kappa$KD model cannot be rejected to any meaningful confidence level. Although only few of the skies presents deviations larger than $3\sigma$, Fig. 8 indicates that the GR mocks contain an additional signal (most likely due to Doppler lensing) which is however comparable than the noise. We thus conclude that no additional sources of redshift-space distortions beyond the Kaiser effect, Doppler contributions and magnification bias can be detected from the angular clustering of all galaxies with spectroscopic redshifts in a Euclid-like survey divided in four equally-populated redshift bins. This, of course, does not prevent the development of dedicated probes to isolate additional contributions, in particular combining photometric and spectroscopic data to define multiple tracers of the large-scale structure along the lines of the forthcoming discussion in Section 3.6.

3.5 Estimating covariances with LIGER

Mock catalogues provide a direct way to estimate the covariance matrix of observables and test the range of validity of idealised theoretical models for the statistical errors and their correlations. The covariance matrix of the power spectrum on large scales is often described in terms of the Gaussian approximation,

[TABLE]

where $\delta_{lm}$ denotes the Kronecker symbol and we have restored the superscripts for the redshift bins, for clarity. As a byproduct of the study presented in the previous section, we use the Euclid mock catalogues to assess the accuracy of the Gaussian approximation. Overall, our numerical estimates are compatible with the diagonal structure of $\Sigma_{lm}^{(ik)(jh)}$. The diagonal elements obtained from the full-sky mocks are in excellent agreement with equation (33). However, we find that the expected scaling with respect to $f_{\rm sky}$ holds true only for low multipoles. Fig. 9 shows how the average standard deviation of $C_{l}^{(ij)}$ in the mocks varies as a function of $f_{\rm sky}$. The thick blue line corresponds to the multipole interval $20<l\leq 40$, while the dashed line represents the bin $140<l\leq 160$. Both curves have been normalised by the corresponding standard deviations measured in the full-sky mocks. The hatched area indicates the error on the ratio estimated bootstrapping the realisations. The black line highlights the theoretical scaling proportional to $f_{\rm sky}^{-1/2}$. Note that the statistical error for $C_{l}^{(ij)}$ at $l\sim 30$ follows this curve for $f_{\rm sky}>0.35$ but rapidly departs from it for smaller sky fractions. Basically, the data cannot optimally constrain the large-scale power when the footprint of the survey covers too small a fraction of the sky. Higher multipoles deviate from the ideal relation for even larger values of $f_{\rm sky}$. For an Euclid-like survey, the statistical error on $C_{l}^{(ij)}$ at $l\sim 150$ is on average 17 per cent larger than expected using the $f_{\rm sky}^{-1/2}$ scaling. All this exemplifies the usefulness of LIGER (and mock catalogues in general) to estimate the size of measurement errors in clustering statistics and warns against using simple approximations outside the range within which they have been accurately tested.

3.6 Detectability of Doppler terms at low redshift

The question we want to address in this section is whether forthcoming probes of the large-scale structure will be able to provide evidence for additional sources of redshift-space distortions beyond the classic Kaiser effect and magnification bias. Although the term proportional to $\alpha$ in equation (26) has been almost invariably neglected in past studies of galaxy clustering, it might become detectable using data from future wide-angle surveys. Recent work based on analytical calculations has concluded that the Doppler terms should leave a measurable imprint on the cross correlations between two galaxy populations (McDonald, 2009; Yoo et al., 2012; Bonvin, Hui & Gaztañaga, 2014) and, possibly, also on the angular clustering of a single population (Raccanelli et al., 2016a). Here we re-examine this case using LIGER to build mock catalogues for the SKA2 survey.

3.6.1 Simulating a galaxy redshift survey with the SKA2

The SKA is an unprecedentedly large and powerful array of radio telescopes that will be built in Australia and South Africa by an international collaboration. The construction will be split into two phases: 10 per cent of the collective area should be in place by 2023 (SKA1) while the full array should follow by 2030 (SKA2). The ‘billion galaxy survey’ conducted with the SKA2 will measure $\sim 10^{9}$ individual galaxy redshifts over 30,000 deg2 using the 21-cm line emission from neutral atomic hydrogen (Maartens et al., 2015). In many senses, this will be the ultimate ground-based redshift survey for cosmology.

We build 125 mock light-cones with the expected characteristics of the billion-galaxy survey. Proceeding as described in Sections 2.3 and 3.1, we populate our N-body simulations with galaxies making sure to reproduce the forecasts for the galaxy number counts presented in Yahya et al. (2015). At low redshifts ($z\lesssim 0.2$) and for flux555Although ${\rm Jy}$ is a unit of flux density for simplicity we refer to it as a flux. limits below $10\,\mu{\rm Jy}$, their predictions for $\bar{n}_{\rm g}$ depend very little on the limiting flux of the survey (implying that $Q\simeq 0$). This might possibly reflect a shortcoming of their fitting formula (which covers a broad redshift range), an imperfection in the HI modeling, or the finite mass resolution of the N-body simulations used to make the forecasts (i.e. fainter galaxies might reside within unresolved dark-matter haloes). Indeed, Yahya et al. (2015) note that observed HI mass functions at low redshifts contain more low-mass objects than found in the forecasts for SKA2. To play safe, we avoid this region of parameter space. We thus use a conservative flux limit of $23\,\mu{\rm Jy}$ to define our main sample (Yahya et al., 2015, use $5.3\,\mu{\rm Jy}$ for their realistic forecasts and $23\,\mu{\rm Jy}$ for the pessimistic ones) and also consider a second galaxy population (the bright sample) with fluxes above $60\,\mu{\rm Jy}$. We derive the corresponding values for $Q(z)$ and ${\cal E}(z)$ after fitting the cumulative number density presented in Yahya et al. (2015) with a third-order polynomial. This is the same approach followed by Camera, Santos & Maartens (2015) and yields consistent results with their revised analysis666In a flux limited survey, $Q$ cannot assume negative values as it would imply that less galaxies are detected when a fainter limiting flux is considered (the presence of an additional bright cut or another selection criterion are necessary to drive $Q<0$). Due to an unfortunate mishap, the original fitting functions presented in Camera, Santos & Maartens (2015) yield $Q<0$ for faint galaxies at low redshifts. Although this mistake does not affect their conclusions, the negative magnification biases have been used by many authors to make forecasts for the SKA2 survey (Camera, Santos & Maartens, 2017). The outcome of our calculations is shown in Fig. 5 together with the functions $b(z)$ which we take directly from Yahya et al. (2015).

To test the impact of the Doppler terms on the clustering of SKA2 galaxies, we build and contrast two sets of mock catalogues. The first includes relativistic effects while the second drops the Doppler terms that are proportional to ${\cal E}$ and $Q$. We omit Doppler magnification by simply ignoring the velocity-dependent terms in equation (12). Further, we neglect the term proportional to ${\cal E}$ in equation (22) so that the weight $w_{\rm g}(\bar{z})$ in the second line of equation (25) is replaced with $w_{\rm g}(z)$. We use the label DS (Doppler suppressed) to indicate the light cones constructed in this way, since it is impossible to isolate the remaining Doppler effects.

3.6.2 Cross spectra and results

General relativistic corrections alter galaxy clustering on large scales with respect to the predictions of the ‘standard model’ including a linear bias and the Kaiser distortions. In particular, they break the symmetry of two-point statistics under the exchange of particles in the pairs. In the distant-observer approximation, the relativistic effects generate odd multipoles in the redshift-space cross-correlation function between two galaxy populations or, equivalently, an imaginary part in the cross spectrum (McDonald, 2009; Yoo et al., 2012; Croft, 2013; Bonvin, Hui & Gaztañaga, 2014; Bonvin, 2014; Raccanelli et al., 2014). In terms of the comoving wavenumber of the perturbations $k$, relativistic corrections to the cross spectra due to Doppler effect and gravitational redshifts are suppressed by a factor ${\cal H}/k$ with respect to the leading standard-model terms (McDonald, 2009). Additional corrections (due to the gravitational potential, the integrated Sachs-Wolfe effect and the Shapiro time delay) are instead suppressed by a factor $({\cal H}/k)^{2}$. Measurements of galaxy clustering on scales comparable with the Hubble radius are therefore necessary to detect them. There is a complication, however. In the standard model, the time evolution of the galaxy populations and wide-angle effects due to the fact we observe on our past light-cone give rise to several anti-symmetric terms with similar amplitudes to the relativistic corrections (Bonvin, Hui & Gaztañaga, 2014). Moreover, dust extinction can further introduce spurious anti-symmetric terms in the galaxy correlation functions (Fang et al., 2011). Finally, from the observational point of view, it is challenging to keep the photometry stable and measure weak clustering signals on very large scales. Therefore, it is still an open question whether the Doppler contribution can be seen. In the rest of this section, we will employ our SKA2 mock light cones to address this issue. In doing this, we will neglect systematic effects due to observational limitations and dust and focus on the feasibility of the experiment from a theoretical point of view.

In order to maximise the chances for a successful outcome, it is important to carefully configure the test we want to perform. There are a few facts to take into consideration: (i) linear peculiar velocities grow bigger at lower redshifts; (ii) $\gamma_{0}$ in $\delta_{\rm g,v}$ is divided by ${\cal H}\chi_{\rm s}$; (iii) the galaxy number density of the samples rapidly decreases for $z\gtrsim 0.2$ thus producing a large evolutionary bias but also increasing noise; (iv) we need to cover enough comoving volume to reduce sample variance. Given all this, we end up considering the interval $0.15<z<0.25$ which we further divide into the bins I: $0.15<z<0.2$ and II: $0.2<z<0.25$. We also make use of the two galaxy populations introduced in Section 3.6.1.

In order to access the scales that are sensitive to the Doppler contribution, we focus on galaxy pairs with very wide angular separations. We thus compute the cross angular power spectrum, $\hat{C}_{l}^{\rm(T_{I}\,B_{II})}$, between the total galaxy sample within the lower redshift bin and the bright subsample within the highest redshift bin. Similarly, we measure $\hat{C}_{l}^{\rm(B_{I}\,T_{II})}$ by considering the bright galaxies in bin I and the full population in bin II. Finally, we consider the difference between the cross spectra:

[TABLE]

Average results over the SKA2 mock light cones are shown in Fig. 10. The left-hand panel displays $\hat{C}_{l}^{\rm(T_{I}\,B_{II})}$ (solid thick line) and $\hat{C}_{l}^{\rm(B_{I}\,T_{II})}$ (dashed thick line) both sampled in bandpowers with $\Delta l=4$. The same quantities derived from the DS mock catalogues are shown in the right-hand panel. The additional Doppler terms included in the GR light cones clearly enhance the difference between the cross spectra at low $l$. They boost $\hat{C}_{l}^{\rm(B_{I}\,T_{II})}$ (i.e. make it more negative) and suppress $\hat{C}_{l}^{\rm(T_{I}\,B_{II})}$. This can be explained as follows. The leading Doppler contribution to the cross spectra on large scales originates from correlating the term $b\delta_{\rm sim}-{\cal H}^{-1}\partial n_{\rm s}^{i}v_{i}/\partial\chi_{\rm s}$ with the Doppler term in equation (26). Therefore, schematically, in terms of single Fourier modes,

[TABLE]

where $f=\mathrm{d}\ln D/\mathrm{d}\ln a$ (with $D$ the linear growth factor of matter perturbations) describes the evolution of the velocity field, while $\mu_{\rm I}$ and $\mu_{\rm II}$ are the cosines between the wavenumber and the two lines of sight in the wide-angle configuration. A peculiarity of our samples is that $b_{\rm B}(z)\simeq b_{\rm T}(z)=b(z)$ to very good accuracy (this is different from Bonvin, Hui & Gaztañaga, 2014, which relies on different linear bias parameters and use three dimensional correlation functions) so that

[TABLE]

Now, the DS mocks only include the geometric distortions and those generated by cosmic acceleration, i.e. $\alpha_{\rm DS}=2+[1-(3/2)\Omega_{\rm m}(z)]\mathcal{H}\chi_{\rm s}$ which is independent of the galaxy sample and only depends on the underlying cosmological model. Thus, $\Delta\hat{C}_{l}\simeq 0$ in the DS case. On the other hand, for the GR mocks we get $\alpha_{\rm B_{I}}=-1.49$, $\alpha_{\rm B_{II}}=-2.42$, $\alpha_{\rm T_{I}}=0.35$ and $\alpha_{\rm T_{II}}=-0.38$ which combine to produce a positive signal. In this case ${\cal Q}$ and ${\cal E}$ give roughly equal contributions.

The shaded areas in Fig. 10 indicate the standard deviation for $\hat{C}_{l}^{\rm(T_{\rm I}\,B_{\rm II})}$ (the scatter for $\hat{C}_{l}^{\rm(B_{\rm I}\,T_{\rm II})}$ is of comparable size). It is evident that the cross spectra extracted from the galaxy survey will be very noisy. In Fig. 11 we show the mean and scatter of $\Delta\hat{C}_{l}$ for the full GR mocks (blue) and the DS ones (red). While the relative error on the single cross-spectra is very large, $\Delta\hat{C}_{l}$ can be measured with much better accuracy (especially for $l<25$). Since both galaxy populations trace the same large-scale structure, most of the noise in the cross-spectra is correlated and thus does not appear in the difference. This exemplifies the advantage of using a multi-tracer approach (McDonald & Seljak, 2009).

We are now ready to investigate whether Doppler effects will be measurable with the SKA2. We first measure $\Delta\hat{C}_{l}$ from the 125 GR mocks and then fit the multipoles in the range $1\leq l\leq 25$ using both the GR and the DS models (we recall that the models are obtained averaging the signal over all the mock light cones). The left-hand side histogram in the left panel of Fig. 12 shows the distribution of the corresponding value for $\Delta\chi^{2}=\chi^{2}_{\rm GR}-\chi^{2}_{\rm DS}$. In parallel, we fit the $\Delta\hat{C}_{l}$ measurements extracted from the DS mocks and plot the corresponding histogram of $\Delta\chi^{2}$ on the right-hand side of the figure. The fact that the two histograms are well separated demonstrates that the SKA2 data should contain enough information to distinguish between the two models. To better quantify this, for each mock light cone we compute the statistical significance with which we can reject the null hypothesis that the GR data are generated by the DS model. The resulting cumulative distribution is plotted in the right panel of Fig. 12 using three different statistical tests (see Section 3.3). In the vast majority of the mock catalogues, the DS model is ruled out at a confidence ranging between 4 and 7$\sigma$ (or with a Bayes factor above 1000 corresponding to decisive evidence). We conclude that the SKA2 should be able to detect the non-standard Doppler contribution to galaxy clustering.

Our results heavily rely upon the multi-tracer technique for the suppression of the statistical noise. After repeating the analysis with a single tracer (using cross correlations between adjacent and narrow redshift bins), we find that only extreme values for the magnification and evolutionary biases (e.g. $Q\gtrsim 10$ and $|{\cal E}|\gtrsim 20$) could lead to a statistically significant detection.

4 Summary & conclusions

There are multiple reasons for which galaxy clustering requires a proper general relativistic description. (i) We observe events lying on our past light cone. (ii) The propagation of light is affected by the presence of inhomogeneities in the matter distribution. (iii) In consequence, galaxy observables (i.e. redshift, flux in some waveband and angular position on the sky) are influenced by the large-scale structure intervening between the source and the observer. However, when we interpret observations we use an unperturbed FRW model to translate redshifts and fluxes into distances and absolute luminosities. This leads to redshift-space distortions, i.e. the reconstructed galaxy density does not coincide with the actual one. The most important source of the discrepancy is the correction due to the peculiar velocity gradient (Kaiser, 1987) but it is long known that there are additional contributions and that they might become significant at large angular separations. Although no consensus has been reached yet, some recent studies based on analytical calculations and on the Fisher information matrix have concluded that signatures of these additional corrections should be detectable with the next generation of wide-angle surveys.

This paper describes the LIGER method, a numerical technique to build mock galaxy catalogues including all general relativistic corrections at linear order in the cosmological perturbations. LIGER post processes the output of a Newtonian simulation and combines its snapshots at constant background time to build the galaxy distribution in comoving redshift space. The LIGER method is composed of three steps: (i) we shift the position of the simulated galaxies from real to redshift space; (ii) we evaluate the corresponding magnification due to gravitational lensing; (iii) we find the intersection of the modified world lines of the galaxies with the unperturbed backward light cone of the observer in redshift space. Steps (i) and (ii) include both local corrections and terms that have been integrated along the line of sight using the Born approximation. Note that standard mock galaxy catalogues generally include only the Kaiser effect for step (i) and do not account for lensing (with the exception of Guimarães, Myers & Shanks, 2005) although magnification maps are routinely built for weak-lensing studies (e.g. Wambsganss, Cen & Ostriker, 1998; Jain, Seljak & White, 2000; White & Hu, 2000; Vale & White, 2003; Fosalba et al., 2008, 2015; Hilbert et al., 2009; Kiessling et al., 2011).

LIGER is characterized by the following features. (i) It can be applied to the output of any Newtonian simulation (N-body or hydrodynamic) independent of the code with which it has been run. (ii) It is fast to execute so that it is computationally feasible to produce very large numbers of mock catalogues for a given survey. (iii) A variant of the standard implementation has been specially designed to work with simulations that cover very large comoving volumes but do not resolve single galaxies. In this case, the galaxy density field is obtained by biasing the dark-matter distribution.

The main applications of LIGER are foreseen for forthcoming wide-angle spectroscopic surveys like DESI777http://desi.lbl.gov, Euclid888http://www.euclid-ec.org, HETDEX999http://hetdex.org, SKAII101010https://www.skatelescope.org, SPHEREx111111http://spherex.caltech.edu, SuMiRe121212http://sumire.ipmu.jp/en/ as well as photometric surveys like J-PAS131313http://j-pas.org and LSST141414https://www.lsst.org. Mock catalogues will be used as forecasting tools and to measure biases and covariance matrices of estimators for several statistics of the large-scale structure. As a proof of concept, we quantify the impact of magnification bias in the angular clustering of galaxies for a Euclid-like survey. Our results show that lensing convergence generates a non-negligible correction for the angular power spectra of galaxies in broad redshift bins and dominates cross-spectra between galaxies at widely separated redshifts (similar conclusions have been reached by Di Dio et al., 2014; Montanari & Durrer, 2015, using analytical calculations). The convergence signal can be detected at 8$\sigma$ significance and this provides the intriguing possibility to measure the lensing potential from the cross spectra (Montanari & Durrer, 2015).

Additional redshift-space distortions on top of the standard Kaiser correction and the weak-lensing convergence are generally small. Therefore, customised techniques need be developed in order to measure their signatures. It is foreseeable that LIGER-based mocks will be key to optimising the design of these probes. In anticipation of these future applications, we have investigated the detectability of several additional Doppler terms using the specifics of the planned ‘billion galaxy survey’ with the SKA2 telescope. Our results show that using two galaxy populations with different flux cuts (T: $f>23\,\mu{\rm Jy}$ and B: $f>60\,\mu{\rm Jy}$) and two consecutive redshift intervals (I: $0.15<z<0.2$ and II: $0.2<z<0.25$), it will be possible to measure a significant Doppler-induced signal. The statistic we use is the difference of the angular cross power spectra $\Delta C_{l}=C_{l}^{\rm(T_{I}\,B_{II})}-C_{l}^{\rm(B_{I}\,T_{II})}$. For multipoles $l\leq 25$, this quantity is dominated by the contribution of the Doppler terms and shows a strikingly reduced variance compared to each of the cross correlations due to the fact that both galaxy populations are biased tracers of the same underlying matter density. Based on our simulations, $\Delta C_{l}$ should be detectable with a signal-to-noise ratio of $\sim 5.5$.

Acknowledgements

We thank Alvise Raccanelli and Stefano Camera for discussions. We acknowledge financial support by the Deutsche Forschungsgemeinschaft through the Transregio 33 ‘The Dark Universe’. During the preparation of this work MB was also partially supported by the Bonn-Cologne Graduate School for Physics and Astronomy.

Appendix A Likelihood-ratio test

Let ${\mathbf{x}}$ be a Gaussian data vector with mean ${\mathbf{m}}$ and covariance matrix $\boldsymbol{\mathsf{\Sigma}}$. Let ${\mathbf{M}}$ be a perfect theoretical model for ${\mathbf{m}}$, i.e. ${\mathbf{M}}={\mathbf{m}}$. The likelihood of ${\mathbf{M}}$ given the data is proportional to the probability of observing the data under the hypothesis that ${\mathbf{M}}$ is true, i.e. ${\cal L}({\mathbf{M}}|{\mathbf{x}})\propto\exp(-\chi^{2}/2)$ with $\chi^{2}=({\mathbf{x}}-{\mathbf{M}})^{\rm T}\cdot\boldsymbol{\mathsf{\Sigma}}^{-1}\cdot({\mathbf{x}}-{\mathbf{M}})$. Let also consider an imperfect or incomplete model ${\mathbf{N}}$ such that ${\mathbf{m}}-{\mathbf{N}}={\mathbf{B}}$. The likelihood ratio between ${\mathbf{M}}$ and ${\mathbf{N}}$ is ${\cal L}({\mathbf{M}}|{\mathbf{x}})/{\cal L}({\mathbf{N}}|{\mathbf{x}})=\exp(-\Delta\chi^{2}/2)$ with

[TABLE]

where we have decomposed the data vector in the signal and noise components, ${\mathbf{x}}={\mathbf{m}}+{\mathbf{e}}={\mathbf{M}}+{\mathbf{e}}$, and used the symmetry of the covariance matrix. The fact that the difference in the log-likelihoods depends linearly on ${\mathbf{e}}$ implies that $\Delta\chi^{2}$ follows a Gaussian distribution over an ensemble of realisations of the data vector. Its expectation is

[TABLE]

Similarly, its variance is

[TABLE]

In a classic likelihood-ratio test for a simple hypothesis, the probability distribution of $\Delta\chi^{2}$ obtained under a model (i.e. assuming that this model perfectly describes ${\mathbf{m}}$) is compared against the corresponding distribution obtained under an alternative model. If the covariance matrix of the data does not depend on the model, the two $\Delta\chi^{2}$ distributions differ only in the sign of their mean values (in fact models $M$ and $N$ are switched in the alternative hypothesis and the sign of $\Delta\chi^{2}$ is reversed). The test rejects one of the models if the two distributions of $\Delta\chi^{2}$ are clearly separated with respect to their intrinsic dispersion. Since the distance between the averages is $2\mu$ and the RMS value of each distribution is $2\sqrt{\mu}$, the ratio $2\mu/(2\sqrt{\mu})=\sqrt{\mu}$ is commonly referred to as the signal-to-noise ratio $S/N$. If, however, the covariance matrix is model dependent, then the comparison should be done between two Gaussian distributions with different mean values $\mu_{M}={\mathbf{B}}^{\rm T}\cdot\boldsymbol{\mathsf{\Sigma}}_{M}^{-1}\cdot{\mathbf{B}}$ and $\mu_{N}=-{\mathbf{B}}^{\rm T}\cdot\boldsymbol{\mathsf{\Sigma}}_{N}^{-1}\cdot{\mathbf{B}}$ as well as variances $4\mu_{M}$ and $4\mu_{N}$, respectively.