© M.-A. Breton et al. 2022

1. Introduction

The observation of distant galaxies provides a wealth of information regarding the nature and content of our Universe. Since galaxies trace the underlying matter density field, their spatial distribution can be used to probe the evolution of the large-scale structure. Moreover, the Doppler effect induced by the peculiar motion of galaxies, also known as redshift-space distortions (RSD; Kaiser 1987), leaves a distinct imprint on the three-dimensional clustering of these distant galaxies, which can in turn be used to gain information on the linear growth rate of structure. The latter is directly sensitive to the theory of gravity (Guzzo et al. 2008).

Galaxy clustering has been used for more than two decades to place constraints on the fσ₈ parameter, that is, the growth rate multiplied by the power spectrum normalisation. No significant deviations from General Relativity were found (e.g., Peacock et al. 2001; Tegmark et al. 2006; Blake et al. 2012; de la Torre et al. 2013; Alam et al. 2017, 2021). Future surveys such as Euclid (Laureijs et al. 2011), the Subaru Prime Focus Spectrograph (Takada et al. 2014), the Square Kilometre Array (Square Kilometre Array Cosmology Science Working Group 2020), or the recently started DESI (DESI Collaboration 2016) will probe the distribution of galaxies with unprecedented accuracy and at higher redshifts (z > 1) than before. In particular, Euclid will provide a 1 − 3% level of precision on fσ₈ between z = 0.9 and z = 1.8 (Laureijs et al. 2011; Majerotto et al. 2012).

Because future galaxy surveys will probe higher redshifts with high accuracy, the theoretical modelling of galaxy clustering in the form of summary statistics such as the two-point correlation function will need to be improved to provide an unbiased estimate of cosmological parameters. More precisely, the current modelling of the observed two-point correlation function (or power spectrum) only accounts for peculiar velocities (RSD), while analytical works have shown in linear theory that the next dominant term at high redshift will be due to gravitational lensing (Yoo et al. 2009; Challinor & Lewis 2011; Bonvin & Durrer 2011), also known as magnification bias (Schneider et al. 1992). Gravitational lensing indeed modifies the apparent angular position of sources due to the deflection of light rays along their path, and it also magnifies their fluxes, so that in magnitude-limited surveys, we can find more or fewer objects in magnified regions, depending on the slope of the luminosity function.

It is interesting to remark that in galaxy-galaxy lensing, the theoretical framework required to incorporate the magnification bias corrections was first devised in Ziour & Hui (2008), but this gained more attention only very recently (Duncan et al. 2014; Ghosh et al. 2018; Thiele et al. 2020; Unruh et al. 2020; Joachimi et al. 2021; Lee et al. 2022). The difference in theoretical modelling between Dark Energy Survey one-year (Prat et al. 2018) and three-year (Abbott et al. 2022) galaxy-galaxy lensing analyses, where only the latter accounts for the effect of lensing magnification, is particularly telling. Similarly, several analytical works have investigated the effect of magnification bias on the three-dimensional two-point correlation function (Matsubara 2000; Hui et al. 2007, 2008), but has not been implemented in observational analysis so far. The main reason for this is that, until now, galaxy redshift surveys have probed low enough redshifts, where gravitational lensing does not have an important impact on the overall three-dimensional distribution of galaxies. However, in the near future, this may no longer be the case. Of particular relevance for the present paper, Jelic-Cizmek et al. (2021) investigated the effect of lensing magnification in spectroscopic surveys and found that neglecting magnification bias leads to an underestimation of the growth rate of structure. This work was restricted to pure analytical models in linear theory and in a Fisher analysis. From a numerical point of view, Elkhashab et al. (2022) assessed the detectability of magnification bias on the power spectrum monopole only for a Euclid-like survey.

The present paper aims at investigating the impact of magnification bias on the determination of the growth rate from clustering in redshift space in a typical observational set up. We use a suite of N-body simulations that accounts for the fully non-linear structure formation and perform realistic galaxy clustering analyses, similarly as in observations, in different regimes of magnification bias. We investigate a minimal model to account for the magnification effect on the multipole moments of the redshift-space correlation function, and we study the accuracy with which the growth rate of structure parameter can be recovered.

The paper is organised as follows. In Sect. 2 we discuss the theoretical modelling of the magnified correlation function multipole moments in redshift space. The numerical data and method we used to perform a likelihood analysis are presented in Sect. 3. Finally, the results are shown in Sect. 4, and we conclude in Sect. 5.

2. Theory

The observed spatial distribution of galaxies in the Universe depends on two main aspects: first, galaxies are biased tracers of the matter density field, which means that in ‘real space’ (i.e. the Universe as it is really), galaxies are located at local matter density peaks formed through gravitational instabilities. Second, the observation of galaxies through their emitting light modifies our perception of their distribution. In particular, their measured redshift can be shifted with respect to their Friedmann-Lemaitre-Robertson-Walker counterpart due to peculiar velocities (RSD) and due to local or integrated gravitational potentials. This leads to distortions when estimating the distance to galaxies given a fiducial cosmology, and therefore modifies their overall perceived spatial distribution. Moreover, light rays follow null geodesics and can therefore be deflected by gravitational potentials along their trajectories, which impacts the observed angular position of sources. In the following, we present how the two dominant effects leading to the observed redshift-space galaxy two-point correlation function, that is, RSD and gravitational lensing magnification, can be modelled.

2.1. Redshift-space distortions

The standard way of modelling clustering in redshift space beyond linear theory (Peebles 1980; Kaiser 1987; Hamilton 1992) is to consider only redshift perturbations due to peculiar velocities. The real-space quasi-linear clustering can be predicted with different flavours of perturbation theory, either Eulerian (Bernardeau et al. 2002; Crocce & Scoccimarro 2006; Taruya et al. 2012, 2013) or Lagrangian (Zel’dovich 1970; Matsubara 2008a,b; Carlson et al. 2013). In addition, the mapping from real to redshift space can be done with different approaches. Common approaches in observational studies include the TNS model (Taruya et al. 2010, 2013) and the streaming model (Scoccimarro 2004; Reid & White 2011).

In this work, we consider the convolution Lagrangian perturbation theory (CLPT; Carlson et al. 2013) to predict real-space quantities. Generally, if the non-linear clustering in real space is well modelled with this approach, Lagrangian theories are not accurate enough on small scales in redshift space (White 2014). We therefore adopt the Gaussian streaming model (Reid & White 2011) for the mapping from real to redshift space. The joint use of CLPT and Gaussian streaming (CLPT-GS) has been shown to provide a good match to data and simulations (Wang et al. 2014) and is routinely used to model observations (e.g., Reid et al. 2012; Samushia et al. 2014; Zarrouk et al. 2018; Bautista et al. 2021).

2.1.1. Convolution Lagrangian perturbation theory

In the Lagrangian framework, the position x of an infinitesimal volume element is given by

$$\begin{array}{c}\hfill \mathit{x}(\mathit{q},t)=\mathit{q}+\mathbf{\Psi }(\mathit{q},t),\end{array}$$(1)

where q is the Lagrangian coordinate (initial position) and Ψ(q) the displacement field. The latter encodes the displacement of any mass element. For conciseness, in the following we omit its time dependence. In Lagrangian theories, the displacement field is assumed to be small, whereas for Eulerian theories, the density contrast is assumed to be small. The displacement field can thus be expanded as a perturbative series, $\mathbf{\Psi }=\sum _{n=1}^{\infty }{\mathbf{\Psi }}^{(n)}$, where the first-order solution is the well-known Zel’dovich approximation (Zel’dovich 1970).

Due to mass conservation, the relation between the density field of a volume element depends on its initial location through

$$\begin{array}{c}\hfill [1+\delta (\mathit{x})]{\mathrm{d}}^3\mathit{x}=[1+\delta (\mathit{q})]{\mathrm{d}}^3\mathit{q},\end{array}$$(2)

where δ is the matter density contrast. We have therefore that

$$\begin{array}{c}\hfill 1+\delta (\mathit{x})={\displaystyle \int }{\mathrm{d}}^3\mathit{q}{\delta }_{\mathrm{D}}(\mathit{x}-\mathit{q}-\mathbf{\Psi }),\end{array}$$(3)

where δ_D is the Dirac delta function. However, the matter density field is not observationally relevant for galaxy clustering, but we are rather interested in biased tracers of the underlying matter field. By assuming a local Lagrangian bias, we can write

$$\begin{array}{c}\hfill 1+{\delta }_g(\mathit{q})=F(\delta (\mathit{q})),\end{array}$$(4)

where δ_g is the density contrast of galaxies and F(δ) is the biasing function smoothed on a given scale. We note that at first order, the linear Eulerian and Lagrangian biases are simply related by $b_E^{(1)}=1+b_L^{(1)}$. In the following, we use the notation $b_n=b_L^{(n)}$, with b₁ = ⟨F′⟩ and b₂ = ⟨F″⟩ the first- and second-order Lagrangian bias parameters given by the expectation value of the first- and second-order derivatives of F(δ) respectively.

The counterpart of Eq. (3) for biased tracers is thus

$$\begin{array}{c}\hfill 1+{\delta }_g(\mathit{x})={\displaystyle \int }{\mathrm{d}}^3\mathit{q}F[\delta (\mathit{q})]{\delta }_{\mathrm{D}}(\mathit{x}-\mathit{q}-\mathbf{\Psi }).\end{array}$$(5)

From this expression, it is possible to derive the real-space galaxy two-point correlation function, mean pairwise velocity, and velocity dispersion as

$$\begin{array}{cc}\hfill 1+\xi (\mathit{r})& =\langle [1+{\delta }_g(\mathit{x})][1+{\delta }_g(\mathit{x}+\mathit{r})]\rangle \hfill \\ \hfill & ={\displaystyle \int }{\mathrm{d}}^3\mathit{q}{M}_0(\mathit{r},\mathit{q}),\hfill \end{array}$$(6)

$$\begin{array}{cc}\hfill v_{12}(\mathit{r})& =\frac{\langle [1+{\delta }_g(\mathit{x})][1+{\delta }_g(\mathit{x}+\mathit{r})][v_z(\mathit{x}+\mathit{r})-v_z(\mathit{x})]\rangle }{\langle [1+{\delta }_g(\mathit{x})][1+{\delta }_g(\mathit{x}+\mathit{r})]\rangle }\hfill \\ \hfill & ={[1+\xi (r)]}^{-1}{\displaystyle \int }{\mathrm{d}}^3\mathit{q}{M}_1(\mathit{r},\mathit{q}),\hfill \end{array}$$(7)

$$\begin{array}{cc}\hfill {\sigma }_{12}^2(\mathit{r})& =\frac{\langle [1+{\delta }_g(\mathit{x})][1+{\delta }_g(\mathit{x}+\mathit{r})]{[v_z(\mathit{x}+\mathit{r})-v_z(\mathit{x})]}^2\rangle }{\langle [1+{\delta }_g(\mathit{x})][1+{\delta }_g(\mathit{x}+\mathit{r})]\rangle }\hfill \\ \hfill & ={[1+\xi (r)]}^{-1}{\displaystyle \int }{\mathrm{d}}^3\mathit{q}{M}_2(\mathit{r},\mathit{q}),\hfill \end{array}$$(8)

where r = |r| is the scale, v_z is the velocity along the line of sight u_z with ${\mathit{v}}_z({\mathit{x}}_2)-v_z({\mathit{x}}_1)=({\dot{\mathit{x}}}_2-{\dot{\mathit{x}}}_1)·{\mathit{u}}_z$, and M₀, M₁ and M₂ are integration kernels containing the parameters of the model (Wang et al. 2014). The main difference between CLPT and other Lagrangian perturbation theories such as in Matsubara (2008b) lies in the resummation employed in the kernels.

2.1.2. Gaussian streaming model

The CLPT predictions for the real-space correlation function and the two first moments of the pairwise velocity distribution, that is, Eqs. (6)–(8), can be used in the Gaussian streaming model (Reid & White 2011) to predict the anisotropic two-point correlation function in redshift space. In this model, the redshift-space correlation function is recovered by convolving the real-space correlation function along the line of sight, with a scale-dependent pairwise velocity distribution taken to be Gaussian. This gives

$$\begin{array}{cc}& 1+{\xi }^s(r_{\perp },r_{\Vert })={\displaystyle \int }\frac{\mathrm{d}y}{\sqrt{2\pi {\sigma }^2(r,\upsilon )}}[1+\xi (r)]\hfill \\ \hfill & \times exp\{-\frac{{[r_{\Vert }-y-\upsilon v_{12}(r)]}^2}{2{\sigma }^2(r,\upsilon )}\},\hfill \end{array}$$(9)

where r_∥ and r_⊥ are the components of the separation vector parallel and perpendicular to the line of sight, $r^2=r_{\perp }^2+y^2$, υ = y/r is the end-point line-of-sight cosine angle definition, and ${\sigma }^2(r,\upsilon )={\sigma }_{12}^2(r,\upsilon )+{\sigma }_v^2$, with σ_v an additional velocity dispersion term that accounts for small-scale random motions in virialised objects. From Eq. (9), it is fairly easy to compute the multipole moments of the correlation function by integrating over the Legendre polynomials. We note that one limitation of the Gaussian streaming model as presented here is that it assumes the distant-observer approximation, meaning that it neglects wide-angle effects (Szalay et al. 1998; Szapudi 2004; Pápai & Szapudi 2008; Raccanelli et al. 2010; Reimberg et al. 2016; Castorina & White 2018; Beutler et al. 2019; Taruya et al. 2020). In our case, we restrict our analysis to relatively high redshifts (z > 1) so that we can safely work under this approximation. Our implementation of the CLPT-GS model is publicly available¹.

2.2. Magnification bias

Beyond the redshift-space distortions induced by peculiar velocities, the next dominant effect at z > 1 that modifies the apparent distribution of galaxies is magnification bias. It originates from the fact that some regions of the sky are magnified (demagnified) due to gravitational lensing, meaning that they contain fewer (more) sources than on average. If we consider a survey of galaxies selected in flux, gravitational lensing will modify the apparent fluxes of the galaxies, so that some of them will enter or exit the sample. This selection effect will impact the observed clustering of galaxies. It is worth noting that magnification bias is not the only effect that arises: additional effects ensue that depend on peculiar velocities or gravitational potential (Challinor & Lewis 2011). However, at high redshift, magnification bias is the dominant effect and is the focus of this work.

2.2.1. Impact on the number counts

The surface brightness of sources, defined as the flux per unit solid angle, is conserved through gravitational lensing. This implies that the apparent size and flux of a source are simultaneously modified as

$$\begin{array}{cc}& S^{\prime }(\mathit{\theta })=\mathrm{\mu }(\mathit{\theta })S(\mathit{\theta }),\hfill \end{array}$$(10)

$$\begin{array}{cc}\hfill & \mathrm{d}{\mathrm{\Omega }}^{\prime }(\mathit{\theta })=\mathrm{\mu }(\mathit{\theta })\mathrm{d}\mathrm{\Omega }(\mathit{\theta }),\hfill \end{array}$$(11)

where S, μ,  and dΩ are the flux, magnification, and solid angle, respectively, in the direction θ on the sky (in the following we omit the angular dependence), and a prime denotes a magnified quantity. The conservation of the number of sources can be written as

$$\begin{array}{c}\hfill n^{\prime }(m^{\prime })\mathrm{d}{m}^{\prime }\mathrm{d}{\mathrm{\Omega }}^{\prime }=n(m)\mathrm{d}m\mathrm{d}\mathrm{\Omega },\end{array}$$(12)

where n(m) is the number density of sources per unit of solid angle and per magnitude interval dm, m = −2.5log(S​)+C is the magnitude, and C is a constant. From this, it is straightforward to infer the magnified magnitude m′=m − 2.5log(μ). The total number of observed sources up to a given magnitude limit m_l is

$$\begin{array}{cc}& {\displaystyle {\int }_{-\infty }^{m_l}{n}^{\prime }(m^{\prime })\mathrm{d}{m}^{\prime }={\mathrm{\mu }}^{-1}{\int }_{-\infty }^{m_l}n(m^{\prime }+2.5log(\mathrm{\mu }))\mathrm{d}{m}^{\prime }}\hfill \end{array}$$(13)

$$\begin{array}{cc}\hfill & ={\mathrm{\mu }}^{-1}{\displaystyle {\int }_{-\infty }^{m_l+2.5log(\mathrm{\mu })}n(m)\mathrm{d}m.}\hfill \end{array}$$(14)

Ultimately, this can be rewritten in terms of cumulative number densities as

$$\begin{array}{c}\hfill n^{\prime }(<m_l)={\mathrm{\mu }}^{-1}n(<m_l+2.5log(\mathrm{\mu })).\end{array}$$(15)

One generally considers a simple model for the luminosity function, a power law such that (Schneider et al. 1992; Broadhurst et al. 1995)

$$\begin{array}{c}\hfill n(<m)\propto {10}^{\mathit{ms}},\end{array}$$(16)

where s is defined as the logarithmic slope of the cumulative distribution and is a property of the target sample. We note that for this specific function, cumulative and differential distributions have the same shape. By inserting Eq. (16) in Eq. (15), we obtain that

$$\begin{array}{cc}& {\mathrm{\Delta }}_{\mathrm{len}}={\mathrm{\mu }}^{2.5s-1}-1\hfill \end{array}$$(17)

$$\begin{array}{cc}\hfill & \approx (5s-2)\kappa ,\hfill \end{array}$$(18)

where Δ_len = n′( < m_l)/n(< m_l)−1 is the perturbation on the number count due to magnification bias in a given direction on the sky, and κ is the lensing convergence. In the second line we have performed a first-order Taylor expansion on the magnification as μ = 1 + 2κ (hereafter referred to as the ‘weak-lensing limit’), with |κ|≪1. Gravitational lensing induces two competing terms on the observed number counts as seen in Eq. (18): the s-dependent term implies that in magnified regions (μ > 1, κ > 0), the flux of sources increases so that we are more likely to find objects (and conversely, the chances of finding objects in demagnified regions are lower). The second term, which is usually referred to as ‘dilution bias’, describes the change in size of solid angles on the sky. Magnified regions occupy more space, so that for a constant density, we should find fewer objects in these regions than on average. These two effects cancel exactly for s = 0.4.

2.2.2. Two-point correlation function correction

To be consistent, the lensing correction associated with magnification bias on the correlation function should in principle be derived using the same theoretical framework as for RSD. However, these developments are beyond the scope of the present paper. We instead propose a simple correction based on linear theory, which can easily be used in addition to any RSD model.

We start from the observed galaxy number counts, which accounts for density, RSD, and lensing perturbations (the full expressions accounting for all the terms at first order in metric perturbations can be found in Yoo et al. 2009; Challinor & Lewis 2011; Bonvin & Durrer 2011),

$$\begin{array}{c}\hfill \mathrm{\Delta }={\mathrm{\Delta }}_{\mathrm{den}}+{\mathrm{\Delta }}_{\mathrm{rsd}}+{\mathrm{\Delta }}_{\mathrm{len}},\end{array}$$(19)

where Δ_den = bδ, b is the Eulerian linear bias, and Δ_rsd = −∂_rv_r/ℋ is the RSD component, where ∂_rv_r and ℋ are the gradient of the velocity field along the line of sight and the conformal Hubble parameter, respectively. The lensing perturbation Δ_len is that of Eq. (18). We note that the decomposition in Eq. (19) is only true at first order since it neglects higher-order lensing correlations.

Since the correlation function can be written as ξ(r) = ⟨Δ(x)Δ(x+r)⟩, the linear correction that comes from the addition of lensing magnification in the number counts is

$$\begin{array}{c}\hfill {\xi }_{\mathrm{corr}}(\mathit{r})={\xi }_{\mathrm{den}-\mathrm{len}}(\mathit{r})+{\xi }_{\mathrm{rsd}-\mathrm{len}}(\mathit{r})+{\xi }_{\mathrm{len}-\mathrm{len}}(\mathit{r}),\end{array}$$(20)

where ξ_A-B(r)≡⟨Δ_A(x)Δ_B(x+r)⟩. The expressions for the different terms in Eq. (20) are derived in Matsubara (2000), Hui et al. (2007, 2008), and in Tansella et al. (2018a,b) for the curved-sky case. Precisely, in the latter case we have

$$\begin{array}{c}\hfill {\xi }_{\mathrm{A}-\mathrm{B}}(\theta ,z_1,z_2)={\displaystyle \int }\frac{\mathrm{d}k}{k}{P}_{\mathrm{R}}(k){\mathcal{Q}}_k^{\text{A}-\text{B}}(\theta ,z_1,z_2),\end{array}$$(21)

where (θ, z₁, z₂) defines the separation vector in observed coordinates², P_R(k) is the primordial matter power spectrum, and the kernels ${\mathcal{Q}}_k^{\text{A}-\text{B}}$ with A − B = {den − len, rsd − len, len − len} read

$$\begin{array}{cc}& {\mathcal{Q}}_k^{\text{den-len}}(\theta ,z_1,z_2)=b(z_1)S_D(z_1)\left(\frac{2-5s}{2{\chi }_2}\right)\hfill \\ \hfill & {\displaystyle {\int }_0^{{\chi }_2}\mathrm{d}\lambda \frac{{\chi }_2-\lambda }{\lambda }{S}_{\varphi +\psi }(\lambda ){\zeta }^{0L}(k{\chi }_1,k\lambda ,\theta ),}\hfill \end{array}$$(22)

$$\begin{array}{cc}\hfill & {\mathcal{Q}}_k^{\text{rsd-len}}(\theta ,z_1,z_2)=\frac{k}{\mathcal{H}(z_1)}{S}_V(z_1)\left(\frac{2-5s}{2{\chi }_2}\right)\hfill \\ \hfill & {\displaystyle {\int }_0^{{\chi }_2}\mathrm{d}\lambda \frac{{\chi }_2-\lambda }{\lambda }{S}_{\varphi +\psi }(\lambda ){\zeta }^{2L}(k{\chi }_1,k\lambda ,\theta ),}\hfill \end{array}$$(23)

$$\begin{array}{cc}\hfill & {\mathcal{Q}}_k^{\text{len-len}}(\theta ,z_1,z_2)=\frac{{(2-5s)}^2}{4{\chi }_1{\chi }_2}\hfill \\ \hfill & {\displaystyle {\int }_0^{{\chi }_1}{\int }_0^{{\chi }_2}\mathrm{d}\lambda \mathrm{d}{\lambda }^{\prime }\frac{({\chi }_1-\lambda )({\chi }_2-{\lambda }^{\prime })}{\lambda {\lambda }^{\prime }}{S}_{\varphi +\psi }(\lambda )}\hfill \\ \hfill & S_{\varphi +\psi }({\lambda }^{\prime }){\zeta }^{\mathit{LL}}(k\lambda ,k{\lambda }^{\prime },\theta ).\hfill \end{array}$$(24)

In these equations, ζ are pure geometrical functions provided in Appendix B of Tansella et al. (2018a), χ_i is the comoving distance to redshift z_i, and S_D, S_V, and S_ϕ + ψ are the scaled transfer functions associated with density, peculiar velocity, and gravitational potentials, respectively. ξ_A-B(θ, z₁, z₂) in Eq. (21) can be written in terms of the separations parallel and perpendicular to the line of sight, r_∥ and r_⊥, using that r_∥ = χ₂ − χ₁, $r_{\perp }=\sqrt{r^2-r_{\Vert }^2}$, and $r=\sqrt{{\chi }_1^2+{\chi }_2^2-2{\chi }_1{\chi }_{2}cos\theta }$.

We implement the magnification bias correction using the COFFE library³ (Tansella et al. 2018b), which directly provides ξ_den-len, ξ_rsd-len, ξ_len-len in bins of (r_⊥, r_∥) using curved-sky linear theory and given an input linear power spectrum. Jelic-Cizmek (2021) noted that although there is no large differences between the curved-sky and flat-sky prescriptions at the scales of interest for us in general, that is r ≲ 150 h⁻¹ Mpc, the ξ_den-len component is quite sensitive to the adopted prescription. We therefore adopted the full curved-sky implementation.

Our theoretical model for the redshift-space correlation function therefore consists of the CLPT-GS prediction for non-linear RSD and of the curved-sky linear theory prediction for the additional lensing magnification correction. Formally, the anisotropic correlation function model is given by

$$\begin{array}{c}\hfill {\xi }_{\mathrm{model}}(r_{\perp },r_{\Vert })={\xi }_{\mathrm{CLPT}-\mathrm{GS}}(r_{\perp },r_{\Vert })+{\xi }_{\mathrm{corr}}(r_{\perp },r_{\Vert }).\end{array}$$(25)

A final step involves evaluating ξ_model at coordinates (r, ν) using that $r=\sqrt{r_{\perp }^2+r_{\Vert }^2}$ and ν = r_∥/r, and computing associated multipole moments as

$$\begin{array}{c}\hfill {\xi }_{\ell }^{\mathrm{model}}(r)=\frac{2\ell +1}{2}{\displaystyle {\int }_{-1}^1{\xi }_{\mathrm{model}}(r,\nu ){\mathcal{L}}_{\ell }(\nu )\mathrm{d}\nu ,}\end{array}$$(26)

where ℒ_ℓ is the Legendre polynomial of order ℓ.

We show the correlation function multipole moments computed with our model in Fig. 1. We considered here the matter at z = 1.8 in the ΛCDM model and s = 1.2. Both flat-sky and curved-sky linear lensing prescriptions are presented.

Fig. 1. Multipoles of the correlation function (monopole, quadrupole, and hexadecapole from top to bottom) when accounting for RSD only (red), for RSD and lensing magnification (blue), and for RSD and a flat-sky implementation of the lensing magnification. We consider here a ΛCDM model at z = 1.8 with a galaxy bias equal to unity and s = 1.2.

We first remark that magnification bias adds a positive contribution to the correlation function multipoles. The case s = 0 would also lead to positive contributions for the even multipoles because at this redshift, |ξ_{len − len}|> |ξ_{den − len}|, where |ξ_{len − len}| is strictly positive and the sign of |ξ_{den − len}| depends on the value of s. While this inequality might not hold at lower redshifts, it is valid for the redshift range we consider, that is, z > 1. Second, the full-sky and flat-sky implementations of the lensing correction give very similar results that are indistinguishable, except on the hexadecapole at large comoving separations. Overall, although we use the curved-sky correction in our modelling, the flat-sky approximation should also work in likelihood analyses, since the covariance associated with the hexadecapole should weakly affect the final results with respect to the dominant monopole and quadrupole.

3. Methods

To investigate the impact of magnification bias on clustering in redshift space, we used N-body simulations, which naturally account for the fully non-linear structure formation, and extracted light cones with various magnification bias selections. We then estimated the first three even multipole moments of the two-point correlation function in the light cones in several tomographic redshift bins and ran a Monte Carlo Markov chain (MCMC) likelihood analysis to sample the parameters of the model described in Sect. 2. We present the different methods in this section.

3.1. Datasets

The RayGal simulation suite⁴ (Breton et al. 2019; Rasera et al. 2022) is based on RAMSES (Teyssier 2002; Guillet & Teyssier 2011). These are dark-matter-only N-body simulations containing 4096³ dark matter (DM) particles of mass 1.8 × 10¹⁰ M_⊙ in a volume of 2.625³h⁻³ Gpc³. Both ΛCDM and wCDM versions are available, and associated fiducial cosmological parameters are given in Table 1. The two cosmologies have different Ω_m and σ₈ and therefore different values of fσ₈(z), since f ≈ Ω_m(z)^0.55 in General Relativity (Wang & Steinhardt 1998; Linder & Cahn 2007). This is shown in Fig. 2, where we show the fiducial values of fσ₈ as a function of redshift for the two cosmologies as well as the expectations from Planck Collaboration XIII (2016)ΛCDM best-fitting model assuming General Relativity. Interestingly, the values of fσ₈(z) for the RayGal ΛCDM (wCDM) model are close those from Planck at high (low) redshift. It is worth emphasising the importance of analysing simulations with different cosmologies, since we can analyse them blindly in a fiducial cosmology, as in observations, and see whether unbiased estimates of the growth rate of structure can be recovered.

Fig. 2. Evolution of the growth rate of structure as a function of redshift for the Planck Collaboration XIII (2016) cosmology in black (error bars are shown in grey), as well as the two cosmologies (ΛCDM in purple, wCDM in cyan) of the RayGal simulations until z = 2.

3.1.1. RayGal light cones

Several light cones have been extracted from the RayGal simulations. In the present work, we used light cones with an aperture of 2500 deg² extending to z = 2, which encompasses the redshift range probed by the DESI (DESI Collaboration 2016) and Euclid (Laureijs et al. 2011) surveys. These light cones contain DM particles as well as DM haloes identified with the parallel friend-of-friend algorithm pFoF (Roy et al. 2014), using a linking length of 0.2. We imposed a minimum of 100 particles per halo, which leads to haloes with masses higher than 1.8 × 10¹² M_⊙.

The gravitational-lensing information was computed in the light cones by using the ray-tracing code MAGRATHEA-PATHFINDER (Reverdy 2014; Breton & Reverdy 2022). This code implements an iterative algorithm that finds the null geodesics connecting the observer to each source (Breton et al. 2019), that is, either particles or haloes. This allows the computation of RSD and lensing effects at the same time in a general and accurate way. It is important to emphasize the fact that the treatment of gravitational lensing does not involve the Born approximation, which is often used. Our light cones contained about 1.2 × 10⁷ haloes in both cosmologies, and we ray-traced about 4 × 10⁸ randomly selected particles. Having both haloes and particles enabled us to perform a redshift-space clustering analysis on a biased population for the former (and hence, closer to observations), and for the latter, to carry out a precise study in which the number of matter tracers was maximised.

In Fig. 3 we show the redshift distribution of the halo and particle samples in the ΛCDM light cone, as well as the adopted tomographic redshift bins. The distributions in the wCDM light cone are very similar. The tomographic redshift bins cover a similar redshift range as present and future galaxy cosmological surveys, a regime in which gravitational lensing effects on galaxy clustering start to be significant (at about z > 1). The shape of the redshift distribution increases for particles monotonically, as expected in the case of constant density. However, at about z = 2, N(z) seems to decrease. This edge effect arises because we built our light cones up to z ∼ 2. To avoid any issue, we used a maximum redshift of z_max = 1.95. For haloes, N(z) reaches a maximum at around z = 1.2 and then decreases. This can be explained by the combined effect of the halo formation and limited mass resolution in the simulation. We did not impose any further selection in redshift to avoid discarding too many objects from our samples, and we thus maximised RSD and lensing magnification signals.

Fig. 3. Normalised redshift distribution of the DM particles and haloes, in blue and orange histograms, respectively, within the narrow cone of the ΛCDM RayGal simulation. The red, green, and yellow regions refer to the redshift bins, i.e. z = [1.0, 1.3],[1.3, 1.6], and [1.6, 1.95].

3.1.2. Implementation of the magnification bias

To reproduce different levels of magnification bias, a halo mass-galaxy luminosity relation might be imposed and apparent magnified fluxes or magnitudes might be estimated, from which then selections might be made. While it is clearly the appropriate methodology when constructing most realistic mock catalogues, our goal is to investigate the effect of the magnification bias on galaxy clustering in a general way, independent of the properties of any target sample. Hence, we found that the easiest and most efficient way to mimic the effect of magnification bias is to directly use the magnification of sources. Objects may be selected using a probability function proportional to μ^2.5s (see also Sect. 2.2.1). The advantages of this approach is that fewer sources are discarded, and the approach does not depend on the mass resolution of the simulation. However, it depends on some normalisation and discards more sources at high values of s. We therefore instead chose to weight each source by μ^2.5s. This allowed us to account for the effect of lensing magnification while keeping all the objects in our sample. We chose a weight equal to μ^2.5s (in addition to the observed angular positions in the galaxy clustering analysis) because our ray-tracing code computes the distortion matrix along the null geodesics that connect the observer to each source. In this case, the weak-lensing statistics in our sample is not the same as if we were using the Born approximation. Our averaging procedure indeed performs a ‘source averaging’ (Kibble & Lieu 2005; Bonvin et al. 2015; Kaiser & Peacock 2016; Breton & Fleury 2021). On average, light rays propagate in under-dense regions due to their path on the real null geodesics, which leads to a negative mean convergence. Had we used the Born approximation instead, our full sample would have been ‘unlensed’ and would have led to a vanishing mean convergence. In this case, to implement the effect of magnification in the mock catalogue, we would rather have used the comoving angular positions and applied a weight equal to μ^{2.5s − 1} to each source.

We note that in our theoretical modelling of the magnification bias, we assumed the weak-lensing limit, meaning that the convergence and magnification are small. We therefore used Eq. (18) instead of the exact Eq. (17) for the lensed number count. Figure 4 shows the impact of this approximation on the averaged convergence as a function of redshift. We first note that the mean convergence in the unlensed case (i.e. when s = 0.4) is very close to zero, demonstrating the validity of our method to implement the effect of lensing magnification. For s = 0.2, the difference between the exact and weak-lensing solution is very small. For higher values, the discrepancy clearly grows. The difference between the exact and weak-lensing solutions increases with s, so that for s = 1.2, it reaches ∼40% at z = 1.9. The question is whether this is a problem and how it affects the correlation function, since analytical prescriptions only work in the weak-lensing limit. The multipoles of the correlation function in Fig. A.1 are indeed different for s = 1.2, but not by the 40% of Fig. 4. This suggests that the first-order lensing correction should be used with caution in the weak-lensing limit, as its validity depends on redshift and on the value of s.

Fig. 4. Mean convergence as a function of redshift. The solid and dashed lines refer to the use of weights from the exact and approximate solutions, respectively. We calculated the convergence using the ΛCDM halo catalogue. The dark blue, blue, cyan, green, orange, and brown lines refer to the mean magnified convergence estimated from an unlensed sample with a factor s = 0.2, 0.4, 0.6, 0.8, 1.0, and 1.2, respectively.

3.2. Cosmological analysis

We now turn to the galaxy clustering analysis of the different samples described in Sect. 3.1. The total number of elements in all the studied cases are summarised in Table 2. Although RayGal simulations provide a full redshift decomposition at first order in metric perturbations, we only focus on the impact of lensing magnification beyond RSD and therefore ignore the more subtle effects that could affect the observed redshift. This means that we only perturb the redshift with the Doppler effect induced by peculiar velocities. The unlensed case corresponds to using the comoving angles, while the s = 0 case corresponds to using observed angles. For s > 0, we took observed angles and weight each source by μ^2.5s, as discussed in Sect. 3.1.2. In any configuration, the number of elements in Table 2 is the same (modulo tiny differences between observed and comoving angles due to the footprint). For particles, we selected a random sub-sample of the 4 × 10⁸ initial particles in our light cones, as the calculation of the anisotropic correlation is very expensive computationally.

3.2.1. Anisotropic two-point correlation function

The estimation of three-dimensional clustering necessitates the assumption of a fiducial cosmology to convert angular positions and redshifts into comoving separations. We assumed the ΛCDM cosmology of RayGal (Table 1) as our fiducial cosmology, which we used to analyse all the samples in Table 2, including those extracted from the wCDM light cone. The anisotropic correlation function was estimated with the Landy-Szalay estimator (Landy & Szalay 1993) as

$$\begin{array}{c}\hfill {\xi }_{\mathrm{LS}}(r,\nu )=\frac{DD(r,\nu )-2DR(r,\nu )+RR(r,\nu )}{RR(r,\nu )},\end{array}$$(27)

where DD, DR, and RR stand for data-data, data-random, and random-random pairs (weighted and normalised by the total number of elements), respectively. We used CORRFUNC (Sinha & Garrison 2020) to count the anisotropic number of pairs in bins of comoving separation r and cosine angle ν. The random samples contain 50 (20) times more objects than haloes (particles). We assigned redshifts in the random catalogues using the shuffling method, which consists of randomly picking redshifts from the data catalogue. Eventually, the multipole moments of the correlation function were obtained from

$$\begin{array}{c}\hfill {\xi }_{\ell }(r)=(2\ell +1){\displaystyle \sum _{\nu =0}^{\nu =1}\xi (r,\nu ){\mathcal{L}}_{\ell }(\nu )\mathrm{\Delta }\nu .}\end{array}$$(28)

We are only interested in the first non-vanishing even multipole moments of the correlation function: the monopole (ℓ = 0), the quadrupole (ℓ = 2), and the hexadecapole (ℓ = 4). We considered a scale range 27.5 to 127.5 h⁻¹Mpc with bins of 5 h⁻¹Mpc, and 200 bins in ν.

3.2.2. Covariance matrices

The covariance matrices on single multipole correlation function measurements are estimated analytically assuming Gaussianity as described in Grieb et al. (2016). Particularly, the covariance matrix between correlation function multipoles ℓ₁ and ℓ₂, and between scales r₁ and r₂, is

$$\begin{array}{c}\hfill C_{{\ell }_1,{\ell }_2}(r_1,r_2)=\frac{i^{{\ell }_1+{\ell }_2}}{2{\pi }^2}{\displaystyle {\int }_0^{\infty }{k}^2{\sigma }_{{\ell }_1,{\ell }_2}^2(k){\overline{j}}_{{\ell }_1}(k{r}_1){\overline{j}}_{{\ell }_2}(k{r}_2)\mathrm{d}k,}\end{array}$$(29)

where ${\overline{j}}_{\ell }$ are the bin-averaged spherical Bessel functions and ${\sigma }_{{\mathcal{l}}_1,{\mathcal{l}}_2}^2$ are the per-mode covariance multipole moments, both given in Grieb et al. (2016). The latter function is an integral over the anisotropic power spectrum, which was set here to the corresponding best-fitting RSD model to measurements.

3.2.3. Likelihood analysis

We performed a likelihood analysis of the measured monopole, quadrupole, and hexadecapole correlation functions in each sample. The likelihood ℒ is defined as

$$\begin{array}{c}\hfill -2ln\mathcal{L}(\mathit{\vartheta })={\displaystyle \sum _{i,j}^{N_p}{\mathbf{\Delta }}_i(\mathit{\vartheta })C_{\mathit{ij}}^{-1}{\mathbf{\Delta }}_j(\mathit{\vartheta }),}\end{array}$$(30)

where ϑ is the vector of parameters, Δ is the data-model difference vector, N_p is the total number of data points, and C is the covariance matrix. The model we used had six free parameters: $\vartheta =\{f,b_1,b_2,{\sigma }_v^2,{\alpha }_{\parallel },{\alpha }_{\perp }\}$, which correspond to the growth rate, first- and second-order Lagrangian bias parameters, squared velocity dispersion, and two dilation parameters that account for Alcock-Paczynski distortions, respectively. We note that we varied f and not directly fσ₈ because the model takes as input the linear power spectrum associated with the ΛCDM simulation at the redshift of interest. We therefore used a fiducial value of σ₈ to compute the theoretical prediction. We cannot let fσ₈ free mainly because the value of σ₈ is degenerate with the growth factor D₊(z). Within linear theory, this is not a problem as σ₈ can factor out. Within the framework of CLPT, we cannot do this because the non-linear part of the power spectrum is redshift dependent. In any case, although we let f free in the likelihood analysis, we eventually extracted fσ₈ and compared it to the fiducial value.

We used the dilation parameters along the parallel and transverse direction to account for any change in cosmology with respect to the fiducial one. Formally, these two parameters enter our formalism as a multiplicative factor on the scales in the anisotropic correlation function: ξ(α_⊥r_⊥, α_∥r_∥). If the fiducial cosmology is that of the data, we expect (α_⊥, α_∥) = (1, 1). Otherwise, these are given by the ratios

$$\begin{array}{cc}& {\alpha }_{\perp }=\frac{D_M(z)/r_d}{D_M^{\mathrm{fid}}(z)/r_d^{\mathrm{fid}}},\hfill \end{array}$$(31)

$$\begin{array}{cc}\hfill & {\alpha }_{\Vert }=\frac{D_H(z)/r_d}{D_H^{\mathrm{fid}}(z)/r_d^{\mathrm{fid}}},\hfill \end{array}$$(32)

where the superscript ‘fid’ refers to an estimation using the fiducial cosmology, D_M(z) = (1 + z)D_A(z) with D_A the angular diameter distance, D_H(z) = c/H(z), and r_d is the sound horizon at the drag epoch. We note that z refers here to the effective redshift of the sample. For the four redshift bins in our case, these are z = 0.857, 1.151, 1.448, and 1.769 for haloes and z = 0.86, 1.155, 1.453, and 1.777 for particles.

As argued in Sánchez (2020), using σ₈ might not be ideal for galaxy clustering analysis, since it relies on a scale given in h⁻¹ Mpc units that is cosmology dependent. However, as noted in Bautista et al. (2021), σ₈ can still be used as long as we take the changes due to dilation parameters into account. Instead of computing σ₈ at 8 h⁻¹ Mpc from a linear power spectrum at z = 0, this quantity should be evaluated at 8αh⁻¹ Mpc, where $\alpha ={\alpha }_{\perp }^{2/3}{\alpha }_{\parallel }^{1/3}$ is the isotropic dilation parameter. In the following, our estimates of fσ₈ implicitly account for this correction, while the fiducial values are those of Fig. 2.

We performed three types of analysis. First, a standard analysis that only accounted for RSD using the CLPT-GS model. This is similar to what is routinely done in observational studies of galaxy clustering. Second, an analysis that included the magnification bias correction in the model and where s was fixed to its fiducial value. Although it is possible to estimate s from the data itself, the standard method consists of fitting the local slope at the faint-end of the luminosity function, which might not be accurate when the selection function is complex (von Wietersheim-Kramsta et al. 2021). We therefore also propose an alternative analysis that includes the magnification bias correction in the model, but allows the s parameter to vary. Finally, we used the flat priors displayed in Table 3 (where only the last analysis allows s to vary). For each sample, we produced MCMCs using the EMCEE algorithm (Foreman-Mackey et al. 2013) with 100 walkers and 50 000 steps per walker. The walkers were initialised around the fiducial values (and a first-order Lagrangian bias of 0 and 1.4 for particles and haloes, respectively), to which we added a random shift using a Gaussian distribution with zero mean and variance unity, multiplied by 0.1. The burn-in phase depended on the auto-correlation time, and chains were thinned by removing highly correlated steps (the thinning was mostly used for convenience when making the plots, as we checked that it did not change the final results). Furthermore, the size of every chain was longer than 50 times the auto-correlation length. We also verified that we were able to recover the same results for the analysis of DM particles using MultiNest⁵ (Feroz & Hobson 2008; Feroz et al. 2009, 2019). The chains were finally analysed with GETDIST (Lewis 2019) to obtain the final parameter constraints.

4. Results

In this section, we present the results of the redshift-space clustering analysis with and without magnification bias for the particle and halo samples. We first discuss the recovered value fσ₈ and its bias with respect to the fiducial value. We then turn to an overview of the results for the other five (or six when s is free) parameters of the model.

4.1. Bias on fσ₈

4.1.1. Particles

We first consider the case in which sources are DM particles, as these should give the clearest trend on the impact of magnification bias. The full results for the ΛCDM case are shown in Fig. 5.

Fig. 5. Relative difference on fσ8 with respect to its fiducial value as a function of redshift bins and sample configuration for ΛCDM particles. The data points in blue, green, and red refer to different analyses where we only account for RSD, add magnification bias with s fixed and s free, respectively (see also Sect. 3.2.3). The central values are those of the overall best fit, while error bars are those of the 68% confidence interval on the marginalised distribution. The dashed blue lines show the value for the unmagnified case. The grey shades highlight the ±5% and 10% ranges.

In the standard analysis (RSD only), the top row for the unlensed sample shows that we recover the expected value of fσ₈ within or close to 1σ statistical uncertainty. When lensing magnification is incorporated in the data, it clearly has almost no impact in the first redshift bin at z = 1.0 − 1.3, but strongly shifts the estimate of fσ₈ at higher redshifts. More precisely, magnification bias leads to an underestimation of the growth rate if it is not modelled. The shift on fσ₈ does not strictly increase with s, which makes sense as the most dominant term in the lensing correction is ξ_lens-lens, which scales as (5s − 2)² (see Sect. 2.2). This means that s = 0 should lead to a similar effect as s = 0.8, while s = 0.4 completely cancels the lensing effect. Then, the shift on the growth rate monotonically increases for s > 0.8. At z = 1.3 − 1.6 and z = 1.6 − 1.95, the shift on fσ₈ for s = 1.2 reaches approximately 9% and 11%, respectively, with respect to the unlensed case, while it is 5% and 6% for s = 1. In the highest redshift bin, we find that the growth rate is not recovered within 1σ statistical uncertainty for s ≳ 1.

When magnification bias is accounted for in the modelling and s is known, we obtain an unbiased estimate of the growth rate with respect to the unmagnified case, except for s = 1.2 in the highest redshift bin. This shows that our linear-theory-based lensing correction is good enough overall to correct for the effect of magnification bias. However, for very high values of s at high redshift, the weak-lensing limit that was used to compute the correction might no longer hold (see also the difference in the quadrupole in Fig. A.1, where the theoretical prediction seems to underestimate the signal). Therefore, accounting for higher-order terms may be important in the future for the most extreme configurations. In Fig. 6 we show the relative difference on fσ₈ with respect to the unmagnified case when the magnification implementation in the data is done in the weak-lensing limit. As described in Sects. 2.2.1 and 3.1.2, to account for this approximation, we need to use a weight equal to μ(1 + (5s − 2)κ) instead of μ^2.5s, the prefactor μ is used to correct for the dilation effect already present in the data. Within the weak-lensing limit (in purple), the relative difference increases only weakly with s and reaches a lower value at s = 1.2 (compared with the other cases), and it remains within 1σ with respect to the expected value. The difference between the green and purple data points only comes from the weak-lensing limit, and it reaches roughly 8% at s = 1.2. This shows that it is important to consider higher-order modelling for high-redshift samples with high s values, beyond the first-order weak-lensing limit.

Fig. 6. Relative difference on fσ8 with respect to the ‘no magnification’ case of Fig. 5. The results for the standard analysis and fixed s of Fig. 5 are shown in blue and green, and the results of a fixed s analysis are shown in purple. For this last analsyis, we implemented magnification in data under the weak-lensing limit, i.e. we used a weight equal to μ(1 + (5s − 2)κ) instead of μ2.5s (see also Sect. 3.1.2).

Lastly, when s is free in the fit, the error bars on the parameters are larger, which is expected since there is more freedom in the model. Letting s free might not be the best strategy in the full-shape RSD analysis. Estimated values of s appear to be better suited. Nonetheless, the best-fit values obtained when letting s free are very close to the case where s is fixed, meaning that we are able to estimate this nuisance parameter without any strong bias.

Overall, we find that depending on the sample selection, galaxy clustering surveys probing redshifts above z ∼ 1.3 will need to properly include lensing magnification as part of the theoretical model to recover unbiased estimates of the growth rate of structure. This will in turn enable them to test gravity while keeping theoretical systematic errors under control.

In Fig. 7 we show the same analysis, but for the wCDM model. The results are extremely similar to those of Fig. 5, which is reassuring since here we compare the estimated fσ₈ with the fiducial values in the wCDM cosmology. This means that, although the full analysis is performed using a ΛCDM fiducial model, we are still able to recover almost unbiased constraints on the cosmology of the data (see also Fig. 2 for the difference in fσ₈ between the two cosmological models). The error bars for s = 0.2 and s = 0.6 in the highest redshift bin are very large. The fact that these configurations give similar results is expected since the leading magnification bias term scales as (5s − 2)². Furthermore, it seems that for these particular configurations, there is a multi-modal behaviour in the joint posterior probability for the growth rate and second-order Lagrangian bias parameters, leading to larger error bars. In the last two redshift bins, the shift on fσ₈ reaches approximately 12% and 10% for s = 1.2 when the magnification bias is not accounted for. As previously, knowing the true value of s in the modelling gives an unbiased estimate of fσ₈, except for the last redshift bin for s = 0.8, which seems slightly overestimated. This shows that although the lensing correction computed with the ΛCDM model is not exact for the wCDM simulation as seen in Appendix A, it still includes most of the effect. However, there can be large discrepancies at high redshift for high values of s for the lensing contribution with respect to the ΛCDM model (see Fig. A.1), and it is therefore not clear that using a correction from the fiducial model is enough when the cosmologies are very different. Even in this case, however, it seems that in our data, adding the magnification bias correction at fixed s gives a much better agreement to the fiducial value of fσ₈. In the future, it will thus be important to properly model the magnification bias term, as it could lead to additional bias.

Fig. 7. Same as Fig. 5, but for wCDM.

4.1.2. Haloes

We now turn to the analysis of DM haloes. While DM particles are useful since they allow characterising the effect of magnification bias with great precision, analysing haloes is more observationally relevant since, just like galaxies, they are biased tracers of the matter density field. This is particularly important because the relative effect of magnification bias on the total number counts scales as (5s − 2)/b, where b is the linear Eulerian galaxy bias. This means that we expect lensing to have less impact on samples with large bias. In our case, the best-fit RSD-only models for the unmagnified halo samples in the ΛCDM cosmology gives roughly b = 2.1, 2.5, and 3 in the different redshift bins, from low to high redshifts. We show the results in Fig. 8. Qualitatively, the results are similar to those of Fig. 5. We find that the values of the growth rate tend to be slightly overestimated, even when the sample is not magnified. This can be attributed to residual theoretical uncertainties of the RSD model and to the modest volume probed by our light cones, which can lead to non-negligible sample variance effects. Nonetheless, we can study and discuss the relative differences in the estimated parameters for the various cases. As previously, we find that the central values of fσ₈ are shifted towards lower values when lensing is not accounted for in the modelling. The main difference between DM particles and haloes is that for the latter, the effect of lensing is not clear as the magnified and unmagnified cases agree within 1σ error. In the highest redshift bin, the best-fit value for s = 1.2 is shifted by roughly 3% with respect to the unmagnified case. While it is still large, it does not reach the 11% seen for the particles.

Fig. 8. Same as Fig. 5, but for haloes.

The results displayed in Fig. 9 for wCDM haloes are very similar. We find a discrepancy of 3.5% and 5.5% in the highest redshift bin for s = 1.0, 1.2 with respect to the unmagnified case. While larger than for ΛCDM haloes, it still does not reach the 10% difference that we have for particles.

Fig. 9. Same as Fig. 8, but for wCDM.

Overall, we find that the effect of magnification bias is lower for haloes than for particles, which is expected due to galaxy bias. Moreover, the relative difference on fσ₈ is suppressed by a factor close to b. This means that very biased target samples might not be significantly subject to lensing effects at the considered redshifts in the present work.

4.2. Other parameters

In this section we study the impact of magnification bias on the parameters of our model beyond the growth rate of structure. We consider ΛCDM particles in two cases: s = 1.0, as it is a realistic case for future spectroscopic surveys (Lepori et al. 2022), and s = 1.2 where the effect of magnification bias is the most significant. We present the posterior distributions for s = 1.0 in Fig. 10. When lensing is not accounted for, the estimate of the growth rate is biased, while the lensing correction (when s is fixed or free) allows us to recover the fiducial value within an error of 1σ. These results are approximately the same as those of Fig. 5, except that we show f (the free parameter of our model) and not fσ₈. Since our fiducial cosmology is already that of the data, there should not be any qualitative difference, in the sense that we only need to multiply f by the fiducial σ₈ (because the dilation parameters should be very close to unity in principle). For DM particles, the first-order Lagrangian bias b₁ should be equal to zero, and this is indeed what we obtain. However, not accounting for lensing magnification leads to an overestimation of the bias. This comes from the fact that lensing adds a positive contribution to the multipole moments (as shown in Fig. 1), which is counter-balanced in the likelihood analyses by higher values of the galaxy bias. Moreover, the (f, b₁) joint posterior probability clearly shows an anti-correlation, where higher values of the linear bias induce lower values of f. This comes from the fact that higher values of the galaxy bias lead to a larger amplitude of the quadrupole. In order to fit the data, the likelihood analysis therefore converges towards a lower value of f.

Fig. 10. Posterior distributions for the parameters of our model for the ΛCDM particle sample with s = 1.0 at z = 1.6 − 1.95 (see also the full results in Fig. 5).

For the growth rate and galaxy bias parameters, the error bars increase when s is let free with respect to the case when s is fixed. However, this is not the case for the other nuisance parameters, that is, b₂, ${\sigma }_v^2$, and dilation parameters, where the two analyses give similar marginal distributions. Regarding the estimation of s, we see a remarkable agreement between the expected and best-fit values. This shows that, in principle, we should be able to accurately recover s through a galaxy clustering analysis. Furthermore, we notice that the marginal distribution of s exhibits a multi-modal behaviour around s = 0. This comes from the fact that the leading correction term is proportional to (5s − 2)², and there are two plausible possibilities when s ≲ 0.8. This is a further evidence that it might not be optimal to let s free with a flat uninformative prior.

Finally, we find that in none of the cases, (α_⊥, α_∥) = (1, 1) is recovered exactly within an error of 1σ error, although our results are very close. This might be related to the fact that we have one particular realisation of the density field in a limited volume, which is subject to sample variance. We note that in the RSD-only case, the estimated α_⊥ and α_∥ parameters depart from fiducial values more than when lensing is accounted for. This means that in this case, the likelihood analysis prefers other cosmological models. This might further impact the estimation of the dilation-corrected σ₈, as described in Sect. 3.2.3. However, this additional bias seems low for the redshift bins that we considered. In the case of the second-order Lagrangian bias b₂, the fiducial value is only recovered when lensing is accounted for, while the inferred values of the velocity dispersion parameter ${\sigma }_v^2$ are very similar in all cases. Nonetheless, these two parameters are considered as nuisance parameters over which we marginalised.

Lastly, we present the posterior distributions for s = 1.2 in Fig. 11. The results are similar to those of Fig. 10, except that here we do not recover the value of f within an error of 1σ, even when we account for magnification bias as seen in Fig. 5. As in Fig. 10, we remark that we are able to recover s accurately, even in this configuration, in which the weak-lensing limit might start to hold no longer. We also note that there is no multi-modal behaviour, in contrast to the case with s = 1.0. This is because s is large enough so that there is no ambiguity between the different solutions at s < 0.8. This shows that, generally, adding the lensing contribution in the modelling allows a considerably better recovery of the fiducial parameters.

Fig. 11. Same as Fig. 10, but with s = 1.2.

5. Conclusion

We have studied the impact of lensing magnification on the three-dimensional galaxy clustering in redshift space and how the growth rate of structure is affected when lensing magnification is not properly accounted for in the modelling. This effect, commonly referred to as magnification bias, comes from the fact that gravitational lensing modifies the apparent angular positions of sources and magnifies their fluxes. Therefore, magnitude-limited surveys are sensitive to it as it generates additional apparent spatial correlations in the data. We performed an exhaustive galaxy clustering analysis of the multipole moments of the correlation function, where the theoretical model relied on the convolution Lagrangian perturbation theory and Gaussian streaming model for non-linear RSD and on a curved-sky, linear-theory-based lensing correction. We compared this model to high-resolution N-body simulations with two different cosmologies, ΛCDM and wCDM, where the effects of RSD and magnification bias were implemented in simulated catalogues of DM particles and haloes.

The main results of our paper are the following. First, magnification bias only becomes relevant at z > 1.3. Second, not accounting for magnification bias in the modelling gives a biased estimate of the growth rate. We found that analysing a magnified sample with a theoretical model that only contains RSD leads to an overestimation of the galaxy bias, and more importantly, to an underestimation of the growth rate. Depending on the redshift and slope of the galaxy luminosity function, the best-fit value of fσ₈ can be shifted by more than 10%. It appears that using the linear-theory lensing correction allows the recovery of unbiased estimates of the growth rate in most cases when the lensing contribution is small if s is accurately known a priori. We also showed that for the cases with high redshift, large s, and low bias, it could be important to model the lensing contribution beyond the first-order weak-lensing limit. If s is unknown, it might be tempting to keep it as a free parameter of the model. While this allows the recovery of an unbiased estimate of s, it significantly increases the uncertainty on fσ₈, which we might refrain from. This shows that s should not be let free (or at least, no flat uninformative prior should be used on that parameter) and a way should be found instead to estimate it differently. Finally, we find that we were able to recover the same results for a wCDM simulation, although it was analysed with a fiducial ΛCDM model. This is encouraging and shows that in our case, the modelling is reliable even if we do not know the underlying cosmology of the data exactly. However, further investigations might be needed to estimate the impact of the cosmological model on the lensing correction as it could in principle be another source of bias.

Gravitational lensing will play an important role in galaxy clustering analyses in future high-redshift surveys, and its effect will have to be implemented in order to accurately recover cosmological parameters. This effect has not been considered in observational studies so far, even when the redshift of the sample was high (Okumura et al. 2016; Zarrouk et al. 2018; Hou et al. 2021). However, this does not necessarily mean that these works were not correct because statistical uncertainties are still significant and the effect of the magnification bias depends on the properties of the galaxy sample under scrutiny. An interesting prospect for future works would be to consistently incorporate the effect of gravitational lensing within the theoretical framework of non-linear RSD. Nonetheless, we showed that in our case, a hybrid model of non-linear RSD and linear lensing is in most cases sufficient for the accuracy of our simulated data.

Acknowledgments

This work was granted access to HPC resources of TGCC/CINES through allocations made by GENCI (Grand Équipement National de Calcul Intensif) under the allocation 2020-A0070402287. The project leading to this publication has received funding from Excellence Initiative of Aix-Marseille University – A*MIDEX, a French “Investissements d’Avenir” programme (AMX-19-IET-008 – IPhU).

References

Appendix A: Magnification bias modelling of the correlation function multipoles

In this appendix we provide more details about the effect of the magnification bias on the multipole moments of the correlation function. We consider the DM particles in the highest studied redshift bin at about z = 1.8, which maximises the magnification bias signal. The magnification bias corrections on the monopole, quadrupole, and hexadecapole of the correlation function are shown in Fig. A.1.

Fig. A.1. Absolute difference in the correlation function multipoles for DM particles at the highest redshift bin (i.e. z = 1.6 − 1.95) when accounting for magnification bias with respect to the case with RSD only (left). The red and green points refer to the monopole and quadrupole in ΛCDM (left) and wCDM (right), while the red, green, and black lines show the ΛCDM theoretical prediction for the monopole, quadrupole, and hexadecapole computed with COFFE

First, we note that in any case, that is, for s = 0 or 1.2, ΛCDM or wCDM, the lensing contribution to the correlation function is always positive at this redshift. This is interesting because it confirms the results of Fig. 1 and explains why, if lensing magnification is not incorporated in the model, the likelihood analysis tries to compensate for this lack by increasing the value of the bias parameter (mostly due to the monopole) and therefore lower the value of the growth rate (as described in Section 4) due to the quadrupole.

Secondly, for the case with s = 0, that is when we use observed angles instead of comoving ones, the monopole seems to agree with the theoretical prediction up to 60 h⁻¹Mpc only. This is surprising as we would expect a good agreement at large scales. This difference might be due to the large variance inherent to these scales or to an inaccuracy in the modelling as well as the non-linearity of lensing corrections (see also Hui et al. 2007 for a discussion). We also note the remarkable agreement between data and prediction for the quadrupole (while the hexadecapole seems underestimated in the prediction).

Thirdly, focusing on the case with s = 1.2, which give the most significant trend, the theoretical prediction overestimates the monopole and underestimates the quadrupole and hexadecapole above 80 h⁻¹Mpc. While this discrepancy is not necessarily large, especially for this high value of s, it impacts the estimation of the growth rate (see Fig. 5). This suggests that we should be careful about the first-order solution given within the weak-lensing limit, as higher-order terms might need to be accounted for at higher redshifts for high values of s.

Lastly, the theoretical prediction in the fiducial ΛCDM model does not agree very well with the data from the wCDM simulation in any case. Although the shape is similar between data points and analytical prediction, the amplitude is different (this difference is clear for the quadrupole, where data points are consistently at least a factor two higher than the prediction for both s = 1.2). For more precise studies at higher redshift, it will be important to find a way to account for this cosmology-dependant correction. Nonetheless, even if the present modelling in the fiducial ΛCDM cosmology is not perfect for analysing the wCDM, it is still much better than not accounting for lensing magnification at all.

All Tables

All Figures

	Fig. 1. Multipoles of the correlation function (monopole, quadrupole, and hexadecapole from top to bottom) when accounting for RSD only (red), for RSD and lensing magnification (blue), and for RSD and a flat-sky implementation of the lensing magnification. We consider here a ΛCDM model at z = 1.8 with a galaxy bias equal to unity and s = 1.2.
In the text

	Fig. 2. Evolution of the growth rate of structure as a function of redshift for the Planck Collaboration XIII (2016) cosmology in black (error bars are shown in grey), as well as the two cosmologies (ΛCDM in purple, wCDM in cyan) of the RayGal simulations until z = 2.
In the text

	Fig. 3. Normalised redshift distribution of the DM particles and haloes, in blue and orange histograms, respectively, within the narrow cone of the ΛCDM RayGal simulation. The red, green, and yellow regions refer to the redshift bins, i.e. z = [1.0, 1.3],[1.3, 1.6], and [1.6, 1.95].
In the text

Fig. 4. Mean convergence as a function of redshift. The solid and dashed lines refer to the use of weights from the exact and approximate solutions, respectively. We calculated the convergence using the ΛCDM halo catalogue. The dark blue, blue, cyan, green, orange, and brown lines refer to the mean magnified convergence estimated from an unlensed sample with a factor s = 0.2, 0.4, 0.6, 0.8, 1.0, and 1.2, respectively.

In the text

Fig. 5. Relative difference on fσ8 with respect to its fiducial value as a function of redshift bins and sample configuration for ΛCDM particles. The data points in blue, green, and red refer to different analyses where we only account for RSD, add magnification bias with s fixed and s free, respectively (see also Sect. 3.2.3). The central values are those of the overall best fit, while error bars are those of the 68% confidence interval on the marginalised distribution. The dashed blue lines show the value for the unmagnified case. The grey shades highlight the ±5% and 10% ranges.

In the text

Fig. 6. Relative difference on fσ8 with respect to the ‘no magnification’ case of Fig. 5. The results for the standard analysis and fixed s of Fig. 5 are shown in blue and green, and the results of a fixed s analysis are shown in purple. For this last analsyis, we implemented magnification in data under the weak-lensing limit, i.e. we used a weight equal to μ(1 + (5s − 2)κ) instead of μ2.5s (see also Sect. 3.1.2).

In the text

	Fig. 7. Same as Fig. 5, but for wCDM.
In the text

	Fig. 8. Same as Fig. 5, but for haloes.
In the text

	Fig. 9. Same as Fig. 8, but for wCDM.
In the text

	Fig. 10. Posterior distributions for the parameters of our model for the ΛCDM particle sample with s = 1.0 at z = 1.6 − 1.95 (see also the full results in Fig. 5).
In the text

	Fig. 11. Same as Fig. 10, but with s = 1.2.
In the text

Fig. A.1. Absolute difference in the correlation function multipoles for DM particles at the highest redshift bin (i.e. z = 1.6 − 1.95) when accounting for magnification bias with respect to the case with RSD only (left). The red and green points refer to the monopole and quadrupole in ΛCDM (left) and wCDM (right), while the red, green, and black lines show the ΛCDM theoretical prediction for the monopole, quadrupole, and hexadecapole computed with COFFE

In the text