© The Authors 2025

1. Introduction

The primary science goal of the recently launched Euclid space mission (Euclid Collaboration: Mellier et al. 2025) is to test whether a cosmological constant can be ruled out as the driver of the accelerated expansion of the Universe. To that end, Euclid is carrying out a wide-angle survey covering nearly 15 000 deg² of the extragalactic sky. The Euclid mission is optimised for the combination of two cosmological probes – weak gravitational lensing and galaxy clustering – and relies on two instruments. The visual imager (VIS, Euclid Collaboration: Borlaff et al. 2022) operates in the 550 to 900 nm pass-band and produces high-quality galaxy images to perform measurements of galaxy shapes (cosmic shear). The near-infrared spectrometer and photometer (NISP, Maciaszek et al. 2022; Euclid Collaboration: Schirmer et al. 2022) carries out imaging photometry (as an input for the estimation of photometric redshifts) and slitless spectroscopy to precisely measure the redshift of the Hα emission line in the range 0.9 < z < 1.8. Over six years of observations, the Euclid wide spectroscopic survey (EWSS) will measure the redshifts and positions of nearly 30 million emission-line galaxies, while its photometric counterpart will measure the positions and shapes of approximately 1.5 billion galaxies.

Galaxy clustering (the only cosmological probe discussed in this paper) sets constraints on the cause of the accelerated expansion of the Universe in two ways. First, the expansion history of the Universe can be reconstructed by locating the characteristic scale imprinted by baryonic acoustic oscillations on the galaxy power spectrum (or the two-point correlation function, 2PCF) as a function of redshift (e.g. Cole et al. 2005; Eisenstein et al. 2005). Second, the growth rate of structure can be determined by studying the anisotropy of the clustering signal (e.g. Peacock et al. 2001). The latter option is generally known as the study of redshift-space distortions (RSDs), which arise when galaxy redshifts are mapped into distances by assuming an unperturbed Friedmann-Lemaître-Robertson-Walker (FLRW) metric. The observed redshift of a galaxy does not coincide with its cosmological component. The dominant (but not sole) correction is due to the relative peculiar velocity between the galaxy and the observer along the line of sight. In a seminal work, Kaiser (1987) used linear perturbation theory to calculate how peculiar velocities distort the galaxy power spectrum (in the Cartesian Fourier basis, see also Hamilton & Culhane 1996; Hamilton 1998, 2000, for extensions to configuration space). The derivation relies on a simplifying assumption that the size of the surveyed region is negligibly small compared to its distance from the observer, so that the lines of sight to all galaxies are effectively parallel. This is nowadays known as the ‘global plane-parallel’ (GPP) approximation and implies that the galaxy power spectrum depends on the cosine of the angle between the wavevector and the fixed line of sight. When decomposed in Legendre polynomials, this functional dependence only includes multipoles of degree 0, 2, and 4.

The EWSS provides us with the opportunity to study galaxy clustering on unprecedentedly large scales. This possibility, however, brings forth new challenges. First of all, galaxy pairs with large angular separations contribute to the clustering signal and might cause systematic deviations of the observations from models based on the GPP approximation. Wide-angle effects have been investigated both for the galaxy 2PCF (see e.g. Matsubara 2000a; Szapudi 2004; Pápai & Szapudi 2008; Raccanelli et al. 2010; Samushia et al. 2012; Raccanelli et al. 2018) and the power spectrum (see e.g. Zaroubi & Hoffman 1996; Reimberg et al. 2016; Castorina & White 2018; Castorina et al. 2019). The current leading opinion is that accounting for a variable line-of-sight using the ‘local plane-parallel’ (LPP) approach (e.g. Beutler et al. 2014; Wilson et al. 2017) is sufficient for conducting cosmological studies with surveys of similar size and depth to the EWSS (Castorina & Di Dio 2022; Noorikuhani & Scoccimarro 2023). An alternative approach that fully accounts for wide-angle effects is to forgo the plane-wave expansion of the Fourier-space overdensity and instead adopt the ‘spherical Fourier-Bessel’ formalism (e.g. Heavens & Taylor 1995; Hamilton 2005; Yoo & Desjacques 2013). The second complication arises from the fact that Kaiser’s RSDs should be complemented with additional corrections. The light bundles from distant galaxies to us propagate through the inhomogeneous Universe and are thus subject to effects like gravitational lensing or aberration. Hence, the observed galaxy positions on the sky, redshifts, and fluxes differ from their analogues obtained in the corresponding unperturbed FLRW model. However, we construct maps of the galaxy distribution by assuming such a homogeneous model in order to convert the observed properties into three-dimensional positions and luminosities. This step introduces a number of artefacts in the galaxy overdensity field (Yoo et al. 2009; Bonvin & Durrer 2011; Challinor & Lewis 2011; Jeong et al. 2012) that we refer to as ‘relativistic¹ RSDs’ and that are also known in the literature as ‘relativistic effects’ or ‘projection effects’ (since we observe the projection of the actual Universe on our past light cone). These alterations can be studied perturbatively. The leading term coincides with Kaiser’s RSDs due to the peculiar velocity gradient along the line of sight. Nevertheless, there exist several additional corrections that can potentially influence galaxy-clustering statistics on very large scales. A number of investigations have characterised the impact of these terms on the angular power spectrum (Di Dio et al. 2013; Borzyszkowski et al. 2017), the 2PCF (e.g. Bertacca 2015; Raccanelli et al. 2016; Tansella et al. 2018; Bertacca 2020; Jelic-Cizmek et al. 2021; Breton et al. 2022, 2019), and the (3D) power spectrum (Elkhashab et al. 2021; Castorina & Di Dio 2022; Foglieni et al. 2023; Noorikuhani & Scoccimarro 2023).

Finding out which of the relativistic RSDs will leave a detectable imprint on the two-point summary statistics measured from the EWSS is the main goal of this study². We investigated four different suites of mock galaxy catalogues built with the LIGER (LIght cones with GEneral Relativity) method (Borzyszkowski et al. 2017; Elkhashab et al. 2021), which allows us to self-consistently correct the output of Newtonian N-body simulations and introduce relativistic RSDs to linear order in the cosmological perturbations. All the mock catalogues we generated match the survey geometry and selection function of the EWSS but the four kinds we consider differ in the number of RSD terms they include. This helped us to isolate the contributions from various effects (e.g. gravitational lensing or the peculiar velocity of the observer). We made use of popular estimators to measure clustering summary statistics from the mock catalogues and we built unbiased models that exactly account for wide-angle effects by averaging the measurements over a large number of realisations. For the power spectrum and 2PCF multipoles, we adopted the standard statistics (and estimators) that will be used for the cosmological analysis for the EWSS. Finally, we employed the likelihood-ratio test to quantify the signal-to-noise ratio (S/N) with which certain effects can be detected and to determine the fraction of realisations in which models that do not account for these effects could be ruled out with a given statistical significance.

The paper is structured as follows. We introduce the relativistic RSDs in Sect. 2 and our mock EWSS galaxy catalogues in Sect. 3. We describe how we assess the detectability of various relativistic RSDs in Sect. 4. Our results for the angular power spectrum, the multipole moments of the 2PCF, and the multipole moments of the power spectrum are presented in Sects. 5, 6, and 7, respectively. Moreover, in Sect. 8, we compare the power-spectrum multipoles extracted from the mocks to the predictions from Kaiser’s model after accounting for the window function of the survey and the integral constraint. Eventually, in Sect. 9, we summarise our findings and conclude.

Throughout this paper, we adopt Einstein’s summation convention and define the space-time metric tensor to have the signature ( − , + , + , + ). Greek indices refer to space-time components (i.e. run from 0 to 3), while Latin indices label spatial components (i.e. run from 1 to 3). Furthermore, the Dirac delta and the Kronecker delta functions are denoted by the symbols δ_D and δ_K, respectively. Our Fourier-transform convention is $\stackrel{\sim }{f}(\mathit{k})={\displaystyle \int f}(\mathit{x}){\mathrm{e}}^{-i\mathit{k}·\mathit{x}}{\mathrm{d}}^{3}x$. Finally, the symbol c denotes the speed of light in vacuum.

2. Relativistic RSDs

In order to build three-dimensional maps of the galaxy distribution, it is usually assumed that the light bundles emitted by the galaxies propagate in an unperturbed FLRW model universe and that their observed redshift, z_obs, coincides with the cosmological one. This implies that their comoving distance in the so-called ‘redshift space’ is

$$\begin{array}{c}\hfill x={\displaystyle {\int }_0^{z_{\mathrm{obs}}}}\frac{c}{H(z)}\mathrm{d}z,\end{array}$$(1)

where H(z) denotes the Hubble parameter in the model universe as a function of redshift. This procedure, however, neglects the fact that inhomogeneities in the Universe alter the observed redshifts and angular positions of the galaxies. Therefore, the reconstructed galaxy maps in redshift space are not faithful (Sargent & Turner 1977). A number of effects, collectively called redshift-space distortions, artificially shift the reconstructed positions of galaxies in both the radial and tangential directions with respect to their actual (hereafter, real-space) location.

The pioneering work by Kaiser (1987) investigated the relationship between galaxy densities in real and redshift space at linear order in the cosmological perturbations, focusing on the impact of peculiar velocities generated by gravitational instabilities (see also Hamilton & Culhane 1996; Hamilton 1998; Matsubara 2000b). More recently, this subject has been revisited using a fully general-relativistic approach and accounting for additional effects like gravitational lensing, the Sachs–Wolfe effects, and the Shapiro delay (Yoo et al. 2009; Bonvin & Durrer 2011; Challinor & Lewis 2011; Jeong et al. 2012). In this latter case, the goal is to compute the geodesics of photons emitted from a source galaxy in the presence of linear cosmological perturbations. This is sufficient to address the large spatial scales considered in this paper. In the remainder of this section we summarise the main results obtained within the general-relativistic framework.

By definition, the co-ordinates of a distant galaxy in redshift space can be trivially expressed as

$$\begin{array}{c}\hfill x^{\mu }=(c{\eta }_0-x,x\mathit{n}),\end{array}$$(2)

where η₀ is the present-day value of conformal time (i.e. at observation), x denotes the comoving distance from the observer (see Eq. 1), and n is the observed galaxy position (pointing towards the galaxy) on the sky. The mapping between the real- and redshift-space co-ordinates of a galaxy can be generically written as x_r^μ = x^μ + Δx^μ. In order to compute the co-ordinate transformation explicitly, we need to specify a gauge. We express the space-time metric in the Poisson gauge, assuming a flat cosmology while neglecting vector and tensor perturbations³,

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

where Ψ and Φ are the dimensionless Bardeen potentials, η is the conformal time, and a is the cosmic scale factor. From this choice, it follows that (Hui & Greene 2006; Yoo et al. 2009; Bonvin & Durrer 2011; Challinor & Lewis 2011; Jeong et al. 2012)

$$\begin{array}{cc}\hfill \mathrm{\Delta }{x}^0& =\frac{c}{\mathcal{H}}\delta lna,\hfill \end{array}$$(4a)

$$\begin{array}{cc}\hfill \mathrm{\Delta }{x}^i& =-({\mathrm{\Phi }}_{\mathrm{o}}+{\mathrm{\Psi }}_{\mathrm{o}}+\frac{{\mathit{v}}_{\mathrm{e}}}{c}·\mathit{n})x^i-x\frac{v_{\mathrm{o}}^i}{c}-\frac{c}{\mathcal{H}}{n}^i\delta lna\hfill \\ \hfill & +n^i{\displaystyle {\int }_0^x}(x-\stackrel{\sim }{x})\frac{(\mathrm{\Phi }\prime +\mathrm{\Psi }\prime )}{c}\mathrm{d}\stackrel{\sim }{x}-{\displaystyle {\int }_0^x}(x-\stackrel{\sim }{x}){\delta }_{\mathrm{K}}^{\mathit{ij}}{\stackrel{\sim }{\partial }}_j(\mathrm{\Phi }+\mathrm{\Psi })\mathrm{d}\stackrel{\sim }{x}\hfill \\ \hfill & +2n^i{\displaystyle {\int }_0^x}(\mathrm{\Phi }+\mathrm{\Psi })\mathrm{d}\stackrel{\sim }{x},\hfill \end{array}$$(4b)

where

$$\begin{array}{c}\hfill \delta lna:=[\frac{({\mathit{v}}_{\mathrm{e}}-{\mathit{v}}_{\mathrm{o}})}{c}·\mathit{n}-({\mathrm{\Phi }}_{\mathrm{e}}-{\mathrm{\Phi }}_{\mathrm{o}})-{\displaystyle {\int }_0^x}\frac{(\mathrm{\Phi }\prime +\mathrm{\Psi }\prime )}{c}\mathrm{d}\stackrel{\sim }{x}],\end{array}$$(5)

represents the fractional redshift change due to the perturbations; that is, −δz/(1 + z). Here, v denotes peculiar three-velocities, the subscripts ‘e’ and ‘o’ specify whether the functions are evaluated at the source or observer locations, respectively, ${\stackrel{\sim }{\partial }}_j=\partial /\partial {\stackrel{\sim }{x}}^j$, ℋ = a(z_obs) H(z_obs), and the prime superscript denotes the partial derivative w.r.t. η. It is worth mentioning that the equations above assume that the peculiar velocity of a galaxy coincides with that of the matter at the same location; in other words, that there is no velocity bias⁴.

Cosmological perturbations also alter the solid angle under which galaxies are seen by distant observers, thus enhancing, or decreasing their apparent flux (e.g. Broadhurst et al. 1995). In terms of the luminosity distance, D_L, the magnification of a galaxy is defined as

$$\begin{array}{c}\hfill \mathcal{M}={\left(\frac{D_{\mathrm{L}}}{{\overline{D}}_{\mathrm{L}}}\right)}^{-2},\end{array}$$(6)

where ${\overline{D}}_{\mathrm{L}}$ denotes the luminosity distance in the background model universe evaluated at z_obs. At linear order (e.g. Challinor & Lewis 2011; Bertacca 2015),

$$\begin{array}{cc}\hfill \mathcal{M}& =1+2{\mathrm{\Phi }}_{\mathrm{e}}-2(1-\frac{c}{\mathcal{H}x})\delta lna-2\frac{{\mathit{v}}_{\mathrm{o}}}{c}·\mathit{n}\hfill \\ \hfill & +2\kappa -\frac{2}{x}{\displaystyle {\int }_0^x}(\mathrm{\Phi }+\mathrm{\Psi })\mathrm{d}\stackrel{\sim }{x},\hfill \end{array}$$(7)

where the weak lensing convergence is

$$\begin{array}{c}\hfill \kappa :=\frac{1}{2}{\displaystyle {\int }_0^x}(x-\stackrel{\sim }{x})\frac{\stackrel{\sim }{x}}{x}[{\mathrm{\nabla }}^2+{(\mathit{n}·\mathbf{\nabla })}^2-\frac{2}{x}\mathit{n}·\mathbf{\nabla }](\mathrm{\Phi }+\mathrm{\Psi })\mathrm{d}\stackrel{\sim }{x}.\end{array}$$(8)

The next step is to understand how the local number density of galaxies in a survey responds to redshift perturbations and magnification. To first approximation, the EWSS is flux limited as it only selects galaxies above a given observed Hα flux⁵, corresponding to a redshift-dependent luminosity limit, L_lim(z). We indicate by n(L_min, z) the mean number density of the target population of galaxies with luminosity L > L_min at redshift z. Then, the so-called evolution bias,

$$\begin{array}{c}\hfill \mathcal{E}(z)=-{\frac{\partial lnn(L_{\min },z)}{\partial ln(1+z)}|}_{L_{\min }=L_{\lim }(z)},\end{array}$$(9)

quantifies how rapidly the number density of the selected galaxies changes with redshift. Similarly, the magnification bias,

$$\begin{array}{c}\hfill \mathcal{Q}(z)=-{\frac{\partial lnn(L_{\min },z)}{\partial ln{L}_{\min }}|}_{L_{\min }=L_{\lim }(z)},\end{array}$$(10)

gives the slope of the cumulative luminosity function evaluated at the luminosity limit of the survey. Taking into account all linear-order corrections to the ‘observed’ galaxy density (in redshift space), n_g, and its angular average at fixed z_obs, ${\overline{n}}_{\mathrm{g}}$, it is possible to express the overdensity, ${\delta }_{\mathrm{g},\mathrm{s}}:=n_{\mathrm{g}}/{\overline{n}}_{\mathrm{g}}-1$, in terms of the cosmological perturbations (Yoo et al. 2009; Challinor & Lewis 2011; Jeong et al. 2012; Bertacca 2015)

$$\begin{array}{cc}\hfill {\delta }_{\mathrm{g},\mathrm{s}}(\mathit{x})& ={\delta }_{\mathrm{g}}^{\mathrm{com}}-\frac{1}{\mathcal{H}}\frac{\partial ({\mathit{v}}_{\mathrm{e}}·\mathit{n})}{\partial x}+2(\mathcal{Q}-1)\kappa \hfill \\ \hfill & +[\mathcal{E}-2\mathcal{Q}-\frac{\mathcal{H}\prime }{{\mathcal{H}}^2}-\frac{2(1-\mathcal{Q})c}{x\mathcal{H}}]\hfill \\ \hfill & \times [\frac{{\mathit{v}}_{\mathrm{e}}}{c}·\mathit{n}-({\mathrm{\Phi }}_{\mathrm{e}}-{\mathrm{\Phi }}_{\mathrm{o}})-{\displaystyle {\int }_0^x}\frac{(\mathrm{\Phi }\prime +\mathrm{\Psi }\prime )}{c}\mathrm{d}\stackrel{\sim }{x}]\hfill \\ \hfill & -2(1-\mathcal{Q}){\mathrm{\Phi }}_{\mathrm{e}}+{\mathrm{\Psi }}_{\mathrm{e}}+\frac{{\mathrm{\Phi }}_{\mathrm{e}}\prime }{\mathcal{H}}+(3-\mathcal{E})\frac{\mathcal{H}\mathrm{\Theta }}{c^2}\hfill \\ \hfill & +\frac{2(1-\mathcal{Q})}{x}{\displaystyle {\int }_0^x}(\mathrm{\Phi }+\mathrm{\Psi })\mathrm{d}\stackrel{\sim }{x}\hfill \\ \hfill & +[2-\mathcal{E}+\frac{\mathcal{H}\prime }{{\mathcal{H}}^2}+\frac{2(1-\mathcal{Q})c}{x\mathcal{H}}]\frac{{\mathit{v}}_{\mathrm{o}}}{c}·\mathit{n},\hfill \end{array}$$(11)

where Θ is the linear velocity potential (i.e. v = ∇ Θ). We note that the real-space galaxy overdensity, ${\delta }_{\mathrm{g}}^{\mathrm{com}}$, is defined in the synchronous comoving gauge while all the rest is set in the Poisson gauge⁶. At linear order, ${\delta }_{\mathrm{g}}^{\mathrm{com}}$ is related to the underlying matter density fluctuation through the linear bias parameter b; that is, ${\delta }_{\mathrm{g}}^{\mathrm{com}}=b{\delta }_{\mathrm{m}}^{\mathrm{com}}$ (Challinor & Lewis 2011; Jeong et al. 2012).

Equation (11) defines what is meant by ‘relativistic (linear) RSDs’ and forms the starting point for our study. In brief, it says that the galaxy overdensities in real and redshift space differ because of a number of physical effects. The second term on the rhs is the classic Kaiser correction due to the variation in peculiar velocities along the line of sight (Kaiser 1987). The third term is the weak lensing contribution due to volume and magnification corrections, which is expected to have an effect on different clustering statistics on large scales (e.g. Matsubara 2000b; Hui et al. 2007, 2008; Yoo et al. 2009; Challinor & Lewis 2011; Camera et al. 2015; Raccanelli et al. 2016; Borzyszkowski et al. 2017). There are then several additional corrections that depend on the gravitational potentials and the peculiar velocities. Previous studies have shown that those proportional to the peculiar velocity of the observer – that is, the last term in Eq. (11) – could generate observable features in the 2PCF at wide angles (Bertacca 2020) as well as superimpose an oscillatory signal to the power-spectrum monopole at very large scales (Elkhashab et al. 2021, see also Bahr-Kalus et al. 2021), dubbed the finger-of-the-observer effect.

3. Galaxy mocks

In this paper, we use a suite of realistic mock galaxy catalogues to determine which of the corrections appearing in Eq. (11) should be accounted for in the analysis of two-point clustering statistics extracted from the EWSS. The main steps to generate the mock catalogues are described below (for further details, see Elkhashab et al. 2021).

3.1. The LIGER method

LIGER (Borzyszkowski et al. 2017; Elkhashab et al. 2021) is a numerical tool for building mock realisations of the galaxy distribution on the past light cone of an observer⁷. As an input, it takes a Newtonian cosmological simulation. This can be either a hydrodynamic simulation or an N-body simulation combined with a semi-analytic model of galaxy formation to which one applies the selection criteria of a given survey. However, resolving individual galaxies within very large comoving volumes is extremely challenging and time consuming with current software and facilities. Therefore, a special ‘large-box’ mode has been developed in which a Newtonian N-body simulation is used in combination with a set of functions describing the galaxy population under study (i.e. their number density, linear bias parameter, magnification bias, and evolution bias). The key argument underlying this approach is that, at linear order and for a pressureless fluid in a ΛCDM background, Ψ = Φ in the Poisson gauge⁸. The potentials satisfy the standard Poisson equation and can be computed starting from the matter overdensity in the Newtonian simulations, which is equivalent to its counterpart in the synchronous comoving gauge, ${\delta }_{\mathrm{m}}^{\mathrm{com}}$ (for more details, see Sect. 2.1.4 in Borzyszkowski et al. 2017).

We used LIGER’s large-box framework to produce the mock Euclid catalogues. A schematic diagram representing our workflow is presented in Fig. 1. In this case, LIGER computes both the co-ordinate maps (Eq. (4)) and the magnification (Eq. (7)) starting from the real-space position of the particles in the input N-body simulation. Then the code identifies the snapshots within which the backward light cone of the observer intersects the world lines of the particles after adding the displacement given in Eq. (4). We save the intersection position together with the corresponding magnification and redshift change.

Fig. 1. Flow chart of the LIGER method in the ‘large-box’ mode (top row) and of the clustering analysis performed in this paper (bottom row). Inputs and outputs are displayed as red and blue parallelograms, respectively, while processes are shown as green rectangles.

The updated particle positions are used to compute the matter overdensity in redshift space, δ_s. Eventually, the galaxy distribution in redshift space is obtained using

$$\begin{array}{c}\hfill {\delta }_{\mathrm{g},\mathrm{s}}=(b-1){\delta }_{\mathrm{m}}^{\mathrm{com}}+{\delta }_{\mathrm{s}}+\mathcal{E}\delta lna+\mathcal{Q}(\mathcal{M}-1),\end{array}$$(12)

which matches Eq. (11) under two assumptions. Namely, (i) $|\mathcal{H}\mathrm{\Theta }|/c^2\ll |\delta lna|$; (ii) we can neglect the linear perturbation of $\sqrt{-g}$, where g denotes the determinant of the metric tensor⁹. The neglected terms are only relevant at scales larger than the Hubble radius (for more details, see Borzyszkowski et al. 2017). Aside from these minor contributions, Eq. (12) recovers the theoretical results obtained in Eq. (11).

3.2. N-body simulations

We considered a flat ΛCDM background cosmological model based on the results from the Planck mission (Planck Collaboration VI 2020) with the matter density parameter Ω_m = 0.3158, baryon density parameter Ω_b = 0.0508, and dimensionless Hubble constant h = 0.673. We also assumed that primordial scalar perturbations form a Gaussian random field with a linear power spectrum of power-law shape characterised by the spectral index n_s = 0.966 and the amplitude A_s = 2.1 × 10⁻⁹ (defined at the wavenumber k_* = 0.05 Mpc⁻¹). We computed the matter transfer function using the CAMB code (Lewis & Bridle 2002).

In order to encompass the full EWSS within each of our simulation boxes, we studied structure formation within periodic cubic volumes with a comoving side length of L_box = 12 h⁻¹Gpc. As we are only interested in quasi-linear scales, we used second-order Lagrangian perturbation theory (2LPT) to build the dark-matter distribution that forms the input to LIGER. For this step, we applied the MUSIC code (Hahn & Abel 2011) to 1024³ equal-mass particles, which initially form a regular Cartesian grid. Overdensities were computed with the classical cloud-in-cell scheme using the same grid. The gravitational potential was obtained by solving the Poisson equation with spectral methods (Hockney & Eastwood 1988).

We ran 35 independent N-body simulations and extracted four non-overlapping light cones from each of them, resulting in a total of N_mocks = 140 mock catalogues.

3.3. Euclid Hα galaxies

Applying Eq. (12) to Euclid requires knowledge of the functions b, ℰ, and 𝒬 for the Hα galaxies targeted by the EWSS. In the absence of accurate data, we assume that model 3 in Pozzetti et al. (2016) provides an accurate description of the luminosity function, as is suggested by recent observations (Bagley et al. 2020). Considering a flux limit of F_lim = 2 × 10⁻¹⁶ erg cm⁻² s⁻¹ (Euclid Collaboration: Scaramella et al. 2022), we computed ${\overline{n}}_{\mathrm{g}}$ by integrating the luminosity function and assuming a (uniform) completeness factor of 70%. We derived the evolution and magnification bias factors using Eqs. (9) and (10). Finally, for the linear bias coefficient, we adopted the linear relation b(z) = 1.46 + 0.68 (z − 1) obtained by fitting the data from Table 3 in Euclid Collaboration: Blanchard et al. (2020). All our results are presented in Fig. 2.

Fig. 2. Top: Background number density of galaxies in the EWSS as a function of redshift. Bottom: Corresponding evolution, magnification, and linear-bias parameters.

3.4. Building the mock catalogues

For each of the 140 light cones, we built four different mock catalogues that progressively include an increasing number of relativistic RSDs. We first generated the galaxy distribution in real space (hereafter denoted by ℛ ). Second, we included the RSDs due to the peculiar velocities of the distant galaxies – that is, the terms depending on v_e in Eqs. (4) and (7) – setting, however, 𝒬 and ℰ to zero in Eq. (12). We dubbed the corresponding catalogues 𝒱. Next, we considered all relativistic RSDs except those due to the observer peculiar velocity v_o (from now on the 𝒢 mocks). Finally, we included all terms to generate the 𝒪 catalogues. In the latter set, we assumed that v_o coincides with the peculiar velocity of the Sun as derived from the CMB dipole (Planck Collaboration I 2020).

To produce catalogues of discrete galaxies, we proceeded as follows. Based on δ_g, s and ${\overline{n}}_{\mathrm{g}}$, we first computed the expected number of galaxies, N_g, in each volume element of the light cone. We then drew from a Poisson distribution with mean N_g and randomly distributed the corresponding number of galaxies within the cell.

After taking all these steps, we obtained a full-sky galaxy catalogue covering the redshift range 0.9 < z < 1.8. In order to mimic the expected angular distribution of the EWSS, we masked 20° around the Galactic and Ecliptic planes as shown in Fig. 3¹⁰. We measured clustering statistics in four tomographic redshift bins with boundaries z ∈ {0.9, 1.1, 1.3, 1.5, 1.8} and, particularly for the angular power spectrum, in a broader bin that covers the entire redshift range covered by the mock catalogues. The binning strategy used in this work has been chosen to reproduce that in Euclid Collaboration: Blanchard et al. (2020).

Fig. 3. Projected galaxy number counts for one of the 𝒱 mocks in the z ∈ (0.9, 1.1) redshift bin. The black lines enclose the region used in Sect. 8.2.1 as an example of a simply connected domain.

3.5. Random catalogues

Estimating two-point statistics requires generating unclustered distributions of points with the same angular footprint and radial selection function as the actual galaxy data (see Sect. 6.1 and Sect. 7.1). We built a ‘random catalogue’ for each of our mocks in three steps. We first measured the mean galaxy number counts within radial shells of 20 h⁻¹ Mpc width. We then interpolated the results with a cubic spline to obtain the cumulative redshift distribution of the galaxies. Finally, we used the inverse transform method to pick a redshift for the unclustered points, which are also assigned a random line-of-sight direction. The random catalogues contain five times more galaxies than the original light cones. These catalogues are smaller than those typically employed in power-spectrum estimation. Increasing the size of the random catalogue reduces the shot noise contribution to the statistical error of the power spectrum (Feldman et al. 1994). However, at the large scales considered in this work, the statistical error is dominated by the sample variance. Thus, using a smaller catalogue has minimal impact on our results while significantly reducing computational overhead.

4. Statistical methods

Given a clustering statistic, S, we want to understand whether the contribution of specific relativistic RSDs to the measured signal is detectable or not with the EWSS. For instance, the impact of the peculiar velocity of the observer can be quantified by comparing the clustering statistic extracted from our 𝒪 and 𝒢 mock catalogues. Similarly, by comparing the 𝒪 and 𝒱 light cones we can also study the relevance of the weak lensing contribution.

We denote with D_i^a the n-dimensional (column) data vector containing the measurements of the clustering statistics in a particular mock (characterised by the index, i) of type a ∈ {ℛ , 𝒱, 𝒢, 𝒪}. By taking the expectation over the 140 realisations, we computed the mean signal and the covariance matrix in the measurements,

$$\begin{array}{cc}\hfill {\mathit{\mu }}_a& =\mathbb{E}[{\mathit{D}}_i^a],\hfill \end{array}$$(13)

$$\begin{array}{cc}\hfill {\mathsf{C}}_a& =\mathbb{E}[({\mathit{D}}_i^a-{\mathit{\mu }}_a){({\mathit{D}}_i^a-{\mathit{\mu }}_a)}^{\mathrm{T}}].\hfill \end{array}$$(14)

Based on the Fisher-information matrix, we can then estimate the S/N for the detection of the RSDs that are not included in the a mocks by using (see e.g. Sect. 3.3.3 in Borzyszkowski et al. 2017)

$$\begin{array}{c}\hfill {(\mathrm{S}/\mathrm{N})}^2={({\mathit{\mu }}_{\mathcal{O}}-{\mathit{\mu }}_a)}^{\mathrm{T}}{\mathsf{C}}_{\mathcal{O}}^{-1}({\mathit{\mu }}_{\mathcal{O}}-{\mathit{\mu }}_a)=:m_0,\end{array}$$(15)

where we corrected for the bias of the inverse covariance matrix due to the finite number of realisations used to estimate it (Kaufman 1967; Hartlap et al. 2007).

Classical hypothesis testing based on the likelihood function provides another possibility to quantify the detectability of the different RSD terms. Assuming Gaussian errors, the likelihood that the dataset, D_i^𝒪, is drawn from a model, M_a, with signal μ_a is

$$\begin{array}{c}\hfill L(M_a|{\mathit{D}}_i^{\mathcal{O}})\propto \frac{exp[-{({\mathit{D}}_i^{\mathcal{O}}-{\mathit{\mu }}_a)}^{\mathrm{T}}{\mathsf{C}}_a^{-1}({\mathit{D}}_i^{\mathcal{O}}-{\mathit{\mu }}_a)/2]}{{(2\pi )}^{n/2}det{\mathsf{C}}_a}.\end{array}$$(16)

We use the term ‘model’ to indicate a prediction for S with fixed values of the parameters that describe the galaxy population and the underlying cosmology. No fitting of model parameters is considered here. Basically, a model corresponds to an infinite ensemble of mocks all including the same RSD terms (for instance, the 𝒢 mocks) and is described by the corresponding signal and noise covariance.

We now formulate the null hypothesis, ℋ₀, that the Euclid data are a realisation of model M_a that does not include all RSD terms present in M_𝒪. We want to test this assumption against the alternative hypothesis, ℋ₁, that the data are drawn from M_𝒪. The Neyman–Pearson lemma states that the likelihood-ratio test statistic

$$\begin{array}{c}\hfill {\lambda }_i=2ln\left[\frac{L(M_a|{\mathit{D}}_i^{\mathcal{O}})}{L(M_{\mathcal{O}}|{\mathit{D}}_i^{\mathcal{O}})}\right]\end{array}$$(17)

provides the most powerful test for two simple hypotheses (i.e. with fixed model parameters). The null hypothesis is rejected with confidence level α if λ_i < ω_α, where ω_α is a real number such that the probability 𝒫(λ_i < ω_α|ℋ₀) = 1 − α. Under ℋ₀, 𝒫(λ_i|ℋ₀) is Gaussian, with a mean of m₀ ≥ 0 and variance of s₀² = 4m₀ (see Appendix A in Borzyszkowski et al. 2017). Adopting a 95% confidence level, we thus obtain ${\omega }_{95}=m_0-3.29\sqrt{m_0}$. Similarly, under ℋ₁, λ_i follows a Gaussian distribution with a mean of m = −(μ_𝒪 − μ_a)^TC_a⁻¹(μ_𝒪 − μ_a)≤0 and variance of s² = 4|m|. Therefore, ℋ₀ is rejected in a fraction,

$$\begin{array}{c}\hfill f_{95}=\frac{1}{2}[1+\mathrm{erf}\left(\frac{{\omega }_{95}-m}{2\sqrt{2|m|}}\right)],\end{array}$$(18)

of the realisations. If we neglect the small difference between the covariance matrices, we find that m = −m₀ and s = s₀ = 4m₀¹¹. In this case, the separation between m and m₀ expressed in units of the standard deviation of the distributions is $(m_0-m)/(2\sqrt{m_0})=\sqrt{m_0}$, which coincides with the S/N given in Eq. (15). The fraction of realisations in which the data reject ℋ₀ is thus $f_{95}=\{1+\mathrm{erf}[(\sqrt{m_0}-1.645)/\sqrt{2}]\}/2$. We note that S/N = 1 gives f₉₅ = 0.259, while f₉₅ = 0.5 and 0.9 correspond to S/N = 1.645 and 2.93, respectively.

In the remainder of this paper, we apply these statistical tests to our mock Euclid light cones. As was anticipated at the beginning of this section, by comparing the 𝒪 and 𝒢 sets (hereafter 𝒪–𝒢 test) we assess the detectability of the RSDs generated by the peculiar velocity of the Sun. The goal here is to determine if and how often we manage to reject the null hypothesis that the 𝒪 data are drawn from an ensemble with v_o = 0 . In addition, we apply the likelihood-ratio test to the 𝒪 and 𝒱 mocks (hereafter 𝒪–𝒱 test). In this case, we aim to quantify whether the EWSS can reject the null hypothesis that the influence of the integrated terms (dominated by the weak lensing contribution) on the clustering signal is negligible.

In order to perform these tests, we directly computed histograms of the values of λ_i derived from the 140 𝒪 mock catalogues. In all cases, we only used the C_𝒪 covariance. This operation mimics what is usually done in the analysis of actual surveys, whereby the covariance is estimated from suites of mock catalogues and not changed with the theoretical models for which the likelihood is evaluated. We computed the S/N as

$$\begin{array}{c}\hfill \text{S/N}=2{\widehat{m}}_0/{\widehat{s}}_0,\end{array}$$(19)

where the hat denotes estimates derived from the 140 mocks and determine f₉₅ by counting the number of realisations in which λ_i < ω₉₅. In order to facilitate the understanding of the likelihood-ratio test for two simple hypotheses, we present an illustrative example in Fig. 4.

Fig. 4. Elucidatory example of the statistical test outlined in Sect. 4. In the top panel, the blue histogram on the right-hand side shows the distribution of the likelihood-ratio test statistic λ evaluated from the mock catalogues that do not include all relativistic RSDs (i.e. under the null hypothesis, ℋ0, which we try to rule out using the observed data). The red histogram on the left-hand side, instead, displays the distribution of λ in the mock light cones that account for all effects (i.e. under the alternative hypothesis, ℋ1). The solid and dashed curves represent Gaussian models for the histograms, as is described in the main text. The S/N is a measure of the separation between the two histograms in units of their RMS scatter. The shaded region highlights the realisations in which ℋ0 is ruled out at the 95% confidence level. The bottom panel only differs from the top one in the fact that the covariance matrix C𝒪 has been used to compute all likelihood functions. We use this approximation in the remainder of this paper.

5. Angular power spectrum

In this section, we investigate the impact of relativistic RSDs on the angular power spectrum, C_ℓ.

5.1. Estimator

We made use of the Hierarchical Equal Area isoLatitude Pixelisation (HEALPix; Górski et al. 2005; Zonca et al. 2019) algorithm¹² to partition the sky into N_pix = 12  ×  (1024)² pixels. After dividing our mock light cones into multiple redshift bins, we computed the projected number counts of galaxies in each bin (indicated by the superscript i) and pixel, 𝒩_gⁱ(Ω). A sample sky map with the Euclid mask is shown in Fig. 3 using a cylindrical Cartesian co-ordinate system. The projected density contrast is then

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\mathrm{g}}^i(\mathrm{\Omega })=\frac{{\mathcal{N}}_{\mathrm{g}}^i(\mathrm{\Omega })}{{\overline{\mathcal{N}}}_{\mathrm{g}}^i}-1,\end{array}$$(20)

where ${\overline{\mathcal{N}}}_{\mathrm{g}}^i$ denotes the average of 𝒩_gⁱ(Ω) over the solid angle subtended by the survey.

We expanded the projected overdensities in spherical harmonics,

$$\begin{array}{c}\hfill a_{\ell m}^i={\displaystyle \int {\mathrm{\Sigma }}_{\mathrm{g}}^i}(\mathrm{\Omega })Y_{\ell m}^{\ast }(\mathrm{\Omega }){\mathrm{d}}^2\mathrm{\Omega },\end{array}$$(21)

and measured the angular auto- and cross-power spectra using the pseudo-C_ℓ (PCL) estimator (Peebles 1973; Loureiro et al. 2019)

$$\begin{array}{c}\hfill {\widehat{C}}_{\ell }^{ij}=\frac{1}{w_{\mathrm{p}}^2(2\ell +1)f_{\mathrm{sky}}}{\displaystyle \sum _{m=-\ell }^{\ell }}{a}_{\ell m}^i{a_{\ell m}^j}^{\ast }-\frac{{\delta }_{\mathrm{K}ij}}{{\overline{\mathcal{N}}}_{\mathrm{g}}^i},\end{array}$$(22)

where f_sky denotes the fraction of the sky covered by the survey and w_p is a correction factor due to the finite pixelisation of the sphere (see the HEALPix documentation for more details). In Appendix A, we present a validation test of our pipeline for producing the mock catalogues and measuring the angular power spectra.

5.2. Results

5.2.1. Peculiar velocity of the observer

In Eq. (11), a dipolar pattern is superimposed on the galaxy density contrast whenever the observer is not comoving with the cosmic expansion traced by matter (see also Gibelyou & Huterer 2012). Figure 5 displays how this kinematic dipole alters the angular power spectrum. Shown are the average and root-mean-square (RMS) scatter over the 140 mocks of the difference between the 𝒪 and 𝒢 spectra, which is $\mathrm{\Delta }{\widehat{C}}_{\ell }:={\widehat{C}}_{\ell ,\mathcal{O}}-{\widehat{C}}_{\ell ,\mathcal{G}}$, for two different redshift bins. It is evident that, in a full-sky survey (red symbols), only the dipole (ℓ = 1) is affected by v_o. On the other hand, the correction spreads over all the odd multipoles (ℓ = 3, 5, …) in our Euclid mocks due to the partial sky coverage of the catalogues. Expanding the projected overdensity in spherical harmonics (Eq. (21)) on a partial sky, where the base functions are no longer orthogonal, results in mode-mixing in the estimated angular power spectrum (Peebles 1973). Consequently, the dipole signal leaks into higher odd multipoles.

Fig. 5. Difference between the angular power spectra of the 𝒪 and 𝒢 mocks for a full-sky survey (red) and the EWSS (blue). The symbols indicate the mean signal, while the error bars show the standard deviation over the 140 realisations.

In order to quantify how well the modification due to v_o can be detected with the EWSS, we applied the 𝒪–𝒢 test introduced in Sect. 4 to the first ten multipoles. The resultant S/N values for the auto- and cross-spectra are shown in Table 1. For the tomographic redshift bins, the excess clustering induced by v_o can be barely identified with a S/N of order one. This characteristic signature becomes much more discernible in the auto-correlation function evaluated after projecting the galaxies from the broad redshift interval z ∈ (0.9, 1.8). In this case, we obtain S/N = 2.1.

5.2.2. Weak lensing

Figure 6 displays the angular auto- and cross-power spectra extracted from the 𝒱 and 𝒪 mocks and re-binned with Δℓ = 5. The symbols show the average signal over the 140 realisations and the shaded region indicates the RMS scatter for the 𝒪 measurements. Since 2LPT at z ≃ 1 underestimates the non-linear matter power spectrum for wavenumbers k ≳ 0.05 h Mpc⁻¹ (e.g. Taruya et al. 2018), we only consider harmonics of degree ℓ ≤ 120. The positioning of the panels is as in Tables 1 and 2. The cross-spectra between non-overlapping redshift bins are consistent with zero for the 𝒱 mocks and show a positive clustering signal for the 𝒪 catalogues. This difference is due to the integral terms in Eq. (11), in particular to the dominant weak lensing contribution. The auto-spectra also show enhanced clustering for the 𝒪 mocks. This is particularly evident in the bins that include galaxies with the highest redshifts.

The S/N values obtained with the 𝒪–𝒱 test are reported in Table 2. The detection of the lensing term is highly significant in the cross-correlations of well separated tomographic redshift bins and in all statistics involving the wide bin z ∈ (0.9, 1.8).

6. Two-point correlation function

The 2PCF is one of the most employed summary statistics to extract cosmological information from the large-scale structure of the Universe. In terms of the galaxy overdensity, it can be defined as (Peebles 1980)

$$\begin{array}{c}\hfill \langle {\delta }_{\mathrm{g}}({\mathit{x}}_1){\delta }_{\mathrm{g}}({\mathit{x}}_2)\rangle ={\xi }_{\mathrm{g}}({\mathit{x}}_1,{\mathit{x}}_2),\end{array}$$(23)

where the brackets denote the average over an ensemble of realisations. In real space, we expect that δ_g is a statistically homogeneous and isotropic random field, so that ξ_g depends only on the magnitude of the separation between the points at which it is evaluated. Redshift-space distortions, however, break the translational and rotational symmetry of the 2PCF (as several terms in Eq. (11) depend on the line-of-sight direction with respect to the observer) into an azimuthal symmetry with respect to the line of sight. The 2PCF in redshift space thus depends on the shape and size of the (possibly non-Euclidean) triangle formed by the observer and the galaxy pair (Szalay et al. 1998; Matsubara 2000a).

Fig. 6. Mean auto-and cross-angular power spectra extracted from the 𝒪 (red stars) and 𝒱 (green circles) mocks. The shaded region indicates the RMS scatter of the 𝒪 spectra. The panels are ordered as the entries in Tables 1 and 2. The cross-spectra of the tomographic redshift bins are multiplied by ten to improve the readability of the figure.

Fig. 7. Mean ℓ = 0, 2, and 4 multipoles of the 2PCF measured from the ℛ (black squares), 𝒱 (green circles), 𝒢 (blue crosses), and 𝒪 (red stars) mocks in the four tomographic redshift bins. The shaded areas highlight the RMS scatter among the 𝒪 light cones.

Fig. 8. PDF of the likelihood-ratio test statistic forthe 𝒪–𝒱 test applied to the multipoles of the 2PCF.

6.1. Estimator

We used the Landy–Szalay (LS) estimator (Landy & Szalay 1993) implemented in the official Euclid code (Euclid Collaboration: De la Torre et al. in prep.), which adopts the midpoint co-ordinate system

$$\begin{array}{c}\hfill \mathit{r}={\mathit{x}}_2-{\mathit{x}}_1,{\mathit{x}}_{\mathrm{m}}=\frac{{\mathit{x}}_2+{\mathit{x}}_1}{2},\end{array}$$(24)

and defines the pair-orientation angle, φ, through

$$\begin{array}{c}\hfill \mu =cos\phi =\widehat{\mathit{r}}·{\widehat{\mathit{x}}}_{\mathrm{m}},\end{array}$$(25)

where the hat denotes the unit vector; in other words, $\widehat{\mathit{r}}:=\mathit{r}/r$¹³. After averaging over x_m within the volume of interest, the code outputs an estimate of the 2PCF as a function of the galaxy separation and orientation with respect to the local line of sight, $\widehat{\xi }(r,\mu )$. It also computes the 2PCF Legendre multipoles

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

where ℒ_ℓ(μ) denotes the Legendre polynomial of degree, ℓ. In what follows, we only consider the monopole (ℓ = 0), quadrupole (ℓ = 2), and hexadecapole (ℓ = 4) moments. These are the only non-vanishing moments in Kaiser’s GPP model and dominate the signal in general.

We originally estimated them in 500 equally spaced bins covering the range r ∈ [20 h⁻¹ Mpc,  Δx], where Δx denotes the comoving radial width of each redshift bin. This choice prevents the clustering signal from being dominated by galaxy pairs with particular angular separations, and thus limits window-function effects. We re-binned our results in different ways depending on our applications. For instance, we used ten equally separated bins in the range r ∈ [35 h⁻¹ Mpc,  Δx] for each multipole to perform the 𝒪–𝒢 and 𝒪–𝒱 tests.

6.2. Results

Figure 7 shows the mean multipoles of the 2PCF obtained from the different sets of mock catalogues. The shaded areas indicate the scatter for the 𝒪 light cones. The real-space monopole moment is positive at small separations, presents the baryonic-oscillation feature at r ≃ 100 h⁻¹ Mpc, and crosses zero at r ≃ 125 h⁻¹ Mpc, while ${\widehat{\xi }}_2$ and ${\widehat{\xi }}_4$ vanish as expected. Redshift-space distortions enhance the clustering signal in ${\widehat{\xi }}_0$ and generate a negative ${\widehat{\xi }}_2$ and a positive ${\widehat{\xi }}_4$.

6.2.1. Peculiar velocity of the observer

The clustering signal extracted from the 𝒪 (red stars) and 𝒢 (blue crosses) mocks is hardly distinguishable at all scales in all tomographic redshift bins. This visual impression is confirmed by the 𝒪–𝒢 test, which consistently gives values of S/N ≤ 1 (see Table 3) for all redshift bins. At first sight, this appears to be at odds with the results by Bertacca et al. (2020) who predict stronger corrections due to v_o. However, this study pushes the analysis to larger separations than ours and uses the bisector convention to define the line of sight to a galaxy pair (which changes the multipoles, e.g. Raccanelli et al. 2014; Reimberg et al. 2016).

6.2.2. Weak lensing

Comparing the mean signal from the 𝒢 (blue crosses) and 𝒱 (green circles) multipoles, we notice that their difference increases with redshift and the pair separation. The contribution of the integral terms always enhances the clustering signal in the quadrupole and hexadecapole moments, but by an amount that is relatively small compared to the scatter in the measurements.

Performing the 𝒪–𝒱 test, we find that the S/N steadily grows from 0.81 to 2.48 from the lowest to the highest redshift bin (see Fig. 8 and Table 3). In the latter, the likelihood-ratio test manages to reject the velocity-only model for RSDs at the 95% confidence level in 78% of our mock catalogues.

7. Power spectrum

The galaxy power spectrum is the workhorse of cosmological-parameter inference. By analogy with Sect. 6, we introduced the covariance between two Fourier modes of the galaxy overdensity, $⟨{\stackrel{\sim }{\delta }}_{\mathrm{g}}(\mathit{k}){\stackrel{\sim }{\delta }}_{\mathrm{g}}(\mathit{k}\prime )⟩={(2\pi )}^{3}C(\mathit{k},\mathit{k}\prime )$. If, in real space, δ_g(x) is a statistically homogeneous and isotropic field, then only the diagonal part of the covariance does not vanish, and the galaxy power spectrum, P_g(k), can be introduced using the relation $⟨{\stackrel{\sim }{\delta }}_{\mathrm{g}}(\mathit{k}){\stackrel{\sim }{\delta }}_{\mathrm{g}}(\mathit{k}\prime )⟩={(2\pi )}^3{\delta }_{\mathrm{D}}(\mathit{k}+\mathit{k}\prime )P_{\mathrm{g}}(k)$. However, in redshift-space, where statistical homogeneity is lost, the covariance is not diagonal and the definition above does not apply (Zaroubi & Hoffman 1996). A different approach is thus needed.

The ‘local’ power spectrum, P_loc, was obtained by Fourier transforming the 2PCF with respect to r (Scoccimarro 2015),

$$\begin{array}{cc}\hfill P_{\mathrm{loc}}({\mathit{x}}_{\mathrm{m}},\mathit{k})& :={\displaystyle \int \langle {\delta }_{\mathrm{g}}({\mathit{x}}_{\mathrm{m}}+\mathit{r}/2){\delta }_{\mathrm{g}}({\mathit{x}}_{\mathrm{m}}-\mathit{r}/2)\rangle }{\mathrm{e}}^{-\mathrm{i}\mathit{k}·\mathit{r}}{\mathrm{d}}^{3}r\hfill \end{array}$$(27)

$$\begin{array}{cc}& ={\displaystyle \int C}(-\mathit{k}+\mathit{q}/2,\mathit{k}+\mathit{q}/2){\mathrm{e}}^{\mathrm{i}\mathit{q}·{\mathit{x}}_{\mathrm{m}}}{\mathrm{d}}^{3}q.\hfill \end{array}$$(28)

In order to compress the clustering information into a set of functions of the wavenumber k, it is convenient to expand P_loc in Legendre polynomials of $\widehat{\mathit{k}}·{\widehat{\mathit{x}}}_{\mathrm{m}}$ and average over x_m. These operations yield the so-called multipole moments of the power spectrum,

$$\begin{array}{c}\hfill P_{\ell }(k):=(2\ell +1){\displaystyle \int \int P_{\mathrm{loc}}}({\mathit{x}}_{\mathrm{m}},\mathit{k}){\mathcal{L}}_{\ell }(\widehat{\mathit{k}}·{\widehat{\mathit{x}}}_{\mathrm{m}})\frac{{\mathrm{d}}^2{\mathrm{\Omega }}_k}{4\pi }\frac{{\mathrm{d}}^3{x}_{\mathrm{m}}}{V},\end{array}$$(29)

where V denotes the volume under consideration. In this section, we investigate the impact of relativistic RSDs on the power spectrum multipoles measured by the Euclid survey.

7.1. Estimator

In order to measure the multipoles of the power spectrum from a galaxy redshift survey, the ensemble average in Eq. (27) was replaced with a mean over a set of Fourier modes. We used the Yamamoto–Bianchi (Bianchi et al. 2015; Scoccimarro 2015) estimator implemented in the official Euclid code (Euclid Collaboration: Salvalaggio et al. in prep.).

Following Feldman et al. (1994, hereafter FKP), we first built the weighted galaxy overdensity,

$$\begin{array}{cc}\hfill F\left(\mathit{x}\right)& =\frac{w\left(\mathit{x}\right)}{\sqrt{A}}[{\widehat{n}}_{\mathrm{g}}\left(\mathit{x}\right)-\alpha {\widehat{n}}_{\mathrm{r}}(\mathit{x})],\hfill \end{array}$$(30)

by comparing the observed number density of galaxies ${\widehat{n}}_{\mathrm{g}}(\mathit{x})$ and its counterpart in the corresponding random catalogue, ${\widehat{n}}_{\mathrm{r}}(\mathit{x})$, introduced in Sect. 3.5. In order to minimise the variance of the measured multipoles, we adopted the weight function

$$\begin{array}{c}\hfill w(\mathit{x})=\mathcal{I}(\mathit{x}){[1+{\overline{n}}_{\mathrm{g}}(x){\mathcal{P}}_0]}^{-1},\end{array}$$(31)

where ${\overline{n}}_{\mathrm{g}}$ denotes the mean density of galaxies, ℐ(x) is an indicator function that is one inside the volume under study and zero elsewhere, and the parameter 𝒫₀ = 2 × 10⁴ h⁻³ Mpc³ gives an approximate value for the galaxy power spectrum at the scales of interest. The normalisation factor, A, was determined through the integral $A={\displaystyle \int {\mathit{w}}^2}{{\overline{n}}^2}_{\mathrm{g}}{\mathrm{d}}^{3}x$. Finally, the rescaling factor, α, was calculated using

$$\begin{array}{c}\hfill \alpha =\frac{{\displaystyle \int w}(\mathit{x}){\widehat{n}}_{\mathrm{g}}(\mathit{x}){\mathrm{d}}^{3}x}{{\displaystyle \int w}(\mathit{x}){\widehat{n}}_{\mathrm{r}}(\mathit{x}){\mathrm{d}}^{3}x}.\end{array}$$(32)

The multipoles of the galaxy power spectrum with respect to the local line-of-sight direction could be directly estimated by computing (Yamamoto et al. 2000)

$$\begin{array}{cc}\hfill {\widehat{P}}_{\ell }(k)& =(2\ell +1){\displaystyle \int \int \int [F({\mathit{x}}_1)F({\mathit{x}}_2){\mathrm{e}}^{-\mathrm{i}\mathit{k}·({\mathit{x}}_{\mathbf{1}}-{\mathit{x}}_{\mathbf{2}})}}\hfill \\ \hfill & {\mathcal{L}}_{\ell }(\widehat{\mathit{k}}·{\widehat{\mathit{x}}}_{\mathrm{m}}){\mathrm{d}}^3{x}_1{\mathrm{d}}^3{x}_2\frac{{\mathrm{d}}^2{\mathrm{\Omega }}_k}{4\pi }]-{\widehat{P}}_{\ell }^{\mathrm{SN}}(k),\hfill \end{array}$$(33)

where Ω_k denotes the solid angle in Fourier space and the shot noise contribution, ${\widehat{P}}_{\ell }^{\mathrm{SN}}(k)$,. is given by

$$\begin{array}{c}\hfill {\widehat{P}}_{\ell }^{\mathrm{SN}}(k)=\frac{(1+\alpha )}{A}{\displaystyle \int w^2}(\mathit{x}){\widehat{n}}_{\mathrm{g}}(\mathit{x}){\mathcal{L}}_{\ell }(\widehat{\mathit{k}}·{\widehat{\mathit{x}}}_{\mathrm{m}}){\mathrm{d}}^{3}x.\end{array}$$(34)

Yamamoto et al. (2006) noticed that replacing the factor ${\mathcal{L}}_{\ell }(\widehat{\mathit{k}}·{\widehat{\mathit{x}}}_{\mathrm{m}})$ with either ${\mathcal{L}}_{\ell }(\widehat{\mathit{k}}·{\widehat{\mathit{x}}}_1)$ or ${\mathcal{L}}_{\ell }(\widehat{\mathit{k}}·{\widehat{\mathit{x}}}_2)$ in Eq. (33) results in a tremendous speed-up of the estimator. With this substitution – nowadays known as the LPP approximation – in fact, ${\widehat{P}}_{\ell }(k)$ can be written as the product of two Fourier transforms that can be conveniently evaluated using the FFT algorithm (Beutler et al. 2014; Bianchi et al. 2015; Scoccimarro 2015); that is,

$$\begin{array}{c}\hfill {\widehat{P}}_{\ell }(k)=(2\ell +1){\displaystyle \int F_{\ell }}(\mathit{k})F_0(-\mathit{k})\frac{{\mathrm{d}}^2{\mathrm{\Omega }}_k}{4\pi }-{\widehat{P}}_{\ell }^{\mathrm{SN}}(k),\end{array}$$(35)

where

$$\begin{array}{c}\hfill F_{\ell }(\mathit{k}):={\displaystyle \int F}(\mathit{x}){\mathrm{e}}^{-\mathrm{i}\mathit{k}·\mathit{x}}{\mathcal{L}}_{\ell }(\widehat{\mathit{k}}·\widehat{\mathit{x}}){\mathrm{d}}^{3}x.\end{array}$$(36)

The official Euclid code that we used implements the faster estimator. It is worth stressing that Eqs. (33) and (35) define two different statistics, which generate different outputs when applied to wide-angle surveys (e.g. Samushia et al. 2015). For example, it has been demonstrated that the end-point convention distorts the values of the multipoles when compared to those evaluated using the midpoint convention (Reimberg et al. 2016; Castorina & White 2018).

As is shown in Figs. 7 and 9, the effects we investigate here leave their imprints on extremely large scales. For this reason, we computed FFTs within (periodic) cubic boxes with a side length of L_FFT = 16 015 h⁻¹ Mpc. That choice allows us to measure the power-spectrum multipoles down to the fundamental frequency k_F = 4 × 10⁻⁴ h Mpc⁻¹. We originally measured the multipoles with ℓ = 0, 2, and 4 within bins of size Δk = k_F and we re-binned them in different ways depending on usage. For the 𝒪–𝒢 and 𝒪–𝒱 tests, we employed 12 equally spaced bins in the range k ∈ [0.4, 36]×10⁻³ h Mpc⁻¹.

Fig. 9. Mean ℓ = 0, 2, and 4 multipoles of the power spectrum measured from the 𝒱 (green circles), 𝒢 (blue crosses), and 𝒪 (red stars) mocks in the four tomographic redshift bins. The shaded areas highlight the RMS scatter among the 𝒪 light cones.

7.2. Results

The mean multipoles of the power spectrum measured from the different sets of mock catalogues are shown in Fig. 9 together with the RMS scatter from the 𝒪 set. The monopole moment shows the largest clustering amplitude and is measured with a high S/N, particularly for k > 2 × 10⁻² h Mpc⁻¹. At the opposite extreme, the hexadecapole moment is suppressed by an order of magnitude with respect to P₀ and its measurements are very noisy at all scales. The quadrupole moment has intermediate properties between the other two.

7.2.1. Peculiar velocity of the observer

By comparing the 𝒪 and 𝒢 spectra, we observe that the peculiar velocity of the observer modifies all multipoles at extremely large scales. For the monopole, this is consistent with the results presented by Elkhashab et al. (2021) who showed that a non-vanishing v_o adds an oscillatory signal (damped with increasing k) to P₀ with an oscillation frequency that increases with the characteristic redshift of the galaxy population. The coarse k-binning we adopt in this work does not reveal the details of the oscillations that then appear as a large-scale boost of the clustering amplitude in Fig. 9. Similar distortions are clearly noticeable also in P₂ and P₄¹⁴. For the latter, in the first three tomographic bins, the corrected signal becomes negative at the largest scales probed here. The chances to detect these signatures with the EWSS are meagre, however, given the large scatter in the measurements at small wavenumbers. The S/N obtained with the 𝒪–𝒢 test is 1.3 at best (see Table 4) since the corrections due to v_o are localised on very large scales where the measurement noise is large. Still, the enhanced clustering could bias measurements of the local-non-Gaussianity parameter, f_NL, based on P₀.

7.2.2. Weak lensing

The difference between the multipoles extracted from the 𝒢 and 𝒱 mocks is minimal in the first three tomographic redshift bins but becomes more pronounced in the last one, particularly for ℓ = 2 and 4. This is consistent with the expectation that weak gravitational lensing should have a larger impact on the clustering of high-redshift galaxies. The 𝒪–𝒱 test closely mirrors the results we obtained for the 2PCF multipoles: the S/N for the detection of the integrated RSDs is always around one in the first three bins and jumps to 2.4 in the last one, for which the velocity-only model for RSDs is ruled out with 95% confidence in 78% of the realisations.

8. Perturbative models and survey window function

Up to this point, we have estimated the relative importance of various types of RSD in the forthcoming Euclid data by comparing the outputs of our different suites of mock catalogues. This approach deviates from what is usually done to interpret clustering measurements from redshift surveys. The standard procedure is to compare the observed summary statistics to analytical models (based on some flavour of perturbation theory) after accounting for the window function of the survey. In this section, we pursue this approach and assess whether Kaiser’s model for RSDs is accurate enough to describe the large-scale limit of the power-spectrum multipoles that Euclid will measure. We also present some interesting findings about the properties of the Euclid window function and the possibility of using disjoint patches of the sky in the same measurement of ${\widehat{P}}_{\ell }(k)$.

8.1. The Kaiser model

In a seminal paper, Kaiser (1987) presented a theoretical model for the galaxy power spectrum in redshift space based on linear perturbation theory. The model only considers RSDs arising from the peculiar velocity gradient in Eq. (11). It also relies on the GPP approximation, according to which the lines of sight to all galaxies are parallel (as is expected for a small volume that is located at a large distance from the observer and that is thus seen under a narrow solid angle, i.e. under the distant-observer approximation). At a fixed cosmological redshift, it gives

$$\begin{array}{c}\hfill P_{\mathrm{g}}(\mathit{k})={\displaystyle \sum _{\ell =0,2,4}}{\mathcal{F}}_{\ell }{b}^2{D}_{+}^2{P}_{\mathrm{m}}(k){\mathcal{L}}_{\ell }(\widehat{\mathit{k}}·{\widehat{\mathit{x}}}_{\mathrm{g}}),\end{array}$$(37)

where D₊ denotes the linear growth factor for matter perturbations, P_m(k) is their linear power spectrum at z = 0, ${\widehat{\mathit{x}}}_{\mathrm{g}}$ is the line-of-sight direction, and

$$\begin{array}{cc}\hfill {\mathcal{F}}_0& =1+\frac{2}{3}\beta +\frac{1}{5}{\beta }^2,\hfill \end{array}$$(38)

$$\begin{array}{cc}\hfill {\mathcal{F}}_2& =\frac{4}{3}\beta +\frac{4}{7}{\beta }^2,\hfill \end{array}$$(39)

$$\begin{array}{cc}\hfill {\mathcal{F}}_4& =\frac{8}{35}{\beta }^2,\hfill \end{array}$$(40)

in terms of the (redshift-dependent) linear RSD parameter

$$\begin{array}{c}\hfill \beta =\frac{1}{b}{\frac{\mathrm{d}ln{D}_{+}}{\mathrm{d}lna}|}_{a={(1+z)}^{-1}}.\end{array}$$(41)

Fig. 10. Similar to Fig. 8 using the power spectrum multipoles as the data vector.

This result can be generalised to model observations taken on a section of the past light cone (with volume V_s) of an observer, obtaining (to first approximation, e.g. Yamamoto et al. 1999; Pryer et al. 2022)

$$\begin{array}{cc}\hfill P_{\mathrm{g}}& (\mathit{k})={\displaystyle \sum _{\ell =0,2,4}}{\overline{\mathcal{F}}}_{\ell }{P}_{\mathrm{m}}(k){\mathcal{L}}_{\ell }(\widehat{\mathit{k}}·{\widehat{\mathit{x}}}_{\mathrm{g}}),\hfill \end{array}$$(42)

where

$$\begin{array}{c}\hfill {\overline{\mathcal{F}}}_{\ell }=\frac{{\displaystyle {\int }_{V_{\mathrm{s}}}}{\mathcal{F}}_{\ell }{b}^2{D}_{+}^2{\overline{n}}_{\mathrm{g}}^2\mathrm{d}V}{{\displaystyle {\int }_{V_{\mathrm{s}}}}{\overline{n}}_{\mathrm{g}}^2\mathrm{d}V}.\end{array}$$(43)

Going beyond this model requires dropping the GPP approximation and/or accounting for all the RSD terms appearing in Eq. (11). In order to correct the GPP predictions, Castorina & White (2018) propose expanding the multipoles of the wide-angle power spectrum in the parameter (k x_m)⁻¹. In this framework, the additional relativistic RSDs can then be treated perturbatively (Beutler et al. 2019; Castorina & Di Dio 2022; Noorikuhani & Scoccimarro 2023). The expansion in (k x_m)⁻¹, however, might become inaccurate for large angular separations. Another possible approach is to use the ‘spherical-Fourier-Bessel’ formalism that was introduced by Peebles (1973), extended to redshift space by Heavens & Taylor (1995), and applied to survey data in Percival et al. (2004). The inclusion of GR effects in this formalism is discussed in Yoo & Desjacques (2013), Bertacca et al. (2018), and Semenzato et al. (2025).

8.2. The window convolution and integral constraint

In order to compare theoretical models to the multipoles estimated from a survey, one needs to account for the fact that only a finite volume is observed. We can gain some insight into this issue by first considering the FKP power-spectrum estimator, ${\widehat{P}}_{\mathrm{obs}}(\mathit{k})={|\stackrel{\sim }{F}(\mathit{k})|}^2$, derived under the GPP approximation. In this case, we obtain (Peacock 1991)

$$\begin{array}{cc}\hfill P_{\mathrm{obs}}(\mathit{k})={\displaystyle \int P_{\mathrm{g}}}(\mathit{k}\prime )& |\stackrel{\sim }{W}{(\mathit{k}-\mathit{k}\prime )|}^2{\mathrm{d}}^{3}k\prime \hfill \\ \hfill & -\frac{|\stackrel{\sim }{W}{(\mathit{k})|}^2}{|\stackrel{\sim }{W}{(0)|}^2}{\displaystyle \int P_{\mathrm{g}}}(\mathit{k}\prime ){|\stackrel{\sim }{W}(\mathit{k}\prime )|}^2{\mathrm{d}}^{3}k\prime ,\hfill \end{array}$$(44)

where $P_{\mathrm{obs}}(\mathit{k})=⟨{\widehat{P}}_{\mathrm{obs}}(\mathit{k})⟩$, and $\stackrel{\sim }{W}(\mathit{k})$ is the Fourier transform of the survey window function

$$\begin{array}{c}\hfill W(\mathit{x})=\frac{\alpha w(\mathit{x}){\widehat{n}}_{\mathrm{r}}(\mathit{x})}{\sqrt{A}}.\end{array}$$(45)

The first term in Eq. (44) shows that the power-spectrum estimator mixes the contributions from Fourier modes that differ by less than the characteristic width of $\stackrel{\sim }{W}(\mathit{k})$, which is of the order of V^−1/3 (in accordance with the uncertainty relation between conjugate Fourier variables). This has two main consequences: (i) the power spectrum is substantially distorted on large scales and (ii) additional anisotropy is generated on top of RSDs and the Alcock-Paczynski effect because of the spherical asymmetry of the window function. The second term in Eq. (44) gives the so-called (global) integral constraint, which arises from the assumption that the average density within the survey coincides with the actual mean density of the Universe. This term subtracts the actual power at k = 0, which leaks to larger wavenumbers because of the convolution with the window function, and thus enforces that P_obs(k) is zero for k = 0.

The considerations above can be generalised to the estimator for the power-spectrum multipoles introduced in Eq. (35) based on the LPP approximation. This gives (Beutler et al. 2014; Wilson et al. 2017; Beutler & McDonald 2021)

$$\begin{array}{cc}\hfill P_{\mathrm{obs},\ell }(k)=& {\displaystyle \sum _{\ell \prime =0,2,4}}{\displaystyle {\int }_0^{\infty }}{k}^{\prime 2}{\mathcal{W}}_{\ell \ell \prime }(k,k\prime )P_{\mathrm{g},\ell \prime }(k\prime )\mathrm{d}k\prime \hfill \\ \hfill & -\frac{P_{\ell }^W(k)}{P_0^W(0)}{\displaystyle \sum _{\ell \prime =0,2,4}}{\displaystyle {\int }_0^{\infty }}{k}^{\prime 2}{\mathcal{W}}_{\ell \ell \prime }(0,k\prime )P_{\mathrm{g},\ell \prime }(k\prime )\mathrm{d}{k}^{\prime },\hfill \end{array}$$(46)

where P_{g , ℓ′}(k′) is computed under the GPP approximation, $P_{\ell }^W$ denotes the multipole spectral moments of the window function, and 𝒲_ℓℓ′(k, k′) are the elements of the so-called mixing matrix. The latter can be obtained from the multipoles, Q_ℓ(r), of the 2PCF of W(x) using

$$\begin{array}{cc}\hfill {\mathcal{W}}_{\ell \ell \prime }(k,k\prime )=\frac{2{\mathrm{i}}^{\ell \prime -\ell }}{\pi }(2\ell +1)& {\displaystyle \sum _{L=0}^{\infty }}{\left(\begin{array}{ccc}\ell & \ell \prime & L\\ 0& 0& 0\end{array}\right)}^2\hfill \\ \hfill & \times {\displaystyle {\int }_0^{\infty }}{r}^2{j}_{\ell }(kr)j_{\ell \prime }(k\prime r)Q_L(r)\mathrm{d}r,\hfill \end{array}$$(47)

where the brackets denote the Wigner 3j symbols and j_ℓ is the spherical Bessel function of order ℓ. The functions Q_ℓ(r) could be directly estimated by counting pairs in the random catalogue (Wilson et al. 2017). However, it is computationally faster to obtain them using a Hankel transform,

$$\begin{array}{c}\hfill Q_{\ell }(r)=\frac{{\mathrm{i}}^{\ell }}{2{\pi }^2}{\displaystyle {\int }_0^{\infty }}{k}^2{P}_{\ell }^W(k)j_{\ell }(kr)\mathrm{d}k,\end{array}$$(48)

evaluated with the FFTLog method (Hamilton 2000).

8.2.1. Disconnected patches

The EWSS avoids the ecliptic and Galactic planes and is thus composed of four disconnected regions (see Fig. 3). Considering all of them together would allow us to measure galaxy clustering on the largest possible scales. This is what we did in Sect. 7 to test the importance of the different RSD terms. However, this procedure differs from what is regularly done in ground-based surveys where northern and southern Galactic caps are analysed separately since each of them is subject to different angular systematics.

In order to contrast these approaches, in Fig. 11, we compare the functions Q₀,  Q₂, and Q₄ derived for the full Euclid footprint (solid lines) and for the simply connected patch covering 2565 deg² displayed in Fig. 3 (dashed lines). In both cases, we consider the redshift bin z ∈ (1.1, 1.3). For separations of a few h⁻¹ Mpc, the monopole moment approaches one while Q₂ and Q₄ are close to zero, reflecting the fact that Q(r) is nearly constant and isotropic for spatial lags that are well ’contained’ within the survey window function. On the other hand, Q(r) vanishes when r is larger than the maximum distance between two galaxies in the patch. In between these two asymptotic regimes, the 2PCF of the random catalogue becomes highly anisotropic with |Q₂| and Q₄ that assume values larger than Q₀. The difference between the two survey footprints becomes evident on the largest scales: for instance, in the Euclid case, Q₄ becomes negative when 1600 h⁻¹ Mpc ≲ r ≲ 4400 h⁻¹ Mpc and shows a prominent peak for r ≃ 5000 h⁻¹ Mpc due to the presence of separate patches.

Fig. 11. 2PCF multipoles of the window function for the EWSS (solid) in the redshift bin z ∈ (1.1, 1.3) and for the single patch (dashed) shown in Fig. 3 in the same redshift range. The vertical lines indicate the characteristic length scale V1/3 for the two surveys. The bottom panel offers a more detailed view of the large-scale behaviour of the functions.

In Fig. 12, we present all the elements of the resultant mixing matrices for the two different geometries considering the second tomographic bin z ∈ (1.1, 1.3). We fix k = 0.03 h Mpc⁻¹ and plot 𝒲_ℓℓ′(k, k′) as a function of k′. Each panel refers to a specific combination of ℓ and ℓ′. We plot the results for the Euclid footprint in black, while those for the connected patch are shown in red. As was expected, the most significant contributions come from k′=k but mixing takes place within a relatively broad range of wavenumbers. The matrix elements are always positive for ℓ = ℓ′ and show an oscillatory behaviour otherwise. These overall trends are present for both survey footprints. However, the results for the EWSS show high-frequency oscillations that are absent in the single patch case. This modulation can be traced back to presence of the peak at very large separations in the Q_ℓ functions presented in Fig. 11.

Fig. 12. Elements of the mixing matrix for the two surveys discussed in Fig. 11 evaluated at k = 0.03 h Mpc−1.

8.2.2. Comparison with the LIGER mocks

We finally compare the theoretical model that accounts for the survey window and the integral constraint to the spectral multipoles extracted from our mock catalogues. For the sake of clarity, we summarise here how the theoretical prediction is computed: (i) we generated the linear matter power spectrum at z = 0 using the CAMB code; (ii) we utilised Kaiser’s model to obtain the galaxy power spectrum in redshift space and resorted to Eq. (42) in order to get the observed spectra on the past light cone; (iii) we computed the elements of the mixing matrix using the 2PCF multipoles estimated from the random catalogue of the mocks (see Eq. (47)); and (iv) we accounted for the survey window function and the integral constraint by making use of Eq. (46).

In Fig. 13, we compare the model (solid lines) to the measurements extracted from the 𝒱 mocks (symbols with error bars). As a reference, we also plot the theoretical predictions obtained without accounting for the finite-volume effects (dotted lines). Since we obtained similar results in all tomographic bins, we only show them for the interval z ∈ (1.1 , 1.3). The model and the measurements are in extremely good agreement at all wavenumbers, k > k_C = 2 π/V^1/3 (highlighted with a vertical dashed line). On larger scales, the quadrupole and hexadecapole moments show a small discrepancy, possibly due to wide-angle corrections missed by Kaiser’s model. A side note is in order here. The fact that Kaiser’s model nicely matches our 𝒱 mocks does not imply that it accounts for all RSDs generated by peculiar velocities. In fact, the 𝒱 catalogues were constructed by setting the magnification and evolution bias to zero, which automatically cancels some wide-angle effects. Anyway, the remarkable consistency we find validates the accuracy of the LIGER mock catalogues and demonstrates that the mixing-matrix formalism can be reliably used for surveys made of disconnected patches.

Fig. 13. Power-spectrum multipoles for the 𝒱 mocks in the redshift range z ∈ (1.1, 1.3). The mean and RMS scatter over the 140 mocks are displayed with squares and error bars, respectively. The solid and dotted lines show the theoretical predictions obtained with the Kaiser model (see Eq. (42)) with and without accounting for the window function of the EWSS (see Eq. (46)), respectively. The vertical dashed line marks the characteristic wavenumber, kC = 2 π/V1/3.

Finally, we compared the theoretical model with the multipole moments extracted from the 𝒪 mocks, which mimic what will be observed by the actual Euclid mission. This way, we sought to assess the significance with which a theoretical model that accounts only for the velocity-gradient term can be ruled out due to the presence of additional RSDs in the data. Given the close match between the model and the numerical results obtained from the 𝒱 mocks, we obviously expected very similar results to those presented in Sect. 7.2.2. We employed the likelihood-ratio test and only considered wavenumbers k < k_C grouped in 11 bins covering the range k ∈ [3.2, 36]×10⁻³ h Mpc⁻¹ (slightly narrower than in Sect. 7.2.2). The resulting S/Ns for this (𝒪− P_obs) test are [1.3,  1.2,  1.5,  2.9] for the four tomographic bins listed in the same order as in Table 4. As was expected, these results are in good agreement with those deriving from the 𝒪–𝒱 test. This reinforces the conclusion that it is necessary to include at least the weak lensing corrections in order to obtain a good fit to the Euclid power-spectrum multipoles on scales k < 0.03 h Mpc⁻¹.

9. Summary

The EWSS will measure redshifts and angular positions of nearly 30 million emission-line galaxies over a third of the sky and in the redshift range 0.9 ≲ z ≲ 1.8 (Laureijs et al. 2011). If observational systematic effects are kept under control, it will provide us with the opportunity to study galaxy clustering on unprecedentedly large scales. Theoretical considerations suggest that relativistic ‘projection’ effects alter the clustering signal we measure. In fact, galaxy observables (i.e. redshift, fluxes, and sky position) are influenced by the inhomogeneities that light encounters as it propagates from the galaxies to the observer. However, we construct maps of the galaxy distribution by assuming that galaxies populate an unperturbed model universe (for instance, in order to convert redshift into distances). This discrepancy leads to RSDs, which we need to model in order to interpret clustering data. While the leading effect generated by the relative peculiar velocities between the sources and the observer (Kaiser 1987) is regularly accounted for in all studies, additional contributions due to Doppler, gravitational lensing, Sachs–Wolfe, and time-delay type terms have mostly been neglected. The importance of these relativistic corrections is expected to grow as we investigate physical length scales approaching the Hubble radius. In this work, we make forecasts for their impact on different galaxy two-point statistics that will be measured from the EWSS; namely, the angular power spectrum, the multipole moments of the anisotropic 2PCF, and the multipole moments of the anisotropic power spectrum.

We used the LIGER method (summarised schematically in Fig. 1) to build 140 mock galaxy catalogues, each covering the whole footprint (and redshift range) of the EWSS and accounting for all relativistic RSDs up to the first order in the cosmological perturbations (see Fig. 3 for one example shown in projection). We studied the impact of different relativistic RSDs by turning on and off particular effects while producing the galaxy catalogues. In particular, for each of the 140 realisations, we generated the galaxy distribution in real space (ℛ ) and three different versions of its redshift-space counterpart. In the 𝒱 mocks, we basically only accounted for the spatial derivative of the peculiar velocity of the galaxies along the line of sight. In the 𝒢 catalogues, we added all the integral terms, which are dominated by the weak lensing signal (magnification bias and convergence). Finally, the 𝒪 mocks also include the distortions generated by a non-vanishing peculiar velocity of the observer, which we matched to the observed dipole temperature anisotropy of the cosmic microwave background. The 𝒪 catalogues mimic what will be obtained with Euclid using heliocentric redshifts. The ℛ , 𝒱, and 𝒢 mocks are useful tools for our investigation but do not correspond to anything¹⁵ that can be measured in reality. We note that our study automatically includes wide-angle effects induced by correlating pairs with large separations on the celestial sphere. On the other hand, we do not account for variations in the survey depth on the sky (which can generate a spurious clustering signal if not corrected for in the random catalogue) and redshift errors as we concentrate solely on the signal produced by relativistic RSDs.

We used the likelihood-ratio test to estimate the statistical significance with which the Euclid measurements could reject the null hypothesis that certain relativistic projection effects can be neglected in the theoretical models. In particular, given the summary statistics measured from the individual 𝒪 mocks, we have compared the likelihoods of the models obtained by averaging the clustering signal over all of the 𝒱, 𝒢, and 𝒪 catalogues with no tunable parameters. Our results can be summarised as follows.

(i) Angular power spectrum: The peculiar velocity of the observer noticeably increases the amplitude of the dipole (ℓ = 1), octupole (ℓ = 3), and dotriacontapole (ℓ = 5) signal measured in the heliocentric frame (Fig. 5). Statistically, this boost can be measured with a S/N of 2.1 for the redshift bin that encompasses the entire survey (Table 1) The combined effect of magnification bias and weak lensing convergence slightly increases the clustering signal for ℓ ≳ 50. This contribution is particularly evident in the broad redshift bin 0.9 < z < 1.8, where it is detectable with S/N = 5.4. Slightly smaller statistical significance (up to S/N = 4.5) is obtained by cross-correlating distant narrower bins (see Table 2). We note that the radial projection in the angular power spectrum strongly suppresses the contribution of the source peculiar velocities at small scales (large ℓ). Consequently, the weak-lensing contribution becomes the most dominant correction at these scales, in contrast to its behaviour in the 3D statistics.

(ii) Multipole moments of the 2PCF: The peculiar velocity of the observer introduces small modifications to the clustering signals and the likelihood-ratio test shows that the S/N for rejecting a model that assumes a comoving observer is always smaller than one. The weak lensing contribution amplifies the quadrupole and hexadecapole moments of the 2PCF at large scales (r > 150 h⁻¹ Mpc) for z > 1.3 (see Fig. 7). The S/N of the combined contributions due to lensing and v_o to the multipoles of the 2PCF grows with redshift and reaches 2.5 in the most distant tomographic bin we consider (1.5 < z < 1.8). Our results are in excellent agreement with the Fisher information-based S/N estimates provided by Euclid Collaboration: Jelic-Cizmek et al. (2024). These authors show that neglecting lensing magnification in the models systematically shifts the constraints on the cosmological parameters derived from the Euclid multipoles of the 2PCF. In the ΛCDM framework, the resulting bias is of the order of 0.4 – 0.7 standard deviations. Even larger biases affect model-independent estimates of the growth rate of structure.

(iii) Multipole moments of the power spectrum: The peculiar velocity of the observer imprints rather strong characteristic features on very large scales (k < 0.015 h Mpc⁻¹) that are, however, difficult to detect with the EWSS due to the large statistical uncertainties of the measurements. The S/N is always of order one¹⁶ in all the tomographic bins we consider (Table 4). The weak gravitational lensing signal is mostly noticeable in the tomographic slice 1.5 < z < 1.8, where it can be detected with S/N = 2.4. In terms of the S/N of the different effects, considering the multipole moments of the power spectrum provides very similar results to those obtained from the 2PCF.

Since they are obtained by averaging the measurements of a given summary statistic over 140 mock catalogues, by construction, the models considered above are unbiased and perfectly account for wide-angle effects. In actual surveys, however, observational data are fit to theoretical predictions derived with some flavour of perturbation theory and based on assuming the LPP approximation. We thus repeated the likelihood-ratio test for the multipoles of the power spectrum by contrasting the mock-based predictions with a simple analytical model obtained combining the linear matter power spectrum with Kaiser’s model for RSDs. In order to compare it to the mock data accounting for the Euclid window function, we used the mixing-matrix formalism introduced in Eqs. (46) and (47). We find that the elements of the mixing matrix present high-frequency oscillations (Fig. 12) because the EWSS is made of four disconnected patches that we analyse simultaneously. The window-corrected theoretical model very closely matches the mean from the 𝒱 mocks. The striking agreement shown in Fig. 13 simultaneously demonstrates that (i) our mock catalogues are accurate, (ii) the mixing-matrix approach is also robust for surveys composed of multiple patches, and (iii) the power spectrum multipoles extracted from the 𝒱 mocks are not particularly affected by the wide-angle effects sourced by the velocity-gradient term. On the other hand, in the redshift bin 1.5 < z < 1.8, the window-corrected Kaiser model is rejected with S/N = 2.9 by the 𝒪 mocks that include all relativistic RSDs. This result reinforces the conclusion that the weak lensing term cannot be disregarded at high redshift and on large scales.

Acknowledgments

DB and MYE acknowledge support from the COSMOS network through the ASI (Italian Space Agency) Grants 2016-24-H.0, 2016-24-H.1-2018 and 2020-9-HH.0. MYE and DB acknowledge funding from the Italian Ministry of Education, University and Research (MIUR) through the “Dipartimenti di eccellenza” project “Science of the Universe”. 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 Agenzia Spaziale Italiana, the Austrian Forschungsförderungsgesellschaft funded through BMK, the Belgian Science Policy, the Canadian Euclid Consortium, the Deutsches Zentrum für Luft- und Raumfahrt, the DTU Space and the Niels Bohr Institute in Denmark, the French Centre National d’Etudes Spatiales, the Fundação para a Ciência e a Tecnologia, the Hungarian Academy of Sciences, the Ministerio de Ciencia, Innovación y Universidades, the National Aeronautics and Space Administration, the National Astronomical Observatory of Japan, the Netherlandse Onderzoekschool Voor Astronomie, the Norwegian Space Agency, the Research Council of Finland, 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 (https://www.euclid-ec.org).

References

Appendix A: Validation of the LIGER method

The calculation of the galaxy angular power spectrum with relativistic RSDs (but excluding corrections due to ${\mathit{v}}_{\mathrm{o}}$ and considering constant values of b and 𝒬) has been integrated into the CAMB and CLASS codes (Di Dio et al. 2013). Here, we validate the LIGER method by comparing the mean of the angular power spectra estimated from our ℛ  and 𝒢 mock catalogues to the output of CAMB. For simplicity, we perform the comparison for a full-sky survey and we provide the density-weighted average of the b and 𝒬 functions as input to CAMB. Figure A.1 shows that the two methods agree to better than 1% for ℓ ≲ 90 for the tomographic bin z ∈ (1.1, 1.3). The loss of power of the mock catalogues at higher multipoles is irrelevant for our study (which focuses on larger scales) and can be attributed to the coarse spatial resolution of the input 2LPT simulations (see Sect. 3.2). Similar results are obtained for the other redshift bins.

Fig. A.1. Top: Mean angular power spectra estimated from our mock catalogues (squares with RMS error bars) for the ℛ (red) and 𝒢 (blue) suites are compared to the corresponding output from the CAMB code (crosses). All spectra were binned using intervals Δℓ = 5. Bottom: Relative difference between the spectra (error bars here indicate the standard error of the mean).

All Tables

All Figures

	Fig. 1. Flow chart of the LIGER method in the ‘large-box’ mode (top row) and of the clustering analysis performed in this paper (bottom row). Inputs and outputs are displayed as red and blue parallelograms, respectively, while processes are shown as green rectangles.
In the text

	Fig. 2. Top: Background number density of galaxies in the EWSS as a function of redshift. Bottom: Corresponding evolution, magnification, and linear-bias parameters.
In the text

	Fig. 3. Projected galaxy number counts for one of the 𝒱 mocks in the z ∈ (0.9, 1.1) redshift bin. The black lines enclose the region used in Sect. 8.2.1 as an example of a simply connected domain.
In the text

Fig. 4. Elucidatory example of the statistical test outlined in Sect. 4. In the top panel, the blue histogram on the right-hand side shows the distribution of the likelihood-ratio test statistic λ evaluated from the mock catalogues that do not include all relativistic RSDs (i.e. under the null hypothesis, ℋ0, which we try to rule out using the observed data). The red histogram on the left-hand side, instead, displays the distribution of λ in the mock light cones that account for all effects (i.e. under the alternative hypothesis, ℋ1). The solid and dashed curves represent Gaussian models for the histograms, as is described in the main text. The S/N is a measure of the separation between the two histograms in units of their RMS scatter. The shaded region highlights the realisations in which ℋ0 is ruled out at the 95% confidence level. The bottom panel only differs from the top one in the fact that the covariance matrix C𝒪 has been used to compute all likelihood functions. We use this approximation in the remainder of this paper.

In the text

	Fig. 5. Difference between the angular power spectra of the 𝒪 and 𝒢 mocks for a full-sky survey (red) and the EWSS (blue). The symbols indicate the mean signal, while the error bars show the standard deviation over the 140 realisations.
In the text

	Fig. 6. Mean auto-and cross-angular power spectra extracted from the 𝒪 (red stars) and 𝒱 (green circles) mocks. The shaded region indicates the RMS scatter of the 𝒪 spectra. The panels are ordered as the entries in Tables 1 and 2. The cross-spectra of the tomographic redshift bins are multiplied by ten to improve the readability of the figure.
In the text

	Fig. 7. Mean ℓ = 0, 2, and 4 multipoles of the 2PCF measured from the ℛ (black squares), 𝒱 (green circles), 𝒢 (blue crosses), and 𝒪 (red stars) mocks in the four tomographic redshift bins. The shaded areas highlight the RMS scatter among the 𝒪 light cones.
In the text

	Fig. 8. PDF of the likelihood-ratio test statistic forthe 𝒪–𝒱 test applied to the multipoles of the 2PCF.
In the text

	Fig. 9. Mean ℓ = 0, 2, and 4 multipoles of the power spectrum measured from the 𝒱 (green circles), 𝒢 (blue crosses), and 𝒪 (red stars) mocks in the four tomographic redshift bins. The shaded areas highlight the RMS scatter among the 𝒪 light cones.
In the text

	Fig. 10. Similar to Fig. 8 using the power spectrum multipoles as the data vector.
In the text

	Fig. 11. 2PCF multipoles of the window function for the EWSS (solid) in the redshift bin z ∈ (1.1, 1.3) and for the single patch (dashed) shown in Fig. 3 in the same redshift range. The vertical lines indicate the characteristic length scale V1/3 for the two surveys. The bottom panel offers a more detailed view of the large-scale behaviour of the functions.
In the text

	Fig. 12. Elements of the mixing matrix for the two surveys discussed in Fig. 11 evaluated at k = 0.03 h Mpc−1.
In the text

Fig. 13. Power-spectrum multipoles for the 𝒱 mocks in the redshift range z ∈ (1.1, 1.3). The mean and RMS scatter over the 140 mocks are displayed with squares and error bars, respectively. The solid and dotted lines show the theoretical predictions obtained with the Kaiser model (see Eq. (42)) with and without accounting for the window function of the EWSS (see Eq. (46)), respectively. The vertical dashed line marks the characteristic wavenumber, kC = 2 π/V1/3.

In the text

	Fig. A.1. Top: Mean angular power spectra estimated from our mock catalogues (squares with RMS error bars) for the ℛ (red) and 𝒢 (blue) suites are compared to the corresponding output from the CAMB code (crosses). All spectra were binned using intervals Δℓ = 5. Bottom: Relative difference between the spectra (error bars here indicate the standard error of the mean).
In the text