© The Authors 2024

1. Introduction

The European Space Agency’s Euclid satellite mission (Laureijs et al. 2011; Amendola et al. 2018) aims at shedding light on the so-called dark components of the Universe, namely dark matter and dark energy. Dark matter, a mysterious form of matter that does not seem to emit light, yet accounts for more than 80% of the total matter content of the Universe, forms the bulk of the large-scale cosmic structure, upon which galaxies form and evolve. Dark energy is even more elusive, and it is what drives the current accelerated expansion of the Universe, contributing to about 70% of the total cosmic energy budget (see e.g. Bull et al. 2016, for a review of the current concordance cosmological model and the main theoretical challenges it faces). In fact, there is another possible explanation for the effects we ascribe to dark matter and/or dark energy: that the theory we use to analyse the data is incorrect. This approach is called “modified gravity” (e.g. Clifton et al. 2012). Thanks to the extent and exquisite precision of Euclid’s data, we shall soon be able to further test general relativity on scales far from the strong-gravity regime, where it has been tested to supreme precision (see e.g. Cardoso & Pani 2019, for a review of the current status).

Euclid will consist of two primary probes: a catalogue of about 30 million galaxies with spectroscopic redshift information, spanning a redshift range between z = 0.8 and z = 1.8, and a catalogue of 1.5 billion galaxy images with photometric redshifts down to z = 2 (see Laureijs et al. 2011; Amendola et al. 2018, for further details on the specifics of the Euclid surveys). One of the main goals of the Euclid spectroscopic survey is to measure the so-called growth rate, which is very sensitive to the theory of gravity (see for example Alam et al. 2017). However, in order to robustly test alternatives to general relativity, it is crucial to take all of the relevant effects into account in the analysis. One effect that has been overlooked in previous forecasts regarding the performance of the Euclid spectroscopic survey is lensing magnification (Matsubara 2004).

The aim of this paper is to investigate whether lensing magnification has to be included in this analysis. It is well known that lensing magnification has to be included in a photometric survey for the correct estimation of cosmological parameters (Duncan et al. 2014; Cardona et al. 2016; Villa et al. 2018; Lorenz et al. 2018; Unruh et al. 2020; Euclid Collaboration 2022; Mahony et al. 2022; Elvin-Poole et al. 2023). However, as the density and redshift space distortion (RSD) contributions are significantly larger in a spectroscopic survey, one might hope that lensing magnification can be neglected in this case. In this paper we show that this is not the case, and that neglecting lensing can shift the inferred cosmological parameters by up to 0.7σ and can affect the measured growth rate by up to 1σ. We then propose a method for reducing the shifts to an acceptable level. This method consists of adding the lensing magnification signal to the modelling, using fixed cosmological parameters within the Λ cold dark matter (CDM). Of course, in this way the lensing magnification is not exactly correct (since we do not know the theory of gravity nor the cosmological parameters), but we show that this is enough to de-bias the analysis, reducing the shifts to less than 0.1σ.

The paper is structured as follows. In the next section, we present the fluctuations of galaxy number counts within linear perturbation theory, concentrating on the redshift-space two-point correlation function (2PCF). In Sect. 3 we present the relevant quantities from the Euclid Flagship simulations used in this work. In Sect. 4 we explain the methods used in our analysis, Fisher matrix and Markov chain Monte Carlo (MCMC). In Sect. 5 we discuss our results and present the method for de-biasing the analysis, and in Sect. 6 we conclude. Some details and complementary results are presented in the appendices.

In this paper, scalar metric perturbations are described via the gauge-invariant dimensionless Bardeen potentials, Φ and Ψ. In longitudinal gauge the perturbed metric is

$$\begin{array}{c}\hfill \mathrm{d}{s}^2=a^2(\eta )[-(1+2\Psi )c^2\mathrm{d}{\eta }^2+(1-2\Phi ){\delta }_{\mathit{ij}}\mathrm{d}{x}^i\mathrm{d}{x}^j],\end{array}$$(1)

where we use the Einstein summation convention over repeated indices. Here a(η) is the cosmic scale factor evaluated at conformal time η, and c is the speed of light. In the above, as well as in subsequent equations, a prime denotes the derivative with respect to conformal time, and ℋ = a′/a = Ha denotes the conformal Hubble parameter. We normalised the scale factor to 1 (i.e. a₀ = 1) such that ℋ₀ = H₀.

2. Fluctuations of spectroscopic galaxy number counts

2.1. Galaxy number counts

An important observable of the Euclid satellite will be the galaxy number counts, that is, the number of galaxies dN(n, z) detected in a given small redshift bin dz around a redshift z and a small solid angle d∠ around a direction n. Expressing dN(n, z) = n(n, z) dz d∠ in terms of the angular-redshift galaxy density n(n, z) and subtracting the mean,

$$\begin{array}{c}\hfill \overline{n}(z)=\frac{1}{4\pi }{\displaystyle {\int }_{\mathrm{\Omega }}n(\mathit{n},z)\mathrm{d}\mathrm{\Omega },}\end{array}$$(2)

we define the galaxy number count fluctuation as

$$\begin{array}{c}\hfill \mathrm{\Delta }(\mathit{n},z)=\frac{n(\mathit{n},z)-\overline{n}(z)}{\overline{n}(z)}·\end{array}$$(3)

This quantity and its power spectra have been calculated at first order in cosmological perturbation theory in Yoo et al. (2009), Yoo (2010), Bonvin & Durrer (2011), Challinor & Lewis (2011), and Jeong et al. (2012). As it is an observable, the result is gauge invariant. It is not simply given by the density fluctuation on the constant redshift hypersurface, but also contains volume distortions. Most notable of those are the radial volume distortion from peculiar velocities, the so-called RSDs (Kaiser 1987), but also the transversal volume distortion due to weak lensing magnification (Matsubara 2004) and the large-scale relativistic effects identified for the first time in the above references. The final formula, including the effect of evolution bias (Challinor & Lewis 2011; Jeong et al. 2012), is given by (Di Dio et al. 2013)

$$\begin{array}{cc}\hfill \mathrm{\Delta }(\mathit{n},z,m_{\ast })& =b\delta +\frac{1}{\mathcal{H}}{\partial }_r(\mathit{V}·\mathit{n}){|}_{r_{\mathrm{s}}}\hfill \\ \hfill & -\frac{2-5s}{2r_{\mathrm{s}}}{\displaystyle {\int }_0^{r_{\mathrm{s}}}\mathrm{d}r\frac{r_{\mathrm{s}}-r}{r}{\mathrm{\Delta }}_{\mathrm{\Omega }}(\Phi +\Psi )+\dots }\hfill \end{array}$$(4)

Here r = r(z) is the comoving distance evaluated at redshift z, r_s is the comoving distance between the observer and the source, b = b(z, m_*) is the linear bias of galaxies with magnitude below m_*, the magnitude limit of the survey, δ is the gauge-invariant density fluctuation representing the density in comoving gauge, ∂_r is the derivative w.r.t. the comoving distance r, V is the velocity of sources in the longitudinal gauge, and Δ_∠ denotes the angular part of the Laplacian.

The function s = s(z, m_*) is the local count slope needed to determine the magnification bias. The local count slope depends on the magnitude limit m_* of the survey and is given by (see e.g. Challinor & Lewis 2011)

$$\begin{array}{c}\hfill s(z,m_{\ast })=\frac{\partial {log}_{10}\overline{N}(z,m<m_{\ast })}{\partial m_{\ast }}·\end{array}$$(5)

Here $\overline{N}$ is the cumulative number of objects brighter than the magnitude cut m_* (for a magnitude-limited sample). The first line of Eq. (4) corresponds to the standard terms of density fluctuations and RSDs, and the second line is the lensing magnification, which is the subject of the present paper (see Bonvin & Durrer 2011; Challinor & Lewis 2011; Jeong et al. 2012; Di Dio et al. 2013 for details)¹. The modelling of magnification in Eqs. (4) and (5) assumes an idealised magnitude-limited sample. On the other hand, realistically the sample selection may also depend on galaxy size, which is also impacted by magnification, or it may be based on a complex colour-magnitude selection. While these additional effects complicate the estimation of an effective local count slope from real data (see e.g. von Wietersheim-Kramsta et al. 2021), our treatment is adequate enough for our forecast. The dots at the end of Eq. (4) stand for the “large-scale relativistic terms” that we omitted in our analysis. These terms are suppressed by factors λℋ c⁻¹, where λ is the comoving wavelength of the perturbations (see e.g. Di Dio et al. 2013; Jelic-Cizmek et al. 2021; Euclid Collaboration 2022). It is well known that the large-scale relativistic terms are relevant only at very large scales and do not significantly impact the even multipoles of the correlation function on sub-Hubble scales (Lorenz et al. 2018; Yoo et al. 2009; Bonvin & Durrer 2011). Some of these relativistic terms will, however, be detectable by measuring odd multipoles in the correlation of two different tracers (see Bonvin et al. 2014, 2023; Gaztanaga et al. 2017; Lepori et al. 2020; Beutler & Di Dio 2020; Saga et al. 2022). A detailed study of the signal-to-noise ratio of all relativistic effects in simulated mock catalogues adapted to Euclid’s spectroscopic survey will be presented in Euclid Collaboration: Elkhashab et al. (in prep.).

The goal of this paper is to study the impact of lensing magnification on the 2PCF. We note that the impact of lensing magnification on the angular power spectrum, C_ℓ(z₁, z₂), has already been computed and found to be relevant for Euclid’s photometric sample (Euclid Collaboration 2022). However, the angular power spectrum, C_ℓ(z₁, z₂) is not well suited to a survey with spectroscopic resolution of δz ≲ 10⁻³ since we would have to split the redshift interval into more than 1000 bins in order to fully profit from the redshift resolution of a spectroscopic survey. This would not only significantly increase the computational effort, but also lead to large shot noise in the auto-correlation spectra². These are the main reasons that, for spectroscopic surveys, the correlation function is a more promising summary statistic than the angular power spectrum; we therefore need to determine the impact of lensing magnification on this statistic.

2.2. The two-point correlation function

In spectroscopic surveys, there are two standard estimators used to extract information from galaxy number counts: the 2PCF, and its Fourier transform, the power spectrum. In this paper, we concentrate on the correlation function, since lensing magnification can be included in this estimator in a straightforward way. This is not the case for the power spectrum, which requires non-trivial extensions to account for magnification (see e.g. Grimm et al. 2020; Castorina & di Dio 2022).

The 2PCF can be calculated in the curved-sky, that is, without assuming that the two directions n and n′ are parallel:

$$\begin{array}{c}\hfill \xi (\mathit{n},z,{\mathit{n}}^{\prime },z\prime )=⟨\mathrm{\Delta }(\mathit{n},z)\mathrm{\Delta }({\mathit{n}}^{\prime },z^{\prime })⟩.\end{array}$$(6)

The curved-sky density and RSD contributions were first derived in Szalay et al. (1998) and Szapudi (2004). This method, which can be straightforwardly applied to any local contribution of Δ(n, z), cannot be used to calculate contributions from integrated effects like lensing magnification, which is the main subject of this work. The magnification contribution was calculated in Tansella et al. (2018a) using an alternative method proposed in Campagne et al. (2017). The detailed expressions for all contributions can be found in Tansella et al. (2018a); for completeness, we repeat them in Appendix E.1. The expressions for the 2PCF significantly simplify in the flat-sky approximation: in this approximation, one assumes that the two directions n and n′ are parallel, and one neglects the redshift evolution of Δ. In this case, the density and RSD contributions (hereafter called standard terms) in a “thick” redshift bin with mean redshift $\overline{z}$ take the following simple form:

$$\begin{array}{c}\hfill {\xi }^{\mathrm{std}}(d,\overline{z},\mu )={\xi }_0^{\mathrm{std}}(d,\overline{z})+{\xi }_2^{\mathrm{std}}(d,\overline{z})L_2(\mu )+{\xi }_4^{\mathrm{std}}(d,\overline{z})L_4(\mu ),\end{array}$$(7)

where d denotes the comoving separation between the correlated volume elements or “voxels”, $\overline{z}$ is the centre of the bin interval in which the correlation function is measured, μ is the cosine of the angle between the direction of observation n and the vector connecting the voxels, and L_ℓ denotes the Legendre polynomial of order ℓ. The standard multipoles, ${\xi }_{\mathcal{l}}^{\text{std}}$, are given by the well-known expressions

$$\begin{array}{cc}& {\xi }_0^{\mathrm{std}}(d,\overline{z})=[b^2(\overline{z})+\frac{2}{3}b(\overline{z})f(\overline{z})+\frac{1}{5}{f}^2(\overline{z})]{\mu }_0(d,\overline{z}),\hfill \\ \hfill & {\xi }_2^{\mathrm{std}}(d,\overline{z})=-[\frac{4}{3}f(\overline{z})b(\overline{z})+\frac{4}{7}{f}^2(\overline{z})]{\mu }_2(d,\overline{z}),\hfill \\ & {\xi }_4^{\mathrm{std}}(d,\overline{z})=\frac{8}{35}{f}^2(\overline{z}){\mu }_4(d,\overline{z}),\hfill \end{array}$$(8)

where

$$\begin{array}{c}\hfill f(\overline{z}):=\frac{\mathrm{d}ln\delta }{\mathrm{d}lna},\end{array}$$(9)

is the growth rate of structure. The functions ${\mu }_{\ell }(d,\overline{z})$ are given by

$$\begin{array}{c}\hfill {\mu }_{\ell }(d,\overline{z})=\frac{1}{2{\pi }^2}{\displaystyle {\int }_0^{\infty }\mathrm{d}k{k}^2{P}_{\delta \delta }(k,\overline{z})j_{\ell }(kd),}\end{array}$$(10)

where j_ℓ denotes the spherical Bessel function of order ℓ, and $P_{\delta \delta }(k,\overline{z})$ is the linear matter power spectrum in the comoving gauge. We see that, in the flat-sky approximation, density and RSD are fully encoded in the first three even multipoles. The expressions for computing higher-order multipoles without using the flat-sky approximation can be found in Appendix E.1.

The magnification contribution can also be simplified using the flat-sky approximation and the Limber approximation. This has been derived in detail in Tansella et al. (2018a). It contains an infinite series of multipoles,

$$\begin{array}{c}\hfill {\xi }^{\mathrm{magn}}(d,\overline{z},\mu )={\displaystyle \sum _{\ell }{\xi }_{\ell }^{\mathrm{magn}}(d,\overline{z})L_{\ell }(\mu ),}\end{array}$$(11)

with

$$\begin{array}{cc}\hfill {\xi }_{\ell }^{\mathrm{magn}}(d,\overline{z})& =\frac{2\ell +1}{2}\{\frac{3{\Omega }_{\mathrm{m},0}}{2\pi }(1+\overline{z})b(\overline{z})(5s-2)d{\displaystyle {\int }_0^1}\mathrm{d}\mu \mu L_{\ell }(\mu )\hfill \\ \hfill & \times {\displaystyle {\int }_0^{\infty }}\mathrm{d}{k}_{\perp }{k}_{\perp }\frac{H_0^2}{c^2}{P}_{\delta \delta }(k_{\perp },\overline{z})J_0\left(k_{\perp }d\sqrt{1-{\mu }^2}\right)\hfill \\ \hfill & +\frac{9{{\Omega }_{\mathrm{m},0}}^2}{4\pi }{(5s-2)}^2{\displaystyle {\int }_0^{\overline{r}}}\mathrm{d}{r}^{\prime }\frac{{(\overline{r}-r^{\prime })}^{2}r{\prime }^2}{{\overline{r}}^2{a}^2(r^{\prime })}{\displaystyle {\int }_0^1}\mathrm{d}\mu L_{\ell }(\mu )\hfill \\ \hfill & \times {\displaystyle {\int }_0^{\infty }}\mathrm{d}{k}_{\perp }{k}_{\perp }\frac{H_0^4}{c^4}{P}_{\delta \delta }(k_{\perp },z(r^{\prime }))J_0\left(k_{\perp }\frac{r^{\prime }}{\overline{r}}d\sqrt{1-{\mu }^2}\right)\}.\hfill \end{array}$$(12)

Here H₀ = 100 h km s⁻¹ Mpc⁻¹, k_⊥ = k − kμn is the projection of the Fourier space wavevector k to the direction normal to n, and J₀ denotes the Bessel function of zeroth order. The first two lines of Eq. (12) contain the density–magnification correlation: when computing the multipoles, one averages over all orientations of the pair of voxels. For each orientation, the galaxy that is farther away is lensed by the one in the foreground. This effect clearly has a non-trivial dependence on the orientation angle, μ, which enters in the argument of the Bessel function J₀. The last two lines contain the magnification–magnification correlation, due to the fact that both galaxies are lensed by the same foreground inhomogeneities. In all the terms, the functions are evaluated at the mean redshift of the bin $\overline{z}$, $\overline{r}$ denotes the comoving distance at that redshift, and a(r) denotes the scale factor evaluated at comoving distance r. For completeness, in Appendix E.2 we list the semi-analytic expressions that allowed us to efficiently evaluate the flat-sky magnification terms.

In Fig. 1 we show a comparison between the curved-sky expression and the flat-sky approximation for the standard multipoles (left panel) and the magnification multipoles (right panel), in one of the redshift bins of Euclid, $\overline{z}=1.4$. We checked that most of the constraining power comes from standard terms of the monopole and quadrupole below d = 150 Mpc, where the difference between the curved-sky and the flat-sky expressions is less than 0.2% (it reaches 0.7% for the hexadecapole). Similar results are obtained for the other redshift bins; hence, using the flat-sky approximation is very well justified. For the magnification contribution to the monopole, we see that the flat-sky approximation differs from the curved-sky result already at small separation by roughly 5%. However, since the magnification is a sub-dominant contamination to the total signal, a 5% error is perfectly acceptable. We note that the difference between curved-sky and flat-sky actually increases for small separations; as shown in Jelic-Cizmek (2021), this can occur when the dominant contribution to the multipoles is the density-lensing term, for which the accuracy of the flat-sky approximation becomes progressively worse at smaller scales.

Fig. 1. Impact of the flat-sky approximation on the multipoles of the 2PCF. Top: Curved-sky (solid) vs. flat-sky (dashed) multipoles, with contributions from standard terms (left) and from just lensing magnification (right). Lensing magnification is computed using the values of the local count slope given in Table 2, which are the fiducial values assumed in our analysis. Bottom: Their relative difference in percent, taking the curved-sky case as the reference value.

In Fig. 2 we compare the 2PCF with and without lensing magnification. On small scales, lensing magnification is not very important. However, for scales of more than 300 Mpc it can contribute up to 20% to the monopole and the quadrupole and up to 50% or more to the hexadecapole. We see that the further apart the galaxies are, the more significant is the contribution from lensing magnification to their correlation. This is due to the fact that the density correlations quickly decrease with separation, while the lensing magnification correlations do not. The magnification–magnification correlations are indeed integrated all the way from the sources to the observer, and therefore contain contributions from small scales, when the two lines of sight are close to the observer.

Fig. 2. Impact of magnificatiom on the multipoles of the 2PCF. Top: Comparison of multipoles without (solid) and with (dashed) lensing magnification. Bottom: Their relative difference in percent, taking the case with magnification as the reference value.

We next assessed the impact of these multipoles on the analysis of data from Euclid. In particular, we determined the following: first, if the contribution from magnification can improve our measurement of cosmological parameters; and second, to what extent neglecting magnification in the analysis will shift the best-fit value of the parameters, consequently biasing the analysis.

We considered two cases. In the first case, we fixed the cosmological model, and we studied how magnification impacts the parameters of this model. For this case, we studied two models: a minimal ΛCDM model and a dynamical dark energy model. In the second case, we performed a model-independent analysis, that is, we rewrote the standard multipoles in terms of the power spectrum at z_*, using that $P_{\delta \delta }(\overline{z})=P_{\delta \delta }(z_{\ast }){\sigma }_8^2(\overline{z})/{\sigma }_8^2(z_{\ast })$. We chose z_* to be well within in the matter-dominated era before acceleration started. We assumed that, at z_*, general relativity is valid, and that the power spectrum is therefore fully determined by the early Universe parameters that have been measured by the cosmic microwave background (CMB). With this, the functions μ_ℓ depend only on the power spectrum at z_*, while the evolution from z_* to $\overline{z}$ is fully encoded in two functions:

$$\begin{array}{c}\hfill \stackrel{\sim }{f}(\overline{z})=f(\overline{z}){\sigma }_8(\overline{z})\mathrm{and}\stackrel{\sim }{b}(\overline{z})=b(\overline{z}){\sigma }_8(\overline{z}).\end{array}$$(13)

We obtained the multipoles of the standard terms:

$$\begin{array}{cc}& {\xi }_0^{\mathrm{std}}(d,\overline{z})=[{\stackrel{\sim }{b}}^2(\overline{z})+\frac{2}{3}\stackrel{\sim }{b}(\overline{z})\stackrel{\sim }{f}(\overline{z})+\frac{1}{5}{\stackrel{\sim }{f}}^2(\overline{z})]\frac{{\mu }_0(d,z_{\ast })}{{\sigma }_8^2(z_{\ast })},\hfill \\ \hfill & {\xi }_2^{\mathrm{std}}(d,\overline{z})=-[\frac{4}{3}\stackrel{\sim }{f}(\overline{z})\stackrel{\sim }{b}(\overline{z})+\frac{4}{7}{\stackrel{\sim }{f}}^2(\overline{z})]\frac{{\mu }_2(d,z_{\ast })}{{\sigma }_8^2(z_{\ast })},\hfill \\ & {\xi }_4^{\mathrm{std}}(d,\overline{z})=\frac{8}{35}{\stackrel{\sim }{f}}^2(\overline{z})\frac{{\mu }_4(d,z_{\ast })}{{\sigma }_8^2(z_{\ast })}·\hfill \end{array}$$(14)

In this case, the functions μ_ℓ(d, z_*) are considered fixed since they are very well determined by CMB measurements, and the functions $\stackrel{\sim }{f}$ and $\stackrel{\sim }{b}$ are two free functions that depend on the mean redshift $\overline{z}$ of the bins (see also Jelic-Cizmek et al. 2021 for an introduction of this method). These two free functions fully encode any deviations from general relativity at late time. The only approximation that enters here is that we neglected the k dependence of the growth of density, meaning that σ₈ and f depend only on redshift. This assumption can easily be relaxed: it slightly complicates the analysis, since $\stackrel{\sim }{f}$ and $\stackrel{\sim }{b}$ would have to be taken inside the integrals in Eq. (10), but it does not change the procedure.

This model-independent analysis is one of the key goals of the Euclid spectroscopic survey. It is very powerful, since it allows us to measure the growth rate of structure without assuming a particular model of gravity or dark energy. This growth rate can then be compared with the predictions of any model beyond ΛCDM. In the following, we determined how neglecting magnification in the analysis could shift the best-fit values of $\stackrel{\sim }{b}$ and $\stackrel{\sim }{f}$ in each redshift bin. We note that here for simplicity we fixed the cosmological parameters that determine the functions μ_ℓ at early time, z_*, to their fiducial value extracted from Planck data (Planck Collaboration XIV 2016). In practice, one can also let these parameters vary and perform a combined analysis with CMB data.

It is worth mentioning that, by construction, the analysis using multipoles of the correlation function does not account for correlations between different redshift bins: the correlation function is averaged over directions within a given bin, and each redshift bin is considered to be independent. However, the main motivation of measuring the multipoles of the correlation function is to extract the growth rate f, which is encoded in the peculiar velocities of galaxies within linear perturbation theory. Moreover, since the correlations of peculiar velocities quickly decrease with separation, one does not lose a significant amount of information by neglecting cross-correlation between bins. The situation is of course different for magnification, which, as it is an integral along the line of sight, is strongly correlated between different bins. Therefore, we expect that neglecting cross-correlations of different bins will strongly reduce the magnification signal, compared to the angular power spectra used in the analysis of the photometric sample, which is able to account for correlations between the bins (Euclid Collaboration 2022). As such, the shift induced by neglecting lensing magnification is expected to be smaller in the spectroscopic analysis than in the photometric one.

3. Euclid specifics from the Flagship simulation

In order to calculate the linear galaxy bias and local count slope observables for this analysis, which we employed as our fiducial values for the Fisher and MCMC analyses, we used Flagship v1.8.4 galaxy mock samples, whose redshift distribution of the number density, 𝒩(z), is split into 13 equally spaced bins (in redshift), between z = 0.9 and z = 1.8, in real space. As we see from Fig. 3, this allowed us to accurately capture the redshift evolution of the bias and the local count slope, while still having enough galaxies in each bin to obtain a precise measurement. We imposed a cut in the Hα flux F_Hα (in units of erg s⁻¹ cm⁻²), log₁₀[F_Hα/(1 erg s⁻¹ cm⁻²)] > −15.7, that can be transformed to the corresponding AB magnitude limit, m_* = m_AB < −15.75³. The linear galaxy bias and local count slope, estimated from the Flagship simulation as described in this section, are assumed in our analysis to be the “true” values for the Euclid spectroscopic sample, and we used them throughout the paper as fiducial values.

Fig. 3. Galaxy bias (top panel) and the local count slope (lower panel) with linear interpolation (dashed), as used in our analysis, along with their associated error bars, as well as a polynomial fit (solid). For the exact numerical values of the coefficients, refer to Eq. (B.3).

3.1. Linear galaxy bias

The linear galaxy bias is obtained by fitting the curved-sky angular power spectrum of the data to the corresponding fiducial prediction for 𝒩(z) in each redshift bin. The angular power spectrum is extracted using Polspice⁴ with a mask to generate 100 jackknife regions that we used to calculate the covariance matrix. This masked region corresponds to the area of the sky that the Cosmohub’s Flagship v1.8.4 release (Tallada et al. 2020; Carretero et al. 2017) does not cover, which corresponds to seven-eighths of the total sky. The mask is also generated with Polspice by masking the pixels outside the region 0 ° < RA < 90° and 0 ° < Dec < 90°. Then when considering the jackknife regions, we used a k-means clustering algorithm to select 100 regions with roughly the same number of pixels (since each HEALPix⁵ pixel covers the same sky area) inside the unmasked region. For each jackknife resampling, we included in the mask the pixels of a different region in order to exclude it from that jackknife iteration C_ℓ calculation. The prediction is determined using CCL⁶ (Chisari et al. 2019) with the fiducial cosmological parameter values being those used in the Flagship v1.8.4 simulation (see Table 1). The linear scales considered range from ℓ = 50 to an ℓ_max that increases as we go to higher redshifts (as the effective non-linear galaxy bias scale shifts to higher multipoles, i.e. smaller angular scales), starting at ℓ_max ∼ 300 for z = 0.9 to ℓ_max ∼ 500 for z = 1.8. Since CCL does not yet allow us to perform calculations without the Limber approximation, we employed it for all of the scales used to estimate the linear galaxy bias. We set the minimum multipole to ℓ = 50, which is a rather conservative limit in order to avoid large Limber approximation deviations from the theory at any redshift considered for this analysis. On top of that, scales with ℓ < 50 usually have very high error bars due to sample variance so they can be ignored since their statistical weight to estimate the galaxy bias is very low. The maximum scale for each redshift is estimated by comparing the relative ratio between the linear and non-linear matter power spectrum prediction and setting a maximum relative difference of 3%. This maximum difference should be good enough since galaxy clustering is known to follow linear predictions down to smaller scales than dark matter. At the scales used, the Limber approximation adds up to a 4% variation of the predicted angular power spectrum, which translates into an error below 2% on the galaxy bias estimation. The χ² distribution is then calculated for different values of the galaxy bias in relation to the square ratio of the data C_ℓ’s to the prediction C_ℓ’s using only the diagonal values of the jackknife covariance matrix. We estimated the linear galaxy bias as the minimum value of the χ² distribution, and we set the error to the 1σ variance.

3.2. Local count slope

To measure the local count slope from the Flagship catalogues, we computed, in each of the 13 redshift bins, the cumulative number of galaxies $\overline{N}(z,m<m_{\ast })$ at the magnitude limits m_* and m_* ± 0.04, and we computed the logarithmic derivative to obtain s(z, m_*) through Eq. (5). We also generated 100 jackknife regions in order to calculate the variance of the results. In Fig. 3, we see that s increases with redshift; this is due to the fact that a fixed apparent magnitude threshold, m_*, corresponds to a larger intrinsic luminosity threshold L_* at high redshift than at low redshift. This is due to the fact that the slope in m of the Schechter luminosity function, which is assumed here, increases with redshift.

4. Method

This study employed two complementary approaches for forecasting the constraining capabilities of future Euclid data: the Fisher matrix formalism and the MCMC method. The details of each approach are outlined in Sects. 4.1, 4.2, and 4.3, respectively. We used the code coffe⁷, which has been validated against the code CosmoBolognaLib⁸ (for details of the validation, see Appendix A), to compute the multipoles of the 2PCF.

4.1. The Fisher matrix formalism: Cosmological constraints

The Fisher matrix can be defined as the expectation value of the second derivatives of the logarithm of the likelihood under study with respect to the parameters of the model (see e.g. Tegmark 1997):

$$\begin{array}{c}\hfill {\mathsf{F}}_{\alpha \beta }=\langle -\frac{{\partial }^2\ln \mathcal{L}}{\partial {\theta }_{\alpha }\partial {\theta }_{\beta }}\rangle ,\end{array}$$(15)

where α and β label two model parameters θ_α and θ_β.

In the particular case of Gaussian-distributed data, the Fisher matrix is given by

$$\begin{array}{c}\hfill {\mathsf{F}}_{\alpha \beta }=\frac{1}{2}\mathrm{Tr}\left[\frac{\partial \mathsf{C}}{\partial {\theta }_{\alpha }}{\mathsf{C}}^{-1}\frac{\partial \mathsf{C}}{\partial {\theta }_{\beta }}{\mathsf{C}}^{-1}\right]+{\displaystyle \sum _{\mathit{pq}}\frac{\partial {\mathsf{D}}_p}{\partial {\theta }_{\alpha }}{\left({\mathsf{C}}^{-1}\right)}_{\mathit{pq}}\frac{\partial {\mathsf{D}}_q}{\partial {\theta }_{\beta }},}\end{array}$$(16)

where D represents the mean of the data vector and C is the covariance matrix of the data. The trace, Tr, and the sum over the indexes, p and q, stand for the summation over the different elements of the data vector.

In the present analysis we considered the 2PCF as our main observable. Given the set of model parameters {θ_α}, the Fisher matrix for the multipoles of the 2PCF measured in a bin centred in ${\overline{z}}_i$ is

$$\begin{array}{c}\hfill {\mathsf{F}}_{\alpha \beta }^{\mathrm{bin}}({\overline{z}}_i)={\displaystyle \sum _{\mathit{jk}}\sum _{\ell m}\frac{\partial {\xi }_{\ell }(d_j,{\overline{z}}_i)}{\partial {\theta }_{\alpha }}{\mathsf{C}}^{-1}[{\xi }_{\ell }^j,{\xi }_m^k]({\overline{z}}_i)\frac{\partial {\xi }_m(d_k,{\overline{z}}_i)}{\partial {\theta }_{\beta }},}\end{array}$$(17)

where the sum runs over the voxel separations {d_j, d_k} as well as the even multipoles ℓ, m = 0, 2, 4 and ${\xi }_{\mathcal{l}}^j\equiv {\xi }_{\mathcal{l}}(d_j)$. We note that (angular) power spectra observables follow a Wishart distribution if fluctuations are Gaussian. In this case, the Fisher analysis gives a better approximation if we consider only the second term in Eq. (16) (see e.g. Carron 2013; Bellomo et al. 2020). In the following we assume that the same is true for the multipoles of the correlation function. The binned covariance of the 2PCF multipoles at mean redshift ${\overline{z}}_i$, denoted with C, is computed following the Gaussian theoretical model described in Grieb et al. (2016) and Hall & Bonvin (2017). The cosmic variance contribution includes only the density and RSDs, while magnification is neglected. This is a good approximation, since the covariance is a four-point function, which contains a sum over all possible separations between pairs of pixels. Hence, even at a large separation, the covariance is dominated by correlations at small scales, where the density and RSD strongly dominate over magnification. The shot noise contribution is estimated from the number densities reported in Euclid Collaboration (2020; henceforth referred to as EP:VII), Table 3. The full expression for the covariance can be found in Appendix C.

Following EP:VII, we neglected the cross-correlations between redshift bins. Thus, the full Fisher matrix is

$$\begin{array}{c}\hfill {\mathsf{F}}_{\alpha \beta }={\displaystyle \sum _{{\overline{z}}_i}{\mathsf{F}}_{\alpha \beta }^{\mathrm{bin}}({\overline{z}}_i).}\end{array}$$(18)

This approximation is justified by the high precision of spectroscopic redshift estimates, which leads to essentially no overlap between redshift bins. In the case of a photometric analysis, cross-correlations between different bins can provide significant information (see e.g. Tutusaus et al. 2020; Euclid Collaboration 2022). The marginalised 1σ errors on the cosmological parameters can then be estimated from the Cramér–Rao bound:

$$\begin{array}{c}\hfill {\sigma }_{\alpha }=\sqrt{{({\mathsf{F}}^{-1})}_{\alpha \alpha }}.\end{array}$$(19)

It is important to mention that, although the Fisher matrix formalism is a powerful forecasting tool, some limitations do exist. A Fisher forecast uses a Gaussian approximation by construction, which can differ from the true posterior if the data are not constraining enough. Furthermore, the signal and covariance may have a strong non-linear dependence on the parameters {θ_α}, in which case the Fisher matrix does not capture all of the information about the likelihood. In order to validate the results of our Fisher formalism, we also performed, for one of the cases, a MCMC analysis to properly sample the posterior of the parameters (see Sect. 4.3). As we will see, we find that the relevance of lensing magnification is well captured by a Fisher forecast. Another drawback of the Fisher formalism worth mentioning is that it only provides forecast uncertainties around a fiducial model. In this analysis we are also interested in the bias on the posteriors because of wrong model assumptions (neglecting magnification). The standard Fisher formalism prevents us from doing this study, but extensions to the formalism can be considered, as described in Sect. 4.2. All Fisher forecasts were computed with the Python package FITK⁹.

4.2. The Fisher matrix formalism: Bias on parameter estimation

The Fisher matrix formalism described above allows us to quantify the gain or loss in constraining power, when magnification is included in the theoretical model for the observed multipoles of the 2PCF. This study can be carried out by simply comparing the Fisher matrix in Eq. (18) and the corresponding marginalised constraints when magnification is neglected or included in the analysis.

Another, actually more important question to address is whether neglecting magnification leads to significant biases (shifts) in the inferred cosmological parameters. In order to answer this question, we followed the approach described in, for example, Taylor et al. (2007) and widely adopted in the literature (see Kitching et al. 2009; Camera et al. 2015; Di Dio et al. 2016; Cardona et al. 2016; Lepori et al. 2020; Jelic-Cizmek et al. 2021; Euclid Collaboration 2022). We extended the parameter space to include the amplitude of magnification, ϵ_L. We can explicitly write the dependence of our model on ϵ_L as follows:

$$\begin{array}{c}\hfill {\xi }_{\ell }(d_j,{\overline{z}}_i)={\xi }_{\ell }^{\mathrm{std}}(d_j,{\overline{z}}_i)+{ϵ}_{\mathrm{L}}{\xi }_{\ell }^{\mathrm{magn}}(d_j,{\overline{z}}_i),\end{array}$$(20)

where ${\xi }_{\mathcal{l}}^{\text{std}}$ represents the standard contributions of density and RSDs, and ${\xi }_{\mathcal{l}}^{\text{magn}}$ is the magnification contribution, which includes the terms magnification × magnification and the cross-correlation of magnification with density and the RSD. The amplitude ϵ_L is not a free parameter but rather a fixed one that is set to either 0 (in a “wrong” model that neglects magnification) or 1 (in a “correct” model that consistently includes magnification). The “wrong” and “correct” models share a common set of parameters, {θ_α}, and their estimation will be biased in the wrong model as a result of the shift in the fixed parameter ϵ_L. Using a Taylor expansion of the likelihood around the wrong model, and truncating the series at the linear order, we obtain the following formula for the biases:

$$\begin{array}{c}\hfill \mathrm{\Delta }({\theta }_{\alpha })={\displaystyle \sum _{\beta }{\left({\mathsf{F}}^{-1}\right)}_{\alpha \beta }{B}_{\beta },}\end{array}$$(21)

where F is the Fisher matrix of the common set of parameters evaluated for the wrong model, and

$$\begin{array}{c}\hfill B_{\beta }={\displaystyle \sum _i\sum _{\mathit{jk}}\sum _{\ell m}{\xi }_{\ell }^{\mathrm{magn}}(d_j,{\overline{z}}_i){\mathsf{C}}^{-1}[{\xi }_{\ell }^j,{\xi }_m^k]({\overline{z}}_i)\frac{\partial {\xi }_{\ell }^{\mathrm{std}}(d_k,{\overline{z}}_i)}{\partial {\theta }_{\beta }}·}\end{array}$$(22)

Equation (22) implicitly assumes that magnification constitutes a small contribution to the observable; therefore, the outcome can be quantitatively trusted only when small values of the biases are found. Nevertheless, large biases are a clear indication that the lensing magnification significantly contributes to the observable and that the starting hypothesis should be rejected. Therefore, it is a good diagnostic to assess whether magnification can be neglected or should be modelled in the analysis.

4.3. Markov chain Monte Carlo

The MCMC method is a standard statistical technique with which we numerically sample the posterior probability starting from a prior probability and assuming a likelihood function (the probability of the data given the hypothesis). An excellent description of the method can be found in Verde (2007).

In our analysis, under the assumption that our data are Gaussian, we sampled the posterior of the following likelihood:

$$\begin{array}{c}\hfill ln(\mathcal{L})=-\frac{1}{2}\mathrm{\Delta }{\mathit{\xi }}^T{\mathsf{C}}^{-1}\mathrm{\Delta }\mathit{\xi },\end{array}$$(23)

where C is the covariance matrix, and Δξ is a vector whose elements are given by

$$\begin{array}{c}\hfill \mathrm{\Delta }{\mathit{\xi }}_{[\ell ,i,j]}={\xi }_{\ell }{(d_j,{\overline{z}}_i)}^{\mathrm{ref}}-{\xi }_{\ell }(d_j,{\overline{z}}_i),\end{array}$$(24)

where the first term is part of a synthetic dataset computed previously used as our “reference” or “fiducial model”, while the second term is computed at each steps of the MCMC, varying the value of free parameters inside the parameters space described by the prior function. To simplify the analysis, we neglected the dependence of the covariance matrix on cosmological parameters, which we fixed to their reference values.

We assumed a flat prior density for each parameter. In order to speed up the convergence of our chains, we assumed as free parameters (with flat priors) ω_m, 0 and ω_b, 0 instead of Ω_m, 0 and Ω_b, 0, and subsequently re-parametrised the chain using the relation ω_i, 0 = Ω_i, 0h².

We used the Python package emcee¹⁰ (Foreman-Mackey et al. 2013) to implement the MCMC. Our sampler was composed of 32 walkers, and we used the “stretch move” ensemble method described in Goodman & Weare (2010). Each walker generates a chain with a number of steps of the order 10⁵ before converging. Our MCMC code is run in parallel using the Python package schwimmbad (Price-Whelan & Foreman-Mackey 2017).

In our analysis we also discarded a number of points as burn-in, given by twice the maximal integrated auto-correlation time, τ, of all the parameters (Goodman & Weare 2010)¹¹. The results of the sampling were then analysed with the Python package GetDist (Lewis 2019).

5. Results

We assessed the impact of neglecting magnification for three different cases: a minimal ΛCDM model, a dynamical dark energy model, and a model-independent analysis measuring the bias and growth rate. For each case, we computed the change in the constraints due to including magnification and the shift in the parameters due to neglecting magnification.

To be consistent with the Flagship simulation, the fiducial cosmology adopted in our analysis is a flat ΛCDM model with no massive neutrino species. The set of parameters varied in the analysis comprises: the present matter and baryon density parameters, respectively Ω_m, 0 and Ω_b, 0; the dimensionless Hubble parameter h; the amplitude of the linear density fluctuations within a sphere of radius 8 h⁻¹ Mpc at present time, σ₈; the spectral index of the primordial matter power spectrum n_s; and the equation of state for the dark energy component {w₀, w_a}, which parametrise the time evolution of the dark energy equation of state parameter as

$$\begin{array}{c}\hfill w(z)=w_0+w_a\frac{z}{1+z}·\end{array}$$(25)

This model is also known as Chevallier–Polarski–Linder parametrisation (Chevallier & Polarski 2001; Linder 2003). The fiducial values of the cosmological parameters used in the analysis are reported in Table 1. They correspond to the w₀w_aCDM parameters used in the Euclid Flagship simulation (see Sect. 3).

In addition to these cosmological parameters, we introduced nuisance parameters and marginalised over them; in particular, the linear galaxy bias in each redshift bin, {b_i}, i = 1, …, N_bins, are included as nuisance parameters. We modelled them as constant within each redshift bin, and we estimated their fiducial values using the Flagship simulation, v1.8.4, as described in Sect. 3. We list the values of the nuisance parameters used in each redshift bin as well as the expected density of emitters in Table 2. The impact of magnification on the cosmological parameters may depend on the model chosen to describe our Universe. We therefore run our analysis for two different cosmological models and comment on the difference between the results when relevant. We considered (i) a minimal flat ΛCDM model with five free parameters, {Ω_m, 0, Ω_b, 0, h, n_s, σ₈} plus nuisance parameters., and (ii) a flat dynamical dark energy model with seven free parameters, {Ω_m, 0, Ω_b, 0, w₀, w_a, h, n_s, σ₈} plus nuisance parameters.

We also included a cosmology-independent analysis, where as free parameters we considered the modified galaxy bias and the modified growth rate in each redshift bin, $\stackrel{\sim }{f}(z)$ and $\stackrel{\sim }{b}(z)$, defined in Eq. (13). To be conservative, as well as to reduce the impact of non-linearities, we only considered separations between d_min = 40 Mpc and d_max = 385 Mpc in each redshift bin, and, unless specified otherwise, multipoles ℓ ∈ {0, 2, 4}. We used voxels of size L_p = 5 Mpc. We checked that reducing them further to 2.5 Mpc does not improve the constraints anymore, due to shot noise, which saturates the signal-to-noise ratio for too small voxel sizes.

5.1. Lensing magnification signal-to-noise

As an estimate of the impact of lensing magnification, we first computed the signal-to-noise ratio (S/N) of the lensing contribution to the multipoles of the 2PCF, which we define as

$$\begin{array}{c}\hfill S/N({\overline{z}}_i):={\displaystyle \sqrt{\sum _{j,k,\ell ,m}{\xi }_{\ell }^{\mathrm{magn}}(d_j,{\overline{z}}_i){\mathsf{C}}^{-1}[{\xi }_{\ell }^j,{\xi }_m^k]({\overline{z}}_i){\xi }_m^{\mathrm{magn}}(d_k,{\overline{z}}_i)},}\end{array}$$(26)

where the sum goes over all pairs of voxels and all even multipoles taken into consideration. The results are shown in Fig. 4 and Table 3. Since lensing magnification also contributes to multipoles larger than the hexadecapole, we show the S/N for two cases: ℓ_m = 4 and ℓ_m = 6, where ℓ_m denotes the highest multipole used in the analysis. As we can see, the S/N is smallest in the lowest redshift bin, and increases as we go to higher redshifts. This is a consequence of two effects: first, the local count slope for Euclid, s(z), increases with increasing redshift (see Table 2 and Fig. 3)¹². Second, the lensing magnification term is an integrated effect and, as such, has the largest impact at high redshifts.

Fig. 4. S/N of the lensing magnification for two different scenarios, with dmin = 40 Mpc. The horizontal bars denote the widths of the redshift bins.

5.2. Full-shape cosmological analysis

In this section, we focus on the impact of magnification on the full-shape cosmological analysis of the 2PCF, for the Euclid spectroscopic sample. Magnification in principle affects both the best-fit estimation of cosmological parameters and their constraints. In order to quantify the relevance of the effect, we adopted the Fisher formalism described in Sect. 4.1, and we validated the results by comparing them to the outcome of a full MCMC analysis. In order to estimate the impact of magnification on constraints on cosmological parameters, we run two Fisher analyses, one with and one without lensing magnification, from which we estimate the marginalised 1σ errors, and we compared the values in the two cases.

As lensing magnification contains additional independent information, we expect the constraints to improve slightly when including it. In Table 4 we report the improvement in constraints on cosmological parameters when magnification is included in the analysis. For a ΛCDM model, the impact of magnification on the reduction of the error bars is ≲5% for all cosmological parameters. In the dynamical dark energy model w₀w_aCDM, the impact of magnification is slightly larger, as the 1σ error bars in w₀ and w_a are reduced by about 10%. Nevertheless, the improvement due to magnification in the constraining power decreases for other parameters. It is important to note that in this test we are assuming the values of the local count slope to be exactly known. While it is in principle possible to estimate s(z) independently of the cosmological analysis, from the slope of the luminosity distribution of the galaxy sample, this measurement will be affected by several systematics (see for example Hildebrandt 2016). Therefore, we also considered a more pessimistic scenario where we assumed no prior knowledge on the local count slopes and thus we marginalise over the values of s in each redshift bins. In Table 4 we also compared the constraints on cosmological parameters for a model that neglects magnification, and a model that includes the effect, assuming no information on the local counts slope. In this pessimistic setting, the constraints obtained when magnification is included are worse than the ones obtained when magnification is neglected. This reduction in constraining power as compared to a model without magnification is up to 6% for ΛCDM and becomes up to 10% for the w₀w_aCDM parametrisation. This is due to the fact that magnification does not contribute very significantly to the cosmological information that can be extracted from the 2PCF, while the extra nuisance parameters introduced in this model slightly increase the degeneracy between the other parameters included in the analysis. We also note that, in the dynamical dark energy model, magnification mostly affects the cosmological constraints of σ₈, w₀ and w_a, while the remaining model parameters are substantially unaffected. Nevertheless, the impact of magnification on the constraints of cosmological parameters is small; including it leads to changes of at most ±10% in the error bars of parameters within the full-shape analysis.

In Fig. 5 we show a visual comparison of the constraints for the three cases discussed above. We validated these results by running an MCMC analysis for the ΛCDM model, including lensing magnification in the analysis. A direct comparison of the contour plots for the Fisher and MCMC methods can be found in Fig. 6 and Table 4. The Fisher results are accurate at the 20% level (see Tables F.1 and F.2 for the actual values). By comparing the results from the Fisher matrix and MCMC analyses in Fig. 6, it is easy to grasp the impact of the non-Gaussianity of the posterior. Such a non-Gaussianity, especially in the form of a skewness of the distribution, leads to a slight variation of the constraining power w/ and w/o magnification. In particular, this can be appreciated by comparing the last row of Table F.2 to the third one (ΛCDM only), which also causes some changes in sign in some of the constraining power variation (e.g. first and second to last lines of Table 4): the change in the 68% C.L.’s on those parameters switches sign. We note that all of these variations are very small – at the percentage level – and do not affect our main conclusions.

Fig. 5. Comparison of 68% C.L.’s obtained from the Fisher analysis for ΛCDM (left) and w0waCDM (right) with no lensing magnification (blue), with lensing magnification and the local count slope fixed (orange) in each redshift bin, and with magnification and local count slope marginalised (green). We note that for wa, we do not divide by the fiducial as it is zero, and instead we show the absolute error. For the corresponding data, see Table 4.

Fig. 6. 68% (inner) and 95% (outer) 2D confidence regions and 1D posteriors for ΛCDM with marginalisation over galaxy biases, using MCMC (blue) and the Fisher analysis (orange), with only contributions from just standard terms (dashed) and with standard terms plus lensing magnification (solid). The dashed black lines denote fiducial values of the cosmological parameters. The label “no magn” refers to an analysis that does not include magnification in the model despite the effect being present in the data. The label “w/magn” incorporates magnification both in the model and in the data.

We also investigated the effect of magnification on the accuracy of best-fit estimation of cosmological parameters. To accomplish this, we employed two distinct techniques, namely, Fisher analysis and MCMC analysis. In the Fisher analysis, we computed the biases in the best-fit estimates by employing Eq. (21). In the MCMC analysis, we generated synthetic data based on a model that takes magnification into account and fitted the data using two theoretical predictions: one that includes magnification, and the other that neglects it. The differences in the best-fit parameters in these two cases provide the shifts induced by ignoring magnification in our modelling. In Table 4 (bottom block), we report the values of the shifts obtained with the Fisher analysis, for the ΛCDM and w₀w_aCDM parametrisations. In both cases, we find shifts below 1σ. For the ΛCDM analysis, the best-fit estimate is biased at the level of ∼0.5–0.7σ for all cosmological parameters. The impact is less relevant in the w₀w_aCDM model, mainly due to the worse constraints on cosmological parameters. The largest shifts in this case are found for n_s (∼0.5σ) and σ₈ (∼0.2σ). In Table 4 (bottom block) we also report the MCMC result, generated only for ΛCDM, where the shifts are larger. We find that while the shift found in the MCMC analysis is in most cases slightly larger, Fisher and MCMC forecast give consistent values of the shifts, both in terms of amplitude and direction. This provides an important check of the validity of the Fisher analysis.

One may wonder if shifts of less than 1σ are something we should worry about. This means after all that the shifts are hidden in the uncertainty of the measurements. The goal of Euclid is however to achieve an analysis where the sum of all systematic effects is below 0.3σ. In this context, our analysis shows that including magnification in the modelling is necessary (see Fig. 7).

Fig. 7. 68% (inner) and 95% (outer) 2D confidence regions and 1D posteriors for w0waCDM with marginalisation over galaxy biases for the Fisher analysis, with contributions from standard terms (dashed) or standard terms plus lensing magnification (solid). Dashed lack lines denote the fiducial values. For corresponding values of the constraints and shift, refer to Table 4. The label “no magn” refers to an analysis that does not include magnification in the model despite the effect being present in the data. The label “w/magn” incorporates magnification both in the model and in the data.

5.3. Estimation of the growth rate

As is clear from Eq. (14), the quadrupole and hexadecapole of the correlation function are most sensitive to the growth factor f(z). Assuming that μ₄(d, z_*) is given by early Universe parameters determined by CMB observations, one might even use the hexadecapole alone to determine $\stackrel{\sim }{f}(z)$. In practice, however, since the quadrupole and the monopole are much larger and correspondingly measured with much better precision, we use the combined multipoles of the correlations function to estimate both the bias and the growth factor together. In this estimation we assume that the standard cosmological parameters are determined, for example via CMB observations, and we only estimate the unknown functions, $\stackrel{\sim }{b}(z)={\sigma }_8(z)b(z)$ and $\stackrel{\sim }{f}(z)={\sigma }_8(z)f(z)$. Within general relativity we expect f(z)≈Ω_m(z)^0.56 (see Appendix D for the exact expression).

We investigated the effect of lensing magnification on the estimation of the growth rate by fitting the full correlation function including lensing magnification with a model that does not include it. We followed a similar approach as described in Jelic-Cizmek et al. (2021), Breton et al. (2022). In each of the four redshift bins, we varied both the growth rate $\stackrel{\sim }{f}(z_i)={\stackrel{\sim }{f}}_i$ and the bias $\stackrel{\sim }{b}(z_i)={\stackrel{\sim }{b}}_i$. The statistical error in the bias estimated via the Fisher analysis and via an MCMC study is typically of the order of $\sigma ({\stackrel{\sim }{b}}_i)\sim 0.7$% while the error in the growth factor is of the order of $\sigma ({\stackrel{\sim }{f}}_i)\sim 1.2-1.4$% (see Table 5). We note that in this case, adding the magnification in the model would not improve the measurement of ${\stackrel{\sim }{f}}_i$, since the magnification does not depend on these parameters. From Table 5, we see that neglecting lensing magnification in the modelling shifts the best fit values of ${\stackrel{\sim }{b}}_i$ and ${\stackrel{\sim }{f}}_i$ by up to one standard deviation in the highest redshift bin, which is most strongly affected by magnification. But already in bins two and three, neglecting lensing magnification leads to a systematic shift of more than 0.3σ (i.e. above the Euclid target). Comparing the results from the Fisher analysis with those of the MCMC analysis (see Table 5 and Fig. 8), we find excellent agreement, both for the predicted constraints and for the shifts, even in the case where the shifts are larger than 1σ. This is on one hand due to the fact that the posteriors are very close to Gaussian, as can be seen from Fig. 8, and on the other hand, the derivatives of the signal with respect to ${\stackrel{\sim }{b}}_i$ and ${\stackrel{\sim }{f}}_i$ (which are used in the Fisher analysis) are trivial, since these parameters are constant coefficients (in each bin) in front of the scale-dependent functions ${\mu }_{\ell }(d,\overline{z})$ (see Eq. (14)).

Fig. 8. 68% (inner) and 95% (outer) 2D confidence regions and 1D posteriors for the $\{\stackrel{\sim }{f},\stackrel{\sim }{b}\}$ parametrisation, using MCMC (blue) and the Fisher analysis (orange), with only contributions from just standard terms (dashed) and standard terms plus lensing magnification (solid). Dashed black lines denote the fiducial values. For corresponding values of the constraints and shift, refer to Table 5. The label “no magn” refers to an analysis that does not include magnification in the model despite the effect being present in the data. The label “w/magn” incorporates magnification in both the model and the data.

From Fig. 2, we see that the contribution from lensing magnification to all multipoles is positive. This is true at all redshifts, since 5s_i − 2 is positive in all bins (see Table 2). As a consequence, lensing magnification increases the amplitude of the monopole and of the hexadecapole (that is positive) but it reduces the amplitude of the quadrupole (which is negative). Since the constraints come mainly from the monopole and the quadrupole, neglecting lensing magnification in the modelling means therefore that ${\stackrel{\sim }{b}}_i$ and ${\stackrel{\sim }{f}}_i$ are shifted in such a way that ${\stackrel{\sim }{b}}_i^2+2/3{\stackrel{\sim }{b}}_i{\stackrel{\sim }{f}}_i+{\stackrel{\sim }{f}}_i^2/5$ increases, while $4/3{\stackrel{\sim }{b}}_i{\stackrel{\sim }{f}}_i+4/7{\stackrel{\sim }{f}}_i^2$ decreases. This is best achieved by having a negative shift in ${\stackrel{\sim }{f}}_i$ and a positive shift in ${\stackrel{\sim }{b}}_i$. We note that with these shifted values the hexadecapole will not be well fitted, because it would require an increase in ${\stackrel{\sim }{f}}_i$. But since its signal-to-noise ratio is significantly smaller than that of the monopole and the quadrupole, it does not have a significant impact on the analysis.

Such a systematic error in the analysis can certainly not be tolerated. Since the aim of the growth rate analysis is to test the theory of gravity, shifts of more than 1σ in the growth rate would be wrongly interpreted as a detection of modified gravity. However, including lensing magnification in the analysis requires a model, which is exactly what we want to avoid in the growth rate analysis. In the case of the ΛCDM and w₀w_aCDM analyses, the problem is less severe since magnification can be modelled together with density and RSDs. However, including it would significantly enhance the complexity of the computation and slow down the data analysis, especially for parameter estimation using MCMC methods. Below we propose a method for resolving these problems.

5.4. A model for magnification as a cosmology-independent systematic effect

Since lensing magnification is a sub-dominant effect, we can include it in the modelling as a contamination, which does not encode any cosmological information, but that we can model sufficiently well. More precisely, we pre-computed the lensing magnification with fixed cosmological parameters, which were determined, for example, via CMB experiments, and only varied the contributions from density and RSDs in our analysis. We did this in order to remove the bias of cosmological parameters that neglecting magnification can induce. On the other hand, this means that we lose the additional constraining power from lensing; this, however, is not very significant. This “template method”, which uses a fiducial template for the lensing magnification, has also been proposed in Martinelli et al. (2022). Here we tested it both on the ΛCDM and w₀w_aCDM analyses (where it is useful to reduce the computational costs) and on the growth rate analysis. On the growth rate analysis, the template method allows us to preserve the model-independence of the method. As explained before, in the growth rate analysis, the early time cosmology enters at the redshift z_*, before acceleration has started, and it is determined by CMB measurements. The late time evolution is then fully encoded in the parameters ${\stackrel{\sim }{f}}_i$ and ${\stackrel{\sim }{b}}_i$. Any deviations in the laws of gravity would appear as a change in these parameters. Adding the lensing magnification to the signal would spoil the model-independence of the method, since this contribution cannot be easily written in a model-independent way¹³. This would mean that part of the signal is modelled with ${\stackrel{\sim }{f}}_i$ and ${\stackrel{\sim }{b}}_i$, while the other part is modelled in a specific model, for example in ΛCDM. We would then have a mix of parameters, some independent of the theory of gravity, and others specific to ΛCDM. The template method circumvents this problem, by assuming that the lensing magnification is a fixed contribution, independent of cosmological parameters. Of course this is not correct, but we show that the mistake that we make by using this assumption does not introduce any significant shifts in the measurements of the variables ${\stackrel{\sim }{f}}_i$ and ${\stackrel{\sim }{b}}_i$, which we want to constrain with this method. To test this, we included the lensing magnification in the model using wrong cosmological parameters (since in practice we do not know the theory of gravity, nor the value of the true cosmological parameters). More precisely, we used cosmological parameters that are one standard deviation below or above our fiducial values. The ±1σ values are reported in Table 6. We then computed the shifts in the cosmological parameters induced by the fact that the template for the lensing magnification is wrong by ±1 σ.

The shifts in the cosmological parameters for both ΛCDM and w₀w_aCDM are given in Table 7. Comparing with Table 4 we see that the shifts are very significantly reduced with the template method. In w₀w_aCDM, the cosmological parameters never change by more than 0.04σ. In ΛCDM, σ₈ and the galaxy bias are shifted by 0.08σ to 0.13σ when using the template method. We note, however, that we did not vary w₀ and w_a for the template lensing magnification since these parameters are not well determined by CMB data used to obtain the template. The results for the growth rate analysis are presented in Table 8. Again, we see that the template method strongly reduces the shift, which, in the highest redshift bin goes down from 1σ (when magnification is fully neglected) to 0.1σ with the template method. In the other three bins, the shifts are even smaller. We note that lensing magnification, similarly to the lensing effect measured in cosmic shear analysis, is most sensitive to the combination S₈ = σ₈(Ω_m, 0/0.3)^0.5. Since in our test we considered cosmologies ±1σ away from the fiducial, for both σ₈ and Ω_m, 0 at the same time, we show that this method is able to drastically reduce the shift induced by neglecting magnification, even when the model for this effect is significantly inconsistent with the truth.

This template method, where magnification has to be computed only once, is therefore a very promising, inexpensive method for including lensing magnification in the analysis. While in the model-independent analysis of the growth factor magnification has to be modelled as a template, for the full-shape analysis, once the best fit parameters are determined, one will want to include the magnification term with these parameters and run the analysis a second time in order to improve the fit. This iterative method ensures that the result is not sensitive to the initial best fit used to compute the magnification. This is particularly important in light of the current tensions between CMB and large-scale structure constraints. Even though this method has the slight disadvantage that it does not use the information in the lensing magnification to constrain the cosmology, we believe that for a spectroscopic survey, for which magnification is weak, it is the simplest way to avoid the very significant biasing of the results that an analysis that neglects lensing magnification does generate, without significant numerical cost and with a very minor loss of parameter precision (a few percent increase in the error bars).

6. Conclusions

In this paper we have studied the impact of lensing magnification on the spectroscopic survey of Euclid. Lensing magnification is commonly assumed to have a significantly lower impact on the spectroscopic analysis than the photometric one (Euclid Collaboration 2022) due to the fact that: (i) density fluctuations and RSDs have higher amplitudes in a survey with spectroscopic resolution; (ii) the correlation between the different redshift bins are not taken into account in the spectroscopic analysis; and (iii) the multipole expansion used in the spectroscopic analysis removes part of the lensing signal (which is not fully captured by the first three multipoles). Despite this, we find that neglecting magnification leads to significant shifts in the cosmological parameters, by 0.2–0.7 standard deviations. In particular, σ₈, but also Ω_m, 0 and Ω_b, 0, is shifted by more than half a standard deviation.

These shifts become even more significant when we consider the growth rates, which we fitted in an analysis that is independent of the late-time cosmological model. If lensing magnification is neglected, the growth rate is shifted by more than one standard deviation in the highest redshift bin.

From these findings we conclude that it is imperative to include lensing magnification in the data analysis of the Euclid spectroscopic survey. In Appendix E.2 we provide simplified expressions, based on the flat-sky approximation, for the contribution of magnification to the multipoles of the 2PCF (see Jelic-Cizmek 2021, for their derivation). While the cross-correlation of density and magnification can be computed very efficiently, the estimation of the magnification–magnification term is slowed by an integral over the line of sight, which complicates the analysis of the ΛCDM and the w₀w_aCDM models. Even more importantly, the lensing magnification contribution cannot easily be written in a model-independent way. Consequently, including this contribution in the growth rate analysis would spoil the model independence of the method. Since testing the laws of gravity using growth rate analysis is one of the key goals of the spectroscopic analysis, this situation is problematic. Fortunately, we propose a method for solving the problem and reducing the shifts in the parameters to less than 0.1σ while keeping the analysis independent of late-time cosmology. In this so-called “template method”, lensing magnification is calculated in a model with fixed cosmological parameters and simply added to the standard terms. We have shown that if these fixed cosmological parameters deviate by 1σ from the true underlying model, including this slightly wrong contribution from magnification leads to shifts in the inferred cosmological parameters of at most 0.1σ. Of course, one can then improve the analysis by iterating the process.

In our work we have compared our Fisher forecasts with a full MCMC analysis at several stages. We have found that both methods provide consistent results for the parameter shifts, even when they are up to 1σ.

We note that, unlike the previous Euclid forecasts (EP:VII), which are based on the power spectrum in Fourier space, our study presents a configuration-space analysis based on the multipoles of the correlation function. Nevertheless, we expect that our conclusions similarly apply to a Fourier-space analysis. However, modelling the effects of magnification on the multipoles of the power spectrum is less straightforward. First, the contribution of magnification depends on the power spectrum estimator and the survey window function. Second, as demonstrated in Castorina & di Dio (2022), considering the contribution of magnification to the power spectrum multipoles necessitates a prior estimation of the multipoles of the correlation function. Hence, our analysis was conducted directly in configuration space. Even if the covariance matrix has more significant off-diagonal contributions in configuration space, this method has the advantage of being more direct and less survey dependent.

To conclude, we find that including magnification in the analysis does not significantly reduce the error bars on the inferred cosmological parameters. However, we will have to include this in our analysis of the data since otherwise we would fit the data to the wrong physical model, which would bias the inferred cosmological parameters.

Acknowledgments

We acknowledge the use of the HPC cluster Baobab at the University of Geneva for conducting our numerical calculations. G.J.C. acknowledges support from the Swiss National Science Foundation (SNSF), professorship grant (No. 202671). F.S., C.B., R.D. and M.K. acknowledge support from the Swiss National Science Foundation Sinergia Grant CRSII5198674. C.B. acknowledges financial support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (Grant agreement No. 863929; project title “Testing the law of gravity with novel large-scale structure observables”). L.L. is supported by a SNSF professorship grant (No. 202671). P.F. acknowledges support from Ministerio de Ciencia e Innovacion, project PID2019-111317GB-C31, the European Research Executive Agency HORIZON-MSCA-2021-SE-01 Research and Innovation programme under the Marie Skłodowska–Curie grant agreement number 101086388 (LACEGAL). P.F. is also partially supported by the program Unidad de Excelencia María de Maeztu CEX2020-001058-M. C.V. acknowledges an FPI grant from Ministerio de Ciencia e Innovacion, project PID2019-111317GB-C31. This work has made use of CosmoHub. CosmoHub has been developed by the Port d’Informació Científica (PIC), maintained through a collaboration of the Institut de Física d’Altes Energies (IFAE) and the Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas (CIEMAT) and the Institute of Space Sciences (CSIC & IEEC), and was partially funded by the “Plan Estatal de Investigación Científica y Técnica y de Innovación” program of the Spanish government. The Euclid Consortium acknowledges the European Space Agency and a number of agencies and institutes that have supported the development of Euclid, in particular the Academy of Finland, the Agenzia Spaziale Italiana, the Belgian Science Policy, the Canadian Euclid Consortium, the French Centre National d’Études Spatiales, the Deutsches Zentrum für Luft- und Raumfahrt, the Danish Space Research Institute, the Fundação para a Ciência e a Tecnologia, the Ministerio de Ciencia e Innovación, the National Aeronautics and Space Administration, the National Astronomical Observatory of Japan, the Netherlandse Onderzoekschool Voor Astronomie, the Norwegian Space Agency, the Romanian Space Agency, the State Secretariat for Education, Research and Innovation (SERI) at the Swiss Space Office (SSO), and the United Kingdom Space Agency. A complete and detailed list is available on the Euclid web site (http://www.euclid-ec.org).

References

Appendix A: Code validation

The analysis presented in this work is performed using the latest version of the code coffe. In order to make sure that the results of the analysis can be trusted, we first performed a code validation. As a reference, we chose to use the well-established CosmoBolognaLib code from Marulli et al. (2016), which can, among other outputs, compute the redshift-space multipoles of the 2PCF.

The baseline settings used for this code comparison are the same as the ones adopted in EP:VII for the spectroscopic galaxy clustering (GCsp) analysis. In summary:

• The cosmological parameter space is θ = {Ω_m, 0,  Ω_b, 0,  w₀,  w_a,  h,  n_s,  σ₈}, that is, a flat cosmology with dynamical dark energy.

• The galaxy sample is split into four redshift bins, with a galaxy number density as specified in Table 2.

• We included the 4 galaxy bias parameters, one in each redshift bin, as nuisance parameters.

• As we are primarily interested in the validation using linear theory only, we set r_min = 22 Mpc as the smallest separation in each redshift bin.

In Fig. A.1, we present the code comparison. We show the percentage difference between the constraints obtained with the two codes and the mean values of the two results. The top panel refers to 1σ marginalised constraints, while the bottom panel shows the comparison for the un-marginalised constraints. The largest discrepancies between the two codes are ∼2% for the 1σ errors and ∼1% for the un-marginalised constraints. We note that the outcome of the two codes has been compared for several intermediate steps, different settings, and different probe combinations, always leading to an excellent agreement. In particular, we verified that using the covariance from either coffe or CosmoBolognaLib when computing the constraints has no impact on the result; we show a comparison of the two signals for the various redshifts in Fig. A.2. Even though the hexadecapoles show differences up to 10% and larger in the vicinity of the baryon acoustic oscillation peak, due to the small amplitude of this contribution, this is not relevant for parameter estimation. In Fig. A.3 we show the monopole from both coffe and CosmoBolognaLib at z = 1, with a close-up of some points of interest, notably, the baryon acoustic oscillation peak at ∼150 Mpc and the zero crossing at ∼180 Mpc.

Fig. A.1. Percentage difference between coffe and CosmoBolognaLib in terms of the 1σ uncertainties (top panel) and un-marginalised constraints (bottom panel) for the spectroscopic sample of galaxy clustering. This analysis includes four nuisance parameters for the galaxy bias that are marginalised over in the 1σ constraints.

Fig. A.2. Percentage difference between coffe and CosmoBolognaLib in the first three even multipoles of the 2PCF, for various redshifts. The large ‘jump’ in the monopole around r ∼ 180 Mpc is caused by its passage through zero. The black dashed lines denote a 1% threshold.

Fig. A.3. Comparison between the monopole from coffe and CosmoBolognaLib at $\overline{z}=1$. The dashed line denotes the zero crossing.

In order to obtain the Fisher matrix, we need to compute derivatives of the multipoles of the 2PCF with respect to cosmological parameters. As it is not possible to compute them analytically, we resorted to the method of finite differences. Since this method suffers from numerical instabilities, it is necessary to first find the optimal step size that is neither too large (causing the derivative to be too ‘coarse’) nor too small (resulting in errors due to numerical underflow). The final step size used (10⁻³ for all parameters) has proven to be sufficiently accurate at the level of the obtained marginalised constraints for parameters of the w₀w_aCDM model. As shown in Fig. A.4, the difference of the inferred parameters for step sizes 10⁻³ and 10⁻⁴ is always below 2% of the standard deviation.

Fig. A.4. Percentage difference between step sizes 10−3 and 10−4 in the 1σ uncertainties (top panel) and un-marginalised constraints (bottom panel) for the spectroscopic sample of galaxy clustering.

Appendix B: Fitting functions for b(z) and s(z)

For convenience, we provide a polynomial fit (obtained using the Levenberg-Marquardt algorithm) for both the galaxy bias and the local count slope as given below. We set

$$\begin{array}{cc}& b(z)={\stackrel{\sim }{b}}_0+{\stackrel{\sim }{b}}_{1}z+{\stackrel{\sim }{b}}_2{z}^2+{\stackrel{\sim }{b}}_3{z}^3,\hfill \end{array}$$(B.1)

$$\begin{array}{cc}& s(z)={\stackrel{\sim }{s}}_0+{\stackrel{\sim }{s}}_{1}z+{\stackrel{\sim }{s}}_2{z}^2+{\stackrel{\sim }{s}}_3{z}^3,\hfill \end{array}$$(B.2)

with parameters

$$\begin{array}{c}\hfill \begin{array}{cccc}{\stackrel{\sim }{b}}_0=0.853,& {\stackrel{\sim }{b}}_1=0.040,& {\stackrel{\sim }{b}}_2=0.713,& {\stackrel{\sim }{b}}_3=-0.164,\\ {\stackrel{\sim }{s}}_0=1.231,& {\stackrel{\sim }{s}}_1=-1.746,& {\stackrel{\sim }{s}}_2=1.810,& {\stackrel{\sim }{s}}_3=-0.505.\end{array}\end{array}$$(B.3)

In Fig. 3 we compare our best fit with the measurements of the Flagship simulation. In our calculations we do not use these fits, but we present them here for convenience. The Flagship specifics have been estimated for the survey binning as described in Sect. 3, and therefore the fitting functions are adapted to this specific configuration.

Appendix C: Binned covariance of 2PCF multipoles

In this section we report the exact expression for the covariance adopted in this work. We modified the implementation of the flat-sky, Gaussian covariance reported in Tansella et al. (2018b) in order to include a binning of the Bessel functions. The full expression of the covariance can be written in the same form as in Tansella et al. (2018b):

$$\begin{array}{cc}\hfill \mathsf{C}[{\xi }_{\ell }^j,{\xi }_{{\ell }^{\prime }}^k]& =\frac{i^{\ell -{\ell }^{\prime }}}{V}[\frac{1}{\overline{\mathcal{N}}}{\overline{\mathcal{G}}}_{\ell {\ell }^{\prime }}(d_j,d_k,\overline{z}){\displaystyle \sum _{\sigma }}{c}_{\sigma }{\left(\begin{array}{ccc}\ell & {\ell }^{\prime }& \sigma \\ 0& 0& 0\end{array}\right)}^2\hfill \\ \hfill & +{\overline{\mathcal{D}}}_{\ell {\ell }^{\prime }}(d_j,d_k,\overline{z}){\displaystyle \sum _{\sigma }}{\stackrel{\sim }{c}}_{\sigma }{\left(\begin{array}{ccc}\ell & {\ell }^{\prime }& \sigma \\ 0& 0& 0\end{array}\right)}^2\hfill \\ \hfill & +\frac{2\ell +1}{2\pi {\overline{\mathcal{N}}}^2{d}_j^2{L}_p}{\delta }_{\mathit{jk}}{\delta }_{\ell {\ell }^{\prime }}],\hfill \end{array}$$(C.1)

where the $\left(\begin{array}{c}\ell {\ell }^{\prime }\sigma \\ 000\end{array}\right)$ denote Wigner 3j symbols; hence, the sum over σ goes from |ℓ−ℓ′| to ℓ + ℓ′. The three terms in the sum respectively denote the cross-correlation between cosmic variance and Poisson noise, the cosmic variance autocorrelation, and the Poisson noise autocorrelation. Here V denotes the comoving volume of the observed sample, 𝒩 denotes the average comoving number density of sources, and L_p denotes the pixel size, that is, the minimum comoving distance we can resolve. We note that the coefficients $\{c_{\sigma },{\stackrel{\sim }{c}}_{\sigma }\}$ depend only on redshift. Their exact expressions are reported in Tansella et al. (2018b), and we repeat them here for the sake of completeness:

$$\begin{array}{cc}& c_0=b^2+\frac{2}{3}bf+\frac{f^2}{5},\hfill \end{array}$$(C.2)

$$\begin{array}{cc}& c_2=\frac{4}{3}bf+\frac{4}{7}{f}^2,\hfill \end{array}$$(C.3)

$$\begin{array}{cc}& c_4=\frac{8}{35}{f}^2,\hfill \end{array}$$(C.4)

$$\begin{array}{cc}& {\stackrel{\sim }{c}}_0=c_0^2+\frac{c_2^2}{5}+\frac{c_4^2}{9},\hfill \end{array}$$(C.5)

$$\begin{array}{cc}& {\stackrel{\sim }{c}}_2=\frac{2}{7}{c}_2(7c_0+c_2)+\frac{4}{7}{c}_2{c}_4+\frac{100}{693}{c}_4^2,\hfill \end{array}$$(C.6)

$$\begin{array}{cc}& {\stackrel{\sim }{c}}_4=\frac{18}{35}{c}_2^2+2c_0{c}_4+\frac{40}{77}{c}_2{c}_4+\frac{162}{1001}{c}_4^2,\hfill \end{array}$$(C.7)

$$\begin{array}{cc}& {\stackrel{\sim }{c}}_6=\frac{10}{99}{c}_4(9c_2+2c_4),\hfill \end{array}$$(C.8)

$$\begin{array}{cc}& {\stackrel{\sim }{c}}_8=\frac{490}{1287}{c}_4^2.\hfill \end{array}$$(C.9)

The main difference between the original implementation in coffe and the covariance used in this analysis lies in the computation of ${\overline{\mathcal{G}}}_{\ell {\ell }^{\prime }}$ and ${\overline{\mathcal{D}}}_{\ell {\ell }^{\prime }}$. Here these are estimated as integrals of the binned spherical Bessel functions:

$$\begin{array}{c}\hfill {\overline{j}}_{\ell }(k{d}_i):=\frac{4\pi }{V_{d_i}}{\displaystyle {\int }_{d_i-L_p/2}^{d_i+L_p/2}\mathrm{d}s{s}^2{j}_{\ell }(ks),}\end{array}$$(C.10)

where $V_{d_i}=\frac{4\pi }{3}(d_{i,\max }^3-d_{i,\min }^3)$ is the volume of the distance bin around d_i. Thus, we have

$$\begin{array}{cc}& \frac{{\overline{\mathcal{G}}}_{\ell {\ell }^{\prime }}(d_i,d_j,\overline{z})}{(2\ell +1)(2{\ell }^{\prime }+1)}=\frac{2}{{\pi }^2}{\displaystyle {\int }_0^{\infty }\mathrm{d}k{k}^{2}P(k,\overline{z}){\overline{j}}_{\ell }(k{d}_i){\overline{j}}_{{\ell }^{\prime }}(k{d}_j),}\hfill \\ \hfill & \frac{{\overline{\mathcal{D}}}_{\ell {\ell }^{\prime }}(d_i,d_j,\overline{z})}{(2\ell +1)(2{\ell }^{\prime }+1)}=\frac{1}{{\pi }^2}{\int }_0^{\infty }\mathrm{d}k{k}^2{P}^2(k,\overline{z}){\overline{j}}_{\ell }(k{d}_i){\overline{j}}_{{\ell }^{\prime }}(k{d}_j).\hfill \end{array}$$

In fact, it is shown in the literature that the covariance matrix is overestimated when this volume average over the spherical Bessel functions is not applied (see Grieb et al. 2016).

In Fig. C.1 we show the ratio of the diagonal entries of the unbinned and the binned covariance; as we can see, for low separations and large multipoles, the unbinned covariance can be larger than its binned counterpart by more than 20%.

Fig. C.1. Diagonal entries of the ratio of the unbinned and the binned covariance, for various multipoles, as a function of comoving separation, d, at the lowest redshift bin of Euclid.

Appendix D: An analytic expression for the growth rate

In Sect. 2.2 we introduced the growth rate $f(z):=\frac{\mathrm{d}ln\delta }{\mathrm{d}lna}$. Here we derive an analytical expression for f(z) in linear perturbation theory and compare it to the commonly used expression f(z)≈Ω_m(z)^0.56. The linear density fluctuation in comoving gauge is δ(z)≈D₁(z) δ(0). Following Durrer (2020), we obtain

$$\begin{array}{c}\hfill f(z)=\frac{\mathrm{d}ln{D}_1}{\mathrm{d}lna}=-(1+z)\frac{\mathrm{d}ln{D}_1}{\mathrm{d}z},\end{array}$$(D.1)

where D₁ is the linear growth factor. In a ΛCDM universe, D₁ is the growing mode solution of the following equation:

$$\begin{array}{c}\hfill \ddot{D_1}+\mathcal{H}\dot{D_1}=\frac{3}{2}{\mathcal{H}}^2{\Omega }_{\mathrm{m}}(a)D_1,\end{array}$$(D.2)

with

$$\begin{array}{c}\hfill {\Omega }_{\mathrm{m}}(a)=\frac{{\Omega }_{\mathrm{m},0}{a}^{-3}}{{\Omega }_{\mathrm{m},0}{a}^{-3}+(1-{\Omega }_{\mathrm{m},0})}.\end{array}$$(D.3)

We rewrite Eq. (D.2) using a prime to indicate the derivative with respect to lna:

$$\begin{array}{c}\hfill {D^{″}}_1+(2-\frac{3}{2}{\Omega }_{\mathrm{m}}(a)){D^{\prime }}_1=\frac{3}{2}{\Omega }_{\mathrm{m}}(a)D_1.\end{array}$$(D.4)

The analytical solution of Eq. (D.4) in terms of hypergeometric functions (see Abramowitz & Stegun 1972), leads to the following expression for the growth rate:

$$\begin{array}{c}\hfill f(z)=\frac{1}{2}{\Omega }_{\mathrm{m}}(z)[\frac{5}{{}_2{F}_1(\frac{1}{3},1;\frac{11}{6};1-\frac{1}{{\Omega }_{\mathrm{m}}(z)})}-3].\end{array}$$(D.5)

In the literature (see e.g. Durrer 2020; Linder & Cahn 2007), one often finds the following approximation for the growth rate:

$$\begin{array}{c}\hfill f(z)\approx {\Omega }_{\mathrm{m}}^{0.56}(z).\end{array}$$(D.6)

A comparison between eqs. Eq. (D.5) and Eq. (D.6) is shown in Fig. D.1. The approximation is clearly excellent, leading to differences below 1% at all redshifts.

Fig. D.1. An analytic expression for the growth rate. Top: Exact expression for the growth rate, given by Eq. (D.5) (solid blue curve), and the approximate expression, given by Eq. (D.6) (dashed orange curve). Bottom: Relative difference between them, in percentage points. In both cases we assume Ωm, 0 = 0.32.

Appendix E: Expressions for the multipoles of the 2PCF

E.1. Curved-sky expressions

For completeness, below we provide the relevant curved-sky contributions for the 2PCF, which were first derived in Tansella et al. (2018a). We note that, for the sake of brevity, in the following equations we use the coordinates r₁ := r(z₁), r₂ := r(z₂), and θ, which denote the comoving distances at redshifts z₁ and z₂ and the angle at the observer, respectively. Furthermore, we use x_i := x(z_i).

We used the following notation for local terms A and B,

$$\begin{array}{cc}\hfill {\xi }^{\mathit{AB}}(d,\overline{z},\mu )& =D_1(z_1)D_1(z_2)\hfill \\ \hfill & \times {\displaystyle \sum _{\ell ,n}}(X_{\ell }^n{|}_A+X_{\ell }^n{|}_{\mathit{AB}}+X_{\ell }^n{|}_{\mathit{BA}}+X_{\ell }^n{|}_B)I_{\ell }^n(d),\hfill \end{array}$$(E.1)

and the following for the integrated terms,

$$\begin{array}{cc}& {\xi }^{\mathit{AB}}(d,\overline{z},\mu )=Z{|}_A+Z{|}_{\mathit{AB}}+Z{|}_{\mathit{BA}}+Z{|}_B.\hfill \end{array}$$(E.2)

$$\begin{array}{cc}& X_0^0{|}_{\mathrm{den}}=b_1{b}_2,\hfill \end{array}$$(E.3)

$$\begin{array}{cc}& X_0^0{|}_{\mathrm{RSD}}=f_1{f}_2\frac{1+2{cos}^2\theta }{15},\hfill \end{array}$$(E.4)

$$\begin{array}{cc}& X_2^0{|}_{\mathrm{RSD}}=-\frac{f_1{f}_2}{21}[1+11{cos}^2\theta +\frac{18cos\theta ({cos}^2\theta -1)r_1{r}_2}{d^2}],\hfill \end{array}$$(E.5)

$$\begin{array}{cc}& X_4^0{|}_{\mathrm{RSD}}=\frac{f_1{f}_2}{35r^4}\{4(3{cos}^2\theta -1)({r_1}^4+{r_2}^4)\hfill \\ \hfill & +r_1{r}_2(3+{cos}^2\theta )[3(3+{cos}^2\theta )r_1{r}_2\hfill \end{array}$$(E.6)

$$\begin{array}{cc}& -8({r_1}^2+{r_2}^2)cos\theta ]\},\hfill \\ \hfill & X_0^0{|}_{\mathrm{den}-\mathrm{RSD}}=\frac{b_1{f}_2}{3},\hfill \end{array}$$(E.7)