© The Authors 2023

1. Introduction

Euclid (Laureijs et al. 2011) is an upcoming space mission that aims to explore the nature of dark matter and dark energy through observations of galaxy clustering and weak lensing. Clustering statistics will mainly be extracted from the spectroscopic dataset, which will contain accurate spectroscopic redshifts for tens of millions of galaxies, spanning a wide redshift range (0.9 ≤ z < 1.8) and covering almost a third of the sky. The primary analyses of such data usually focus on the galaxy auto-correlation statistics, in particular the two-point correlation function (2PCF) and its Fourier space counterpart, the power spectrum. These analyses allow for the measurement of the characteristic baryon acoustic oscillation (BAO) feature (e.g., Alam et al. 2017, 2021), together with the signatures of galaxy velocities which give rise to redshift-space distortions (RSD; Kaiser 1987). However, non-Gaussianity introduced by non-linear gravitational evolution means that two-point statistics do not fully characterise all of the information available from galaxy surveys. In addition, models based on perturbation theory start to break down at small scales (at r ≲ 25 h⁻¹ Mpc in the 2PCF, or k ≳ 0.15 h Mpc⁻¹ in the power spectrum); therefore, information from these smaller scales is often excluded from analysis to avoid systematic biases in the modelling (Beutler et al. 2017; Satpathy et al. 2017).

Several alternative clustering summary statistics have been used to extract additional information from the galaxy distribution, such as polyspectra (e.g., Verde et al. 2002; Scoccimarro et al. 2001; Gil-Marín et al. 2016; Philcox & Ivanov 2022) and N-point correlation functions (e.g., Wang et al. 2004; Nichol et al. 2006; Marín et al. 2013; Guo et al. 2015; Slepian et al. 2017). Other proposed statistics include marked correlation functions or power spectra (e.g., White 2016; Satpathy et al. 2019; Massara et al. 2023), counts-in-cells statistics (e.g., Yang & Saslaw 2011; Uhlemann et al. 2017), density-dependent clustering (e.g., Tinker 2007; Chiang et al. 2014; Bayer et al. 2021; Bonnaire et al. 2022; Paillas et al. 2023), anisotropic clustering (e.g., Paz et al. 2011; Kazin et al. 2012; Correa et al. 2019), wavelet-based methods (e.g., Valogiannis & Dvorkin 2022a,b), non-linear transformations of the galaxy field (e.g., Neyrinck et al. 2009; Carron & Szapudi 2014; Wolk et al. 2015), and density field reconstruction (e.g., Wang et al. 2022).

With the advent of large galaxy redshift surveys, cosmic voids – large low-density regions in the galaxy distribution – have emerged as interesting probes of cosmology in many contexts, for example using the void size function (e.g., Pisani et al. 2015; Nadathur 2016; Correa et al. 2019; Contarini et al. 2022), gravitational lensing by voids (e.g., Sánchez et al. 2017; Raghunathan et al. 2020; Bonici et al. 2023), secondary cosmic microwave background (CMB) anisotropies (e.g., Granett et al. 2008; Nadathur & Crittenden 2016; Alonso et al. 2018; Kovács et al. 2019), or void clustering (e.g., Chan et al. 2014; Kitaura et al. 2016; Zhao et al. 2022). Readers can refer to Pisani et al. (2019) for a comprehensive list.

The anisotropic distribution of galaxies around voids is a particularly interesting observable (Lavaux & Wandelt 2012) that has been shown to be a valuable source of cosmological information, providing parameter constraints that are highly complementary to those obtained from galaxy clustering (e.g., Nadathur et al. 2019b; Hamaus et al. 2022; Woodfinden et al. 2022). Individual voids have irregular shapes and arbitrary orientations on the sky, but, given a large enough sample of voids, the assumption of statistical isotropy of the Universe (together with some additional assumptions discussed in more detail in Sect. 2) should give rise to a spherically symmetric distribution of galaxies around voids on average (Ryden 1995). Two effects that spoil this symmetry are distortions due to the Alcock–Paczyński (AP) effect (Alcock & Paczynski 1979), arising from the choice of the fiducial cosmological model used to convert observed galaxy redshifts and angles on the sky into distances, and RSD arising from the peculiar velocities of galaxies around these voids.

As voids are underdense environments, the RSD contributions to the anisotropy can be relatively successfully modelled using linear theory (Hamaus et al. 2014; Cai et al. 2016; Nadathur & Percival 2019; Paillas et al. 2021) without needing to exclude small scales. This provides constraints on the growth rate of structure parametrised by fσ₈ or β (e.g., Paz et al. 2013; Hamaus et al. 2016; Hawken et al. 2017, 2020; Nadathur et al. 2019b; Achitouv 2019; Aubert et al. 2022; Woodfinden et al. 2022). More importantly, Hamaus et al. (2015) and Nadathur et al. (2019b) show that the anisotropies from RSD contributions can be easily distinguished from those arising due to the AP effect. This allows the stacked void profile – equivalent to the void-galaxy cross-correlation function (CCF) – to be used as a ‘standard sphere’ for the AP test, as originally advocated by Lavaux & Wandelt (2012). This observable can then be used to measure the AP parameter, which is the ratio of the transverse comoving distance D_M(z) to the Hubble distance D_H(z) as a function of redshift. The constraints obtained on this quantity by applying this method to the BOSS and eBOSS surveys are a factor of 1.7–3.5 more precise than the constraints from galaxy clustering and the BAO in the same data (Hamaus et al. 2016; Nadathur et al. 2019b, 2020b). This precise AP measurement underpins the power of the void-galaxy CCF as a cosmological probe: Nadathur et al. (2020a) showed that the combination of voids and BAO measurements alone provides greater than 10σ evidence for late-time cosmic acceleration, independent of the cosmic microwave background or supernovae.

In this work we aim to forecast the constraints on the AP parameter D_M/D_H and the growth rate fσ₈ that can be obtained from measurements of the void-galaxy CCF by Euclid. This paper is one of a set of four companion papers published on behalf of the Euclid Consortium, which forecast the expected constraints from cosmic voids using different observables. Contarini et al. (2022) explored the use of the void size function, Bonici et al. (2023) investigated the cross-correlation of voids with weak lensing, while this paper and Hamaus et al. (2022) examine different methods of studying the void-galaxy CCF. The primary difference between our work and that of Hamaus et al. (2022) is that we use a velocity-field reconstruction method to approximately remove large-scale galaxy RSD effects before running the void-finding pipeline (Nadathur et al. 2019a) in order to mitigate against systematic errors due to anisotropic selection effects. On the other hand, instead of applying reconstruction, Hamaus et al. (2022) modify terms in the theoretical model used for the fits, and so measure a very different pattern of anisotropy in the CCF.

The structure of this paper is as follows: in Sect. 2 we provide a theoretical description of the void-galaxy cross-correlation function, as well as the AP test. In Sect. 3 we describe the Flagship galaxy mock and the reconstruction method we apply and which enables us to recover the approximate real-space galaxy positions. We also describe the void finder and detail the properties of the void catalogue. In Sect. 4 we describe our pipeline, obtaining the data vectors and translating this into constraints on parameters of interest. Sect. 5 describes the results we obtain, which are then summarised and put into a larger context in Sect. 6.

2. Theory

2.1. The void-galaxy cross-correlation function in redshift space

The void-galaxy CCF ξ represents the excess probability of finding a galaxy at a given separation from a void centre. We use the notation ξ^rr to denote this quantity when the cross-correlation is performed between the real-space positions of galaxies and the positions of voids identified in this real-space galaxy distribution. For convenience, we refer to ξ^rr as the real-space void-galaxy CCF. For the same fixed void centres, the cross-correlation with the redshift-space galaxy positions is then denoted ξ^rs, and will be referred to as the redshift-space CCF. Here r refers to real space, s to redshift space, the first superscript refers to the voids, and the second to the galaxies. We denote the vectors from the void centre to the galaxy in real and redshift space by r and s respectively, and decompose them into their components perpendicular and parallel to the line-of-sight direction, s = s_⊥ + s_||. Thus s_⊥ = r_⊥ while

$$\begin{array}{c}\hfill {\mathit{s}}_{||}={\mathit{r}}_{||}+\frac{{\mathit{v}}_{||}}{\mathit{aH}},\end{array}$$(1)

where a is the scale factor, H the Hubble rate, v_|| is the component of the galaxy peculiar velocity along the line of sight (not the pairwise velocity; see e.g., Massara et al. 2022), and the line of sight is directed through the void centre.

Since, by construction, the number of void-galaxy pairs is conserved under the mapping between real and redshift space, ξ^rs and ξ^rr are related by a convolution with the (position-dependent) velocity distribution:

$$\begin{array}{c}\hfill 1+{\xi }^{\mathrm{rs}}(\mathit{s})={\displaystyle \int }[1+{\xi }^{\mathrm{rr}}(\mathit{r})]P(v_{||},\mathit{r})\mathrm{d}{v}_{||}.\end{array}$$(2)

This is recognisable as the streaming model (e.g., Peebles 1979), which was used to describe the void-galaxy CCF by Paz et al. (2013) and also, for example, by Cai et al. (2016) and Paillas et al. (2021). Although – for reasons explained below – we have chosen to hold fixed the void positions identified in the real-space galaxy field, Eq. (2) would be equally valid if both the CCFs were defined with respect to the void positions in the redshift-space galaxy field, that is if ξ^rs and ξ^rr were replaced by ξ^ss and ξ^sr respectively. Indeed this is the approach taken in a number of other papers (Paz et al. 2013; Hamaus et al. 2016; Achitouv 2017). However, the number of voids in a size-selected sample is not conserved under the RSD mapping (Chuang et al. 2017; Nadathur et al. 2019a; Correa et al. 2022), so Eq. (2) cannot be used to relate CCFs obtained from running void-finding independently in the two spaces, for example ξ^ss and ξ^rr.

We can view the galaxy peculiar velocities as being composed of a term describing a coherent velocity outflow from the void centre, plus a random velocity component along a random direction. Given the assumption of statistical isotropy in real space, when averaged over sufficiently large numbers of voids, the coherent outflow $\mathit{v}(r)=v_{\mathrm{r}}(r)\widehat{\mathit{r}}$ must be directed radially outwards from the void centre and be spherically symmetric about it. Thus the line-of-sight component of the galaxy velocity can be written as

$$\begin{array}{c}\hfill v_{||}(r,{\mu }_{\mathrm{r}})=v_{\mathrm{r}}(r){\mu }_{\mathrm{r}}+{\stackrel{\sim }{v}}_{||},\end{array}$$(3)

where μ_r ≡ cos(θ) is the cosine of the angle θ between r and the line of sight, and ${\stackrel{\sim }{v}}_{||}$ is a zero-mean random variable describing the projection of the stochastic velocity component along the line of sight.

Using Eq. (3) to change variables, Woodfinden et al. (2022) show that Eq. (2) is equivalent to

$$\begin{array}{c}\hfill 1+{\xi }^{\mathrm{rs}}(\mathit{s})={\displaystyle \int }[1+{\xi }^{\mathrm{rr}}(\mathit{r})]J_{\mathit{r},\mathit{s}}P({\stackrel{\sim }{v}}_{||},\mathit{r})\mathrm{d}{\stackrel{\sim }{v}}_{||},\end{array}$$(4)

where

$$\begin{array}{c}\hfill J_{\mathit{r},\mathit{s}}\equiv \left|\frac{\mathrm{d}\mathit{r}}{\mathrm{d}\mathit{s}}\right|={[1+\frac{v_{\mathrm{r}}}{raH}+\frac{r{v}_{\mathrm{r}}^{\prime }-v_{\mathrm{r}}}{raH}{\mu }_{\mathrm{r}}^2]}^{-1}\end{array}$$(5)

is the Jacobian of the transformation from real to redshift space, prime denotes derivatives with respect to r, and $P({\stackrel{\sim }{v}}_{||},\mathit{r})$ now describes the distribution of ${\stackrel{\sim }{v}}_{||}$ about zero. Equation (4) is the same model described by Nadathur & Percival (2019). Several other works (e.g., Cai et al. 2016; Hawken et al. 2020; Aubert et al. 2022) use approximations to this expression, first by assuming a Dirac delta function form for the distribution $P({\stackrel{\sim }{v}}_{||},\mathit{r})$, which reduces the integral equation to

$$\begin{array}{c}\hfill 1+{\xi }^{\mathrm{rs}}(\mathit{s})=[1+{\xi }^{\mathrm{rr}}(\mathit{r})]J_{\mathbf{r},\mathbf{s}},\end{array}$$(6)

and then approximating the square brackets in Eq. (5) by a series expansion, which one can choose to truncate at some desired order. Hamaus et al. (2020) introduce additional nuisance parameters to Eq. (6), and Hamaus et al. (2022) further modify some terms of the Jacobian expansion.

A key assumption for all models based on Eq. (4) is that the coherent velocity term v_r(r) is spherically symmetric and radially directed, which depends crucially on the assumption that in real space there is no preferred orientation direction for the large number of individually asymmetric voids. This follows from statistical isotropy only provided there is no orientation-dependent selection effect introduced in creating the void sample. Nadathur et al. (2019a) argued that exactly such an effect is introduced for voids identified in redshift space, since regions with velocity outflows oriented along the line of sight direction appear to have lower densities in redshift space. This means that elongated voids aligned along the line-of-sight direction are more likely to be selected as void candidates if void-finding is performed on the redshift-space galaxy positions.

This orientation-dependent bias has a large effect on the measured CCFs (Nadathur et al. 2019a; Correa et al. 2022; Paillas et al. 2023). The CCF with real-space galaxy positions, ξ^sr, becomes strongly anisotropic, the coherent velocity outflow v_r no longer has spherical symmetry, and the multipoles of the CCF with redshift-space galaxy positions, ξ^ss, have a very different shape. The loss of spherical symmetry in particular introduces significant challenges for the modelling approach described above. To avoid this problem, in this paper we do not consider the case of voids identified in the redshift-space galaxy field. As the real-space galaxy positions are in general not known, we used the reconstruction technique described in Sect. 3.2 below to approximately recover them before performing void-finding, and used these void positions to approximately measure ξ^rs and ξ^rr.

2.2. Modelling the velocity distribution

Given the stated assumptions, the convolution described by Eq. (4) – or, equivalently, Eq. (2) – is completely general and exact. In order to make specific predictions, however, we need to specify a model for the velocity field around voids to obtain the radial term v_r(r) and the distribution $P({\stackrel{\sim }{v}}_{||},\mathit{r})$.

The radial outflow velocity v_r itself is commonly modelled (e.g., Cai et al. 2016; Nadathur et al. 2019b; Woodfinden et al. 2022) using a a linearised form of the continuity equation (Peebles 1980):

$$\begin{array}{c}\hfill v_{\mathrm{r}}(r)=-\frac{1}{3}faHr\mathrm{\Delta },\end{array}$$(7)

where

$$\begin{array}{c}\hfill \mathrm{\Delta }(r)=\frac{3}{r^3}{\displaystyle {\int }_0^r}\delta (y)y^2\mathrm{d}y\end{array}$$(8)

is the fully non-linear enclosed density profile of dark matter, $\delta (r)=\rho /\overline{\rho }-1$ is the matter overdensity, and f is the linear growth rate. Nadathur et al. (2019b) found that in simulated galaxy catalogues the dark matter profile of voids approximately scales as Δ ∝ σ₈, where σ₈ describes the amplitude of the matter density perturbations on 8 h⁻¹ Mpc scales¹. When combined with Eq. (7), this implies a scaling v_r ∝ fσ₈.

Equation (7) has been found to give a reasonably good description of the matter velocity field around voids by some studies (Hamaus et al. 2014; Nadathur & Percival 2019). However, other works (Achitouv 2017; Paillas et al. 2021) have found that it can be necessary to add correction terms. Moreover, it has been shown by Massara et al. (2022) that for some void finders, when void-finding is performed on sparse tracers of the matter density field such as galaxies, an effective velocity bias can arise between the matter velocity which satisfies Eq. (7), and the galaxy velocity appearing in Eq. (5). This can mean that even if Eq. (7) accurately describes the velocity field for matter, galaxy velocities can differ from this in the inner regions of voids. These considerations will need to be taken into account when improving the velocity modelling beyond that in Eq. (7). It is not our aim to attempt this task in this paper, but only to forecast the constraints from the void-galaxy cross-correlation function that could be obtained from such a model in the future. We therefore applied a template-based method that assumes only that the scaling inferred from Eq. (7) holds. That is, we measured the mean velocity profile $v_{\mathrm{r}}^{\mathrm{sim}}(r)$ from the Flagship simulation as described in Sect. 4.2 below, and assume that the velocity in other cosmological models is given by

$$\begin{array}{c}\hfill \frac{v_{\mathrm{r}}(r)}{aH}=\frac{f{\sigma }_8}{f^{\mathrm{sim}}{\sigma }_8^{\mathrm{sim}}}{\left(\frac{v_{\mathrm{r}}(r)}{aH}\right)}^{\mathrm{sim}}.\end{array}$$(9)

For the velocity distribution $P({\stackrel{\sim }{v}}_{||},\mathit{r})$, we assume a Gaussian form,

$$\begin{array}{c}\hfill P({\stackrel{\sim }{v}}_{||},\mathit{r})=\frac{1}{\sqrt{2\pi {\sigma }_{v_{||}}^2(r)}}exp(-\frac{{\stackrel{\sim }{v}}_{||}^2}{2{\sigma }_{v_{||}}^2(r)}),\end{array}$$(10)

which has been shown to be a good approximation by Nadathur & Percival (2019) and Paillas et al. (2021). As with v_r(r), we assume the width of the dispersion σ_v||(r) is spherically symmetric about the void centre and take a template-based approach to its modelling. We measured a template dispersion profile ${\sigma }_{v_{\left|\right|}}^{\text{sim}}(r)$ from the Flagship simulation (Sect. 4.2). At large distances from the void centre, this profile asymptotes to a constant, ${\sigma }_{\mathrm{v}}^{\mathrm{sim}}\equiv {\sigma }_{v_{||}}{}^{\text{sim}}(r\to \infty )$ (see also Hamaus et al. 2015). To account for changes in other cosmologies, we renormalised this function by an amplitude factor σ_v, taken to be a free parameter to be marginalised over:

$$\begin{array}{c}\hfill {\sigma }_{v_{||}}(r)=\frac{{\sigma }_{\mathrm{v}}}{{\sigma }_{\mathrm{v}}^{\mathrm{sim}}}{\sigma }_{v_{||}}^{\mathrm{sim}}(r).\end{array}$$(11)

In summary, the template-based approach to modelling the velocity field around voids requires the specification of two template functions, $v_{\mathrm{r}}^{\mathrm{sim}}(r)$ and ${\sigma }_{v_{\left|\right|}}^{\text{sim}}(r)$, and depends on two free parameters, fσ₈ and σ_v.

2.3. The real-space correlation function

Although the velocity model has been specified as in Sect. 2.2, the RSD mapping of the void-galaxy correlation described by Eq. (4) only relates ξ^rs to its counterpart in real space, ξ^rr. Though some empirical fitting formulae have been used (Paz et al. 2013; Hamaus et al. 2014; Hawken et al. 2017), models which attempt to predict ξ^rr directly from the cosmological parameters do not exist at present². Previous studies have therefore taken different approaches to obtain the real-space CCF. Hamaus et al. (2020, 2022) apply an inversion algorithm (Pisani et al. 2014) to recover it from the redshift-space CCF projected along the line of sight (implictly assuming that the real-space correlation function is spherically symmetric about the void centre). Nadathur et al. (2019b, 2020b) and Woodfinden et al. (2022) instead measure the average ξ^rr from a large number of mocks that match the galaxy clustering seen in the data sample: this is effectively a template-based method similar to that described for v_r(r) and σ_v||(r) above.

In this paper, we use the reconstruction-based estimators described in Sect. 3.2 to estimate the real-space CCF ξ^rr(r) from the simulation data. While we have argued that in the absence of orientation-dependent selection effects, statistical isotropy should imply spherical symmetry, ξ^rr(r) = ξ^rr(r), when measured in the true background cosmology, this will in general not be true in the presence of the AP distortions discussed in Sect. 2.4, and could also be spoiled by imperfections in the reconstruction-based RSD removal step. Our implementation of Eq. (4) therefore uses the full ξ^rr(r) without imposing an assumption of isotropy.

Our treatment of ξ^rr here differs from that in Nadathur et al. (2019b) and Woodfinden et al. (2022) because, unlike in those works, we currently only have access to the single representative Flagship simulation here, so we cannot construct a template based on the mean of thousands of mocks (this situation will change by the time of the final Euclid data analyses). This means that our estimates of ξ^rr are necessarily much noisier. As it is measured from the same data as the redshift-space correlation function ξ^rs we wish to fit, correlations between the theory model and data vector may also be introduced, and we consider the implications of this in Appendix A.

2.4. Alcock–Paczynski effect

Besides RSD, another source of anisotropies in the observed void-galaxy cross-correlation function is the AP effect. This arises from the conversion of observed galaxy angular positions and redshifts to distances in order to apply the modelling developed above. This conversion requires the assumption of a fiducial cosmological model which may differ from the true background and thus lead to incorrectly inferred distances, introducing anisotropies.

To parameterise this effect, we introduce two AP scaling parameters describing the ratio of distances parallel and perpendicular to the line of sight to those in the fiducial model,

$$\begin{array}{c}\hfill {\alpha }_{\Vert }=\frac{D_{\mathrm{H}}(z)}{D_{\mathrm{H}}^{\mathrm{fid}}(z)},{\alpha }_{\perp }=\frac{D_{\mathrm{M}}(z)}{D_{\mathrm{M}}^{\mathrm{fid}}(z)},\end{array}$$(12)

where D_M(z) is the comoving angular diameter distance and D_H(z) = c/H(z) is the Hubble distance at redshift z. In this model the correlation function scales as

$$\begin{array}{c}\hfill {\xi }^{\mathrm{rs}}(s_{\perp },s_{\Vert })={\xi }^{\mathrm{rs},\mathrm{fid}\left\}\right\}}({\alpha }_{\perp }{s}_{\perp }^{\mathrm{fid}},{\alpha }_{\Vert }{s}_{\Vert }^{\mathrm{fid}}),\end{array}$$(13)

where the superscript ^fid denotes quantities calculated in the fiducial cosmology.

The parameters α_∥ and α_⊥ can also be recast into the following combinations:

$$\begin{array}{c}\hfill \alpha ={\alpha }_{\perp }^{2/3}{\alpha }_{\Vert }^{1/3},\end{array}$$(14)

which describes an isotropic volume dilation, and

$$\begin{array}{c}\hfill ϵ=\frac{{\alpha }_{\perp }}{{\alpha }_{\Vert }},\end{array}$$(15)

describing anisotropic distortions. The dilation parameter α can be measured using a standard ruler, such as the BAO scale (e.g., Blake & Glazebrook 2003; Seo & Eisenstein 2003), by comparing a measured distance scale to one predicted by fundamental physics models. In our analysis, however, the calibrated template functions used depend on the void size cut applied to the sample. This cut is numerically identical to the one applied to the data. However, in general if the template cosmology differs from the truth, due to AP dilation the same numerical cut might not select the same population of voids, leading to a shift in the mean void size. Therefore, we conservatively choose to avoid using the observed void size as an absolute ruler by marginalising over values of α in order to reduce the possible effects of the template cosmology. To achieve this we isotropically rescaled all of the template profiles ξ^rr(r), v_r(r) and σ_v||(r) with α, as described by Nadathur et al. (2019b). This rescaling removes all sensitivity of the model to α. The model then depends only on the AP distortion parameter ϵ, via its effect on the quadrupole and (to a lesser extent) hexadecapole moments.

2.5. Multipole decomposition

The model described above applies generally to the full three-dimensional correlation function, which may be written as ξ^rs(s, μ_s), where s = |s| and μ_s ≡ s_||/s. In practice it is preferable to compress the information into a small number of Legendre multipoles, as

$$\begin{array}{c}\hfill {\xi }_{\ell }^{\mathrm{rs}}(s)=\frac{2\ell +1}{2}{\displaystyle {\int }_{-1}^1}{P}_{\ell }({\mu }_{\mathrm{s}}){\xi }^{\mathrm{rs}}(s,{\mu }_{\mathrm{s}})\mathrm{d}{\mu }_{\mathrm{s}},\end{array}$$(16)

where P_ℓ(μ) are the Legendre polynomials of order ℓ. Symmetry dictates that only the even multipoles are non-zero, and we restrict our attention to the monopole, quadrupole and hexadecapole moments (ℓ = 0, 2, 4) only, as these contain almost all of the measurable information in the correlation.

The real-space cross-correlation function ξ^rr can also be decomposed into multipoles analogously to Eq. (16). As mentioned in Sect. 2.3, a commonly used approximation, justified by the assumption of spherical symmetry in real space, is that only the monopole term ${\xi }_0^{\text{rr}}$ is non-zero. However, this assumption can be spoiled by the AP effect or imperfect reconstruction of real-space galaxy positions, so for generality we keep track of the high-order multipoles of ξ^rr as well.

3. Data

3.1. Galaxy catalogue

We used the Euclid Flagship galaxy mock (Castander et al., in prep.), which is constructed from an N-body simulation with 2 trillion particles of mass m_p ≈ 2×10⁹ h⁻¹ M_⊙. The simulation was run using PKDGRAV3 (Potter et al. 2017) in a box of side L = 3780 h⁻¹ Mpc, with a flat Lambda Cold Dark Matter (ΛCDM) cosmology and matter density Ω_m = 0.319, baryon density Ω_b = 0.049, dark energy density Ω_Λ = 0.681, scalar spectral index n_s = 0.96, Hubble parameter h = H₀/(100 km s⁻¹ Mpc⁻¹) = 0.67, and the RMS value of density fluctuations on 8 h⁻¹ Mpc scales σ₈ = 0.83. These values were chosen to be very close to those obtained by Planck Collaboration VI (2020). The halo catalogue was constructed on the fly, with halos identified using Rockstar (Behroozi et al. 2013), and then populated with a halo occupation distribution (HOD) model which had been calibrated to reproduce observables such as galaxy luminosity and clustering statistics as a function of galaxy luminosity and colour. This was used to create a lightcone mock covering an octant (5157 deg²) of the sky, which was then cut to match the expected observed Hα flux f_Hα > 2 × 10⁻¹⁶ erg s⁻¹ cm⁻², and expected redshift range 0.9 ≤ z < 1.8 for the Euclid spectroscopic galaxy sample (Laureijs et al. 2011; Costille et al. 2018). To simulate the expected sample completeness of Euclid, we randomly downsampled the final catalogue, retaining 60% of all objects.

For each mock galaxy in this catalogue we record both the true redshift (corresponding to the background Hubble expansion) and the observed redshift, which includes the Doppler effect due to the galaxy peculiar velocities. We also consider the effects of redshift errors by adding an additional random component drawn from a Gaussian with RMS σ_z = 0.001³. The final resulting galaxy catalogue contains about 6.5×10⁶ mostly central galaxies. For the analysis described below we divided this sample into four non-overlapping redshift bins as illustrated in Fig. 1, chosen to match those used in Euclid forecasts by Euclid Collaboration (2020). The effective redshift of each bin is generically defined as the weighted average redshifts of void-galaxy pairs in the bin,

$$\begin{array}{c}\hfill z_{\mathrm{eff}}=\frac{{\sum }_{\mathit{ij}}\left(\frac{Z_i+z_j}{2}\right)W_i{w}_j}{{\sum }_{\mathit{ij}}{W}_i{w}_j},\end{array}$$(17)

Fig. 1. Left: distribution of void sizes in each of the four redshift bins of the Flagship mock catalogue. The vertical bars indicate the void size cuts applied to each catalogue, and dashed lines represent the voids which are discarded. Right: comoving number density n(z) of galaxies and voids as function of redshift in Flagship. For voids, n(z) is multiplied by 100 for visibility. The edges of the redshift bins are shown by dashed vertical lines. The decrease in the number of voids near the edges of each bin arises because void-finding is performed on each bin individually as discussed in the text; some void candidates near the redshift edges of the bin are then missed.

where Z_i is the redshift of the ith void centre, z_j is redshift of the jth galaxy, and W_i and w_j are optional weights associated with individual voids and galaxies when computing the CCF (as the Flagship mock does not include observational systematics, we do not use weights, that is they are all set to unity). The effective redshifts of the different bins are z_eff = 1.0, 1.2, 1.4 and 1.64.

3.2. Reconstruction

As discussed in Sect. 2, in order to avoid an orientation-dependent selection bias in the void sample that would invalidate the theory model used, void-finding must be performed on the real-space galaxy distribution without large-scale RSD effects. To achieve this, we used a velocity-field reconstruction method, related to the reconstruction commonly used in BAO analysis (Eisenstein et al. 2007; Padmanabhan et al. 2012), to approximately remove galaxy RSD before the void-finding step (Nadathur et al. 2019a).

Reconstruction was performed using an iterative fast Fourier transform (FFT) algorithm (Burden et al. 2015) which solves the Zeldovich equation in redshift space

$$\begin{array}{c}\hfill \mathrm{\nabla }·\mathbf{\Psi }+\frac{f}{b}\mathrm{\nabla }·[(\mathbf{\Psi }·\widehat{\mathit{r}})\widehat{\mathit{r}}]=-\frac{{\delta }_{\mathrm{g}}}{b},\end{array}$$(18)

for a Lagrangian displacement field Ψ, where f is the growth rate, b is the linear galaxy bias, and δ_g is the galaxy overdensity in redshift space. This step is implemented in the pyrecon package⁴. It is performed on a regular grid, and densities estimated on the grid are first smoothed with a Gaussian kernel of width R_s before solving for the displacement. Given the solution to Eq. (18), individual galaxy positions are shifted by −Ψ_RSD evaluated at their locations, where

$$\begin{array}{c}\hfill {\mathbf{\Psi }}_{\mathrm{RSD}}=-f(\mathbf{\Psi }·\widehat{\mathit{r}})\widehat{\mathit{r}}.\end{array}$$(19)

This shift subtracts the estimated RSD contribution to the observed galaxy redshifts, in order to approximately recover the real-space galaxy field. The algorithm for solving Eq. (18) is the same as that used for BAO reconstruction (e.g., in the eBOSS analyses Gil-Marín et al. 2020; Bautista et al. 2021) with the difference that the shifts applied to the galaxy positions only attempt to correct for large-scale RSD and not for non-linear evolution.

If we define a new variable ψ ≡ bΨ, Eq. (18) can be rewritten as $\mathrm{\nabla }·\mathit{\psi }+\beta \mathrm{\nabla }·[(\mathit{\psi }·\widehat{\mathit{r}})\widehat{\mathit{r}}]=-{\delta }_{\mathrm{g}}$, where β = f/b is the RSD parameter. The shifts applied to the galaxy positions in order to remove RSD are then ${\mathbf{\Psi }}_{\mathrm{RSD}}=-\beta (\mathit{\psi }·\widehat{\mathit{r}})\widehat{\mathit{r}}$, and it is clear that the output of our reconstruction procedure depends only on the parameter β and not on f and b individually⁵. The fiducial growth rate for the Flagship cosmology can be calculated as $f(z)={\mathrm{\Omega }}_{\mathrm{m}}^{\gamma }(z)$ (Linder 2005), where Ω_m is the matter density parameter, and γ = 6/11 ≈ 0.55. To determine the fiducial linear bias b, we measured the galaxy power spectrum in real-space P_g(k), calculated the theoretical matter power spectrum P_m(k) for Flagship using CAMB (Lewis et al. 2000), and fit for b² = P_g(k)/P_m(k), using only large-scale modes at k < 0.07 h Mpc⁻¹.

The fiducial values b^fid and β^fid obtained in this way are reported in Table 1 for each redshift bin. However, the growth rate f is a parameter in any fit we perform, and so we do not want to fix β to this fiducial value, but instead to marginalise over it. To achieve this in a computationally feasible way, we pre-computed the reconstruction on a grid of β values in the range [β^fid−0.15,β^fid+0.15] around the fiducial value. The data vectors obtained from each of these approximations to the real-space galaxy field consist of the monopole, quadrupole and hexadecapole (all ranging from 0 to 120 h⁻¹ Mpc), and are interpolated over when varying β in the analysis.

As the Zeldovich approximation breaks down on small scales, reconstruction also depends on the width of the Gaussian smoothing kernel, R_s, which dictates the smallest scales used. Reconstruction with a smoothing radius that is too small will over-correct the RSD, due to contributions from non-linear small-scale fluctuations, while a smoothing radius that is too big will fail to remove all of them. To determine the optimal smoothing radius to use, we computed the quadrupole of the post-reconstruction galaxy power spectrum, P₂(k), to use as a proxy. In real space this quantity is consistent with zero, while RSD introduce significant anisotropies in redshift space (Fig. 2). In each redshift bin we tested reconstruction with multiple smoothing radii at the fiducial value β^fid, and pick the one that results in a quadrupole which is roughly consistent with zero. This gives R_s = 8 h⁻¹ Mpc for the first two redshift bins, and R_s = 9 h⁻¹ Mpc for the last two.

Fig. 2. Quadrupole moments of the galaxy-galaxy power spectrum in each redshift bin of the Flagship mock, measured in redshift space (yellow), real space (blue) and after applying reconstruction to remove RSD from the galaxy positions (red). Reconstruction was performed using the fiducial values of βfid listed in Table 1 and using smoothing scales Rs of 8 and 9 h−1Mpc for the first two and last two bins, respectively. The similarity of the red and blue lines and their proximity to zero is because reconstruction is successfully removing the large-scale RSD-induced anisotropy in the galaxy positions.

The reconstruction method outlined here cannot perfectly remove all redshift-space effects in the galaxy distribution. We are interested in quantifying the impact of these inaccuracies, as well as those caused by redshift errors, on the void-finding and data analysis pipeline below. We therefore repeated our analysis pipeline in three different scenarios, labelled as: perfect reconstruction, using the true real-space galaxy positions from the simulation; realistic reconstruction, applying reconstruction to the observed redshifts to obtain approximate real-space galaxy positions; and realistic reconstruction with σ_z, applying reconstruction to observed redshifts including redshift errors.

3.3. Void finding

Voids were found in the real space or post-reconstruction galaxy distribution using the voxel algorithm, implemented in REVOLVER⁶ (Nadathur et al. 2019c). voxel is a watershed-based void-finding algorithm that operates on a particle-mesh interpolation of the galaxy density field. In the particle-mesh assignment step, tracers (in our case galaxies) are first placed onto a grid to estimate the local density field, a step which is much faster than commonly used tessellation algorithms, and scales linearly with the number of tracers, 𝒪(N) (as opposed to 𝒪(N^3/2) for tesselation). The mesh size is determined by requiring each cubic grid cell, or voxel, to have a side length $a_{\mathrm{vox}}=0.5{\left(\frac{4\pi \overline{n}}{3}\right)}^{-1/3}$, where $\overline{n}$ is the approximate value for the mean number density of galaxies. As $\overline{n}$ varies by a factor of about 5 over the Flagship redshift range, in order to maintain an appropriate resolution in the highest density bin while maintaining a reasonable number of voxels for computational efficiency, we decided to perform the void-finding separately in each of the four redshift bins, using a different a_vox for each. The voxel algorithm (like all void-finder algorithms) unavoidably misses some voids near the redshift edges of the galaxy sample, so this analysis choice leads to the drops in void numbers around redshifts z = 1.1, 1.3 and 1.5 seen in Fig. 1. An alternative approach could involve the use of overlapping redshift bins, which would require a more delicate treatment of correlation between redshifts. We have found that this does not improve our final constraints enough to justify the complexity of the analysis.

In order to correct for the survey window and redshift-dependent or angular selection effects, this density estimate is normalised by the density of a random unclustered catalogue of points which has the same window and selection effects imposed, projected on the same grid. We see from this that, in addition to being significantly faster, this method of estimating the galaxy density field can then more easily account for complex survey masks and systematic effects than popular tessellation-based density estimation methods used in other void-finding codes, such as those based on ZOBOV (Neyrinck 2008). For our purposes with the Flagship catalogue, we created a random catalogue with 50 times as many points as the galaxy catalogue, covering the Flagship octant and following the same redshift distribution n(z).

The estimated galaxy density field is then smoothed with a Gaussian filter of smoothing length $r_{\mathrm{s}}={\overline{n}}^{-1/3}$, in order to reduce shot-noise fluctuations on smaller scales. Voids are then identified as the sites of local minima in the smoothed density field, and their extents are determined using a watershed method, similar to that of ZOBOV, to create a final catalogue of non-overlapping voids. Each void is composed of a contiguous set of voxels, which define its volume. Individual voids have irregular shapes and orientations, but we define an effective void radius as the radius of a sphere with the same volume,

$$\begin{array}{c}\hfill R_{\mathrm{v}}={\left(\frac{3}{4\pi }{N}_{\mathrm{vox}}{V}_{\mathrm{vox}}\right)}^{\frac{1}{3}},\end{array}$$(20)

where N_vox is the number of voxels in the void and $V_{\text{vox}}=a_{\text{vox}}^3$ is the voxel volume. The centre of each void is identified as the location of the lowest-density voxel within it, being the position of the density minimum. This results in voids with very few galaxies in the centre.

For each reconstruction scenario considered, the full void catalogues obtained in this way contain about 87 000 voids (with small variations in number) across all redshift bins. However, we expect that smaller voids will not be adequately described by our approach to modelling the coherent velocity outflow v_r(r). This is due to a combination of factors. The dynamics of small voids are more likely to be dominated by the presence of neighbouring large density fluctuations so that the isolated void model is inappropriate (see also Schuster et al. 2023); they are more likely to be segments of larger voids that appear artificially small because they extend outside of the survey boundaries; and they are also more likely to be entirely spurious detections arising from shot-noise fluctuations (Cousinou et al. 2019). Therefore, we applied a void size cut to exclude them from the catalogue used for clustering measurements⁷:

$$\begin{array}{c}\hfill R_{\mathrm{v}}>R_{\mathrm{cut}}=1.5{\overline{n}}^{-\frac{1}{3}}.\end{array}$$(21)

Keeping only voids that survive this cut reduces the total numbers by about 40%. The radius cuts and the final sizes of the void catalogues are presented in Table 1.

4. Method

4.1. Measuring the void-galaxy cross-correlation function

We measured correlation functions using the Python wrapper pycorr⁸ for Corrfunc (Sinha & Garrison 2020; Sinha et al. 2019). Corrfunc is a pair-counting engine that counts the number of pairs of voids and galaxies within a given distance and angular separation bin. This was then compared to the equivalent pair counts computed using random unclustered catalogues to obtain an estimate of the CCF using the Landy-Szalay estimator (Landy & Szalay 1993):

$$\begin{array}{c}\hfill \xi (r,\mu )=\frac{{\mathrm{DD}}_{12}(r,\mu )-{\mathrm{DR}}_{12}(r,\mu )-{\mathrm{DR}}_{21}(r,\mu )+{\mathrm{RR}}_{12}(r,\mu )}{{\mathrm{RR}}_{12}(r,\mu )}.\end{array}$$(22)

Here D and R refer to the data and unclustered random catalogues respectively, subscripts 1 and 2 are used to differentiate between voids and galaxies, and each term XY refers to the number of pairs of species X and Y in the (r, μ) bin, normalised relative to the total possible number of such pairs. Thus DD₁₂ refers to normalised number of void-galaxy pairs, DR₁₂ to the normalised number of pairs of voids and galaxy randoms, and so on. We created the unclustered random catalogues as described in Sect. 3.3, such that they have 50 times as many points as the respective data samples, are uniformly sampled from the Flagship octant footprint and follow the redshift distributions of voids and galaxies respectively (Fig. 1). We measured pair counts and correlation functions 30 radial bins over 0 < r < 120 h⁻¹ Mpc and in 200 angular bins over −1 < μ < 1, keeping distances in units of h⁻¹ Mpc rather than rescaling them in units of the void size, as is done in some void publications.

It is important to note that shot-noise fluctuations in the galaxy random catalogue used in the mesh density estimation step during void finding (Sect. 3.3) will be slightly correlated with the positions of identified voids. We therefore used an independent realisation of the galaxy random catalogue (generated using the same method but with a different random seed) for measuring correlations via Eq. (22). If instead the same random catalogue were used for both, small artefacts could be propagated to the correlation ξ in the regions very close to the void centres.

We used the estimator in Eq. (22) with the pair separations determined by the positions of voids and galaxies in different combinations of real and redshift space to measure the different correlation functions relevant to our analysis. For instance, in the idealised case of perfect reconstruction, voids are identified in the true real-space galaxy field, and their cross-correlation with the real and redshift-space galaxy positions give ξ^rr and ξ^rs, respectively. On the other hand, when we use realistic reconstruction we use the superscript p to distinguish positions in this post-reconstruction approximation of the real-space galaxy field from the true real-space positions. Voids are found in the post-reconstruction field and, as before, their positions are held fixed to allow measurements of ξ^ps and ξ^pp with the galaxies in redshift space and approximate real space respectively.

Figure 3 shows the recovered monopole, quadrupole and hexadecapole moments for the true real-space void-galaxy CCF ξ^rr in the first Flagship redshift bin, and its comparison with the equivalent multipoles of ξ^pp in the realistic reconstruction scenarios with and without including redshift errors. Reconstruction recovers the real-space CCF multipoles well, and anisotropies in all cases are negligible, as evidenced by the quadrupole and hexadecapole moments being close to zero. This demonstrates that our void sample selection successfully avoids the problem of orientation-dependent selection bias discussed in Sect. 2.

Fig. 3. Void-galaxy CCF in real space, measured in the Flagship simulation zeff = 1.0 redshift bin. From top to bottom, the panels show the monopole, quadrupole and hexadecapole moments of the CCF. The multipoles of ξrr, corresponding to the perfect reconstruction scenario and measured using the true real-space positions are shown in blue, while light and dark red respectively show the multipoles of ξpp obtained using reconstruction without and with redshift errors. The inset panels show the residuals with respect to the perfect reconstruction case of ${\xi }_{\mathcal{l}}^{\text{rr}}$ for each multipole. In all cases ξ2 and ξ4 are close to zero, indicating that reconstruction successfully removes orientation-dependent selection bias.

4.2. Constructing template functions

As discussed in Sect. 2, we adopt a template-fitting approach to the analysis, and so we require templates for the three functions v_r(r), σ_v||(r) and ξ^rr(r). These were measured from the Flagship simulation data as follows.

To measure the radial velocity profile v_r(r), we measured the mean radial component of galaxy peculiar velocities (given in the simulation box rest frame) in spherical shells around the void centre, where the average was taken over all void-galaxy pairs in that separation bin. We used 30 linearly spaced radial bins, up to r = 120 h⁻¹ Mpc. The mean radial velocity profiles at different redshifts are shown in the top panel of Fig. 4. The velocity is positive over almost the entire radial range, corresponding to an outflow from the void centre as expected. The profile shows a peak in the vicinity of the edge of the void and at larger radii it goes to zero, showing that there is no net bulk motion far outside of the void, as expected.

Fig. 4. Measured radial outflow velocity profile vr(r) (top), and normalised line-of-sight velocity dispersion σv||(r) (bottom) of galaxies around voids found in the four redshift bins of Flagship. Shading represents the 68.3% (1σ) measurement uncertainty coming from Poisson noise. The dispersion amplitude σv for the four bins can be found in Table 1.

We also measured the width of the line-of-sight velocity PDF, σ_v||(r), from the distribution of the galaxy peculiar velocities in the same spherical shells. The full vector form of Eq. (3) can be written as

$$\begin{array}{c}\hfill \mathit{v}(\mathit{r})=v_{\mathrm{r}}(r)\widehat{\mathit{r}}+\stackrel{\sim }{v}\widehat{\mathit{n}},\end{array}$$(23)

where the stochastic velocity component has random magnitude $\stackrel{\sim }{v}$ and is directed along a random direction $\widehat{\mathit{n}}$. The quantity we wish to estimate is ${{\sigma }_{v_{||}}}^2(r)=⟨{\stackrel{\sim }{v}}^2{(\widehat{\mathit{n}}·\widehat{\mathit{r}})}^2⟩$ and from Eq. (23) we obtain (Fiorini et al. 2022)

$$\begin{array}{c}\hfill {\sigma }_{v_{||}}(r)=\frac{\sqrt{{⟨|\mathbf{v}|}^2⟩-v_{\mathrm{r}}^2}}{\sqrt{3}},\end{array}$$(24)

where ⟨|v|²⟩ and v_r(r) are quantities independent of μ_r, and can be measured in radial bins.

The bottom panel of Fig. 4 shows the normalised template functions for σ_v||(r) measured for each of the four redshift bins. The measured amplitudes at large distances, σ_v||, are listed in Table 1, where we see that the dispersion is higher at lower redshift, consistent with the growth of large-scale structure. From the figure we also see that the dispersion decreases in the void interior, where the potential is steep and the bulk outflow motion dominates.

The final template function required is that for ξ^rr(r). In the two realistic reconstruction scenarios we test, we assume that the true real-space information is not available, so that ξ^rr cannot be directly measured. Instead we used the quantity ξ^pp(r), which could be measured using the approximate real-space catalogues, as an estimator for ξ^rr(r). This is the same approach as used by Nadathur et al. (2019b, 2020b) and Woodfinden et al. (2022) except that here we do not have a large number of mock catalogues over which we can average this estimator. This greatly increases the noise in the estimate, which is especially relevant when attempting to determine the dependence of ξ^pp on the reconstruction parameter β. For the current work we therefore decided to fix this estimator to ξ^pp(β^fid) evaluated at the fiducial β value given in Table 1. This approach will be revisited for a future full Euclid analysis, by which time we expect a full complement of mocks to be available.

As mentioned in Sect. 2.3, our general code implementation seeks to preserve information on possible anisotropies in ξ^pp(r) which may be introduced by the AP effect or imperfections in reconstruction, without imposing the assumption of spherical symmetry, ξ^pp(r) = ξ^pp(r). In practice we achieve this by measuring the monopole ${\xi }_0^{\text{pp}}(r)$, quadrupole ${\xi }_2^{\text{pp}}(r)$ and hexadecapole ${\xi }_4^{\text{pp}}(r)$ moments and using them to reconstruct ξ^pp(r). In the idealised perfect reconstruction scenario, things are much simpler and we just use the true real-space information from Flagship and directly measure the multipoles ${\xi }_0^{\text{rr}}(r)$, ${\xi }_2^{\text{rr}}(r)$ and ${\xi }_4^{\text{rr}}(r)$ from the data.

4.3. Model fitting and parameter inference

We compare the model to the measured redshift-space CCF in terms of the compression to multipoles. In each redshift bin, we define the data vector ${\mathit{\xi }}^{\text{s}}\equiv \left({\xi }_0^{\text{rs}},{\xi }_2^{\text{rs}},{\xi }_4^{\text{rs}}\right)$ formed by concatenating the monopole, quadrupole and hexadecapole moments of ξ^rs(s, μ_s), measured in 30 radial bins. In the two realistic reconstruction scenarios where the true real-space information is not used, we use the CCF ξ^ps(s, μ_s) measured using the void positions identified in the post-reconstruction, approximate real-space galaxy distribution to construct ξ^s. In this case the data vector inherits a dependence on the reconstruction parameter β from the void positions.

To determine the uncertainty in this measurement of ξ^s, we used a jackknife resampling method to estimate the covariance matrix. For each of the four redshift bins we divided the Flagship data sample into N_JK = 900 non-overlapping sub-regions of equal volume. We iterated over all sub-boxes, in each case estimating the cross-correlation function ξ^s(k) by removing all the void-galaxy pairs for which both void and galaxy lie in the kth sub-box, and reweighting the ones where either void or galaxy lie within the sub-box using the matched weighting scheme prescribed by Mohammad & Percival (2022) and implemented in pycorr. We then used these realisations to estimate the covariance between the ith and jth bins of the data vector

$$\begin{array}{c}\hfill {\mathsf{C}}_{\mathit{ij}}^{\mathrm{data}}=\frac{N_{\mathrm{JK}}-1}{N_{\mathrm{JK}}}{\displaystyle \sum _{k=1}^{N_{\mathrm{JK}}}}({\mathit{\xi }}_i^{\mathrm{s}(k)}-\overline{{\mathit{\xi }}_i^{\mathrm{s}}})({\mathit{\xi }}_j^{\mathrm{s}(k)}-\overline{{\mathit{\xi }}_j^{\mathrm{s}}}),\end{array}$$(25)

where ${\overline{{\mathit{\xi }}^{\mathrm{s}}}}_i=\frac{1}{N_{\mathrm{JK}}}\sum _{k=1}^{N_{\mathrm{JK}}}{\mathit{\xi }}_i^{\mathrm{s}(k)}$ is the mean of the N_JK jackknife realisations. The correlation structure of the resultant covariance matrix is shown in Fig. 5. Due to the limitations of having only the single Flagship mock, there will be an additional correlation between the model prediction and the measured data vector that has been neglected in this estimate of the covariance but is discussed in Appendix A.

Fig. 5. Normalised covariance matrix Cdata for the data vector from the Flagship zeff = 1.0 redshift bin, estimated using jackknife resampling. Dashed lines differentiate between monopole, quadrupole and hexadecapole blocks of the data vector.

At any point in parameter space, we computed the model prediction for ξ^s as described in Sect. 2 using the public code victor⁹. We assume that the likelihood has a Gaussian form,

$$\begin{array}{c}\hfill log\mathcal{L}=-\frac{1}{2}({\mathit{\xi }}_{\mathrm{theory}}^{\mathrm{s}}-{\mathit{\xi }}_{\mathrm{data}}^{\mathrm{s}}){\mathsf{C}}^{-1}{({\mathit{\xi }}_{\mathrm{theory}}^{\mathrm{s}}-{\mathit{\xi }}_{\mathrm{data}}^{\mathrm{s}})}^{\mathrm{T}}.\end{array}$$(26)

In the realistic reconstruction scenarios we consider, this likelihood evaluation depends on four parameters: explicitly on fσ₈, σ_v|| and ϵ, and implicitly on the reconstruction parameter β via its effect on the recovered void positions in the post-reconstruction field and thus on the measured correlations. In the idealised case of perfect reconstruction there is no β dependence. We sampled the full parameter space using the interface between victor and the Monte Carlo Markov chain (MCMC) sampling code Cobaya¹⁰ (Torrado & Lewis 2021). Repeating reconstruction, void-finding and CCF measurement at each β value along the chain would be too slow for the MCMC so instead we followed the method of previous works (Nadathur et al. 2019b, 2020b) and pre-computed the correlations on a grid of β values around the fiducial β^fid and interpolated the data vector at the time of evaluation of the likelihood¹¹. To avoid extrapolating outside of this grid, the prior on β was limited to the same range as the interpolation; we used uninformative flat priors for all of the other parameters. Additionally, convergence of the MCMC chains is checked using the Gelman-Rubin R − 1 statistic, using the convergence criterion R − 1 < 0.01, and 20% of the initial steps were discarded as burn-in.

5. Results

5.1. Fits to Flagship data

In Fig. 6 we show the measured monopole, quadrupole and hexadecapole moments of the redshift-space CCF ξ^ps and the predictions for the best-fit model, together with the residuals, in each of the redshift bins for the realistic reconstruction case and for the catalogue with redshift errors. It can be seen that the model describes the observed multipoles well, although the hexadecapole moment does not contain much information except in the lowest redshift bin where the density of galaxies and voids is highest. This is the most realistic of the scenarios considered, but a similar quality of fit was obtained for the idealised perfect reconstruction scenario and when neglecting redshift errors.

Fig. 6. Measured multipole moments of the redshift-space void-galaxy CCF, ${\xi }_{\mathcal{l}}^{\text{ps}}$, for realistic reconstruction applied to the catalogue including redshift errors, are shown as the data points. Error bars shown are determined from the diagonal entries of the covariance matrix Cdata. The solid lines show the corresponding best-fit models, and insets show the residuals of the model with respect to the data. The rows correspond to (from top to bottom) the monopole, quadrupole and hexadecapole moments (ℓ = 0, 2, 4). Columns correspond to the different redshift bins, labelled with the effective redshift zeff.

The marginalised one-dimensional posterior constraints on the growth rate fσ₈ and the AP distortion parameter ϵ individually obtained from analysis of the Flagship mock are summarised in Table 2. For the idealised case of perfect reconstruction, we recover unbiased constraints on both fσ₈ and ϵ in all redshift bins. We find a measurement precision on ϵ of about 0.5–0.6% in each individual redshift bin, and the precision on fσ₈ varies between ∼8.5% for the lowest redshift bin and ∼12% in the two highest redshift bins. Using realistic reconstruction and adding redshift errors to the mocks does not significantly degrade these constraints, and the central values recovered are consistent to within the stated uncertainty in each case. The reconstruction method used is therefore not causing a significant loss of information.

The marginalised 1D and 2D posteriors on fσ₈ and ϵ are shown in Fig. 7, where the dashed crosshairs indicate the fiducial values for the Flagship cosmology, which are consistently recovered. Comparing the perfect and realistic reconstruction contours, the largest shift (of ∼1σ) is seen in the lowest redshift bin. The introduction of redshift errors does not cause significant additional shifts in any redshift bin.

Fig. 7. Marginalised 1D and 2D posterior constraints for parameters fσ8 and ϵ from fits to the Flagship void-galaxy CCF measurements in four redshift bins. Contours show the 68.3% (1σ) and 95.5% (2σ) confidence intervals obtained in the case of perfect reconstruction (blue), realistic reconstruction (light red), and realistic reconstruction with added redshift error (dark red). Dashed crosshairs indicate the true values for the Flagship cosmology.

Figure 8 shows the marginalised 1D measurements of the growth rate and ratio of the transverse comoving distance D_M to the Hubble distance D_H, derived from the ϵ constraints, from the Flagship data as a function of redshift. These measurements are compared to the 68.3% (1σ) and 95.5% (2σ) confidence intervals on these quantities derived from extrapolating the CMB constraints from Planck Collaboration VI (2020) down to these redshifts while assuming a ΛCDM cosmological model. Results for ϵ itself are shown in the inset plot.

Fig. 8. Redshift dependence of the measured cosmological parameters obtained from fits to Flagship simulation data. The top panel shows constraints on the ratio of comoving transverse and Hubble distances divided by redshift z, DM/(zDH), which is derived from ϵ (inset). The bottom panel shows measurements of the growth rate fσ8. Coloured data points indicate the different reconstruction scenarios considered, plotted at slightly shifted redshift locations for clarity. The grey shaded bands represent the 68.3% (1σ) and 95.5% (2σ) confidence intervals from fits to Planck Collaboration VI (2020), extrapolated down to these redshifts assuming ΛCDM, and centred on the Flagship cosmology.

5.2. Full Euclid volume

The planned final footprint of Euclid is 15 000 deg², almost three times larger than that of the Flagship mock used in the analysis above. The constraints we expect from using the full Euclid volume will therefore be correspondingly better. A simple scaling of the covariance matrix suggests that it should shrink proportional to the increase in the volume, that is by a factor of 2.91. This corresponds to measurement uncertainties on the binned multipole moments that should be over 40% smaller than those shown in Fig. 6. We first confirmed this scaling with survey volume by computing covariance matrices for eight volume subsections of the Flagship simulation, ranging from 0.125 to 0.875. We confirmed that the terms of C scale roughly linearly with the survey volume. We therefore obtained the expected covariance matrix ${\mathsf{C}}_{\mathit{ij}}^{\mathit{Euclid}}$ for the full Euclid survey by scaling the one obtained in Sect. 4.3 by a factor of 1/2.91.

To test how this reduction in uncertainties in the measurement of the void-galaxy CCF propagates through to the recovered cosmological parameter constraints, we repeated the parameter inference step with this rescaled covariance and using a synthetic data vector generated assuming this level of measurement noise. To do this, we considered the idealised perfect reconstruction scenario. Using the real-space CCF multipoles ${\xi }_{\mathcal{l}}^{\text{rr}}$(r) measured from Flagship, we computed the fiducial model prediction, ξ^s, th, for the redshift-space data vector. We then created a synthetic data vector

$$\begin{array}{c}\hfill {\mathit{\xi }}^{\mathrm{s},\mathrm{m}\mathrm{ock}\left\}\right\}}=\mathsf{L}\mathit{Z}+{\mathit{\xi }}^{\mathrm{s},\mathrm{t}\mathrm{h}\left\}\right\}},\end{array}$$(27)

where LL^T = C^Euclid, L is obtained by a Cholesky decomposition of the scaled covariance matrix C^Euclid, and Z is a vector of independent standard normal random variables. We then repeated the parameter inference of Sect. 4.3 fitting to ξ^s, mock using the rescaled covariance matrix.

The results of this fitting procedure provide the forecast for Euclid constraints from the void-galaxy CCF measurement, and form the main result of this paper. The forecasts are summarised in Table 3 and shown as a function of redshift in Fig. 9, compared to existing constraints from the application of the same void-galaxy method on current data. As expected, we recover the fiducial cosmological parameters up to the 68.3% (1σ) statistical uncertainty in each redshift bin. Comparing Tables 2 and 3 shows that while the forecast uncertainty in ϵ scales approximately as the $1/\sqrt{2.91}$ factor expected from the simple volume scaling, the forecast constraints on fσ₈ are slightly more pessimistic. Our final results suggest that with Euclid voids we will be able to measure the distance ratio D_M/D_H to a precision of 0.3–0.4% in four redshift bins over the range 0.9 ≤ z < 1.8. This is comparable even to the extrapolated constraints from Planck which assume ΛCDM. Voids will also provide a precision of about 5–8% in measurement of the growth rate over these redshifts, a large increase in the statistical power of this method compared to current surveys, and a useful consistency check for results from galaxy clustering.

Fig. 9. Redshift dependence of the forecast measurements that could be obtained from the 15 000 deg2Euclid survey (purple data points). These are generated from a random realisation of the Euclid CCF (Sect. 5.2) and so the central values can be shifted relative to the Flagship results in Fig. 8. Low redshift measurements from current SDSS data obtained using the same analysis method by Nadathur et al. (2020b) and Woodfinden et al. (2022) are included for comparison.

Figure 9 highlights the key point that AP constraints from voids will be competitive with those derived from model-dependent extrapolations to low redshift from the CMB (and will exceed those that can be obtained from galaxy clustering alone). AP measurements therefore present the best potential for information gain from using the void-galaxy CCF. On the other hand, while the constraints on fσ₈ are weaker than those that can be obtained from galaxy clustering alone and do not add much cosmological information, they can serve as useful consistency checks.

5.3. Cosmology forecasts

We now investigate the impact on cosmological model forecasts of the expected measurement precision that can be achieved with Euclid voids, summarised in Table 3. Euclid Collaboration (2020) have presented forecasts for the constraints on model parameters that can be achieved with the primary Euclid probes, namely spectroscopic galaxy clustering, weak lensing, photometric galaxy clustering, and their cross-correlations. In this work, we compare voids to these benchmarks as a separate standalone probe, and do not investigate the combination of voids with other observables (see, for instance, Nadathur et al. 2019b, 2020b; Bonici et al. 2023 for examples of such combinations).

In this scenario, voids probe cosmology through the AP measurement of D_M(z)/D_H(z) and the stringent constraints this places on the background expansion history. In contrast, the growth rate constraints are weaker than those that can be obtained from galaxy clustering. These measurements can be used as tests of consistency of the model in different regimes, but in the absence of any discrepancy they do not add much cosmological information, even when voids are combined with other probes (Nadathur et al. 2020b). When treating voids as a standalone probe as we do here, they do not provide sufficient constraints to be of interest. We therefore neglect the fσ₈ results and focus entirely on the geometrical AP constraints. To do so, we centred the D_M(z)/D_H(z) values at their expected values for the Flagship cosmology and used the measurement uncertainties from Table 3.

Within the flat ΛCDM model scenario, measurement of D_M/D_H at a given redshift directly translates to a constraint on the matter density parameter Ω_m (or, equivalently, on Ω_Λ = 1 − Ω_m). The precision that we can achieve on D_M(z)/D_H(z) with voids corresponds to Ω_m = 0.3183 ± 0.0028 (at the 68.3% confidence level), a relative statistical uncertainty of ±0.009. The corresponding forecast relative uncertainty from Euclid spectroscopic galaxy clustering in the optimistic (pessimistic) forecast scenario is 0.013 (0.021), from weak lensing is 0.012 (0.018), and from the combination of the two is 0.006 (0.009) (Euclid Collaboration 2020). Thus for this model and parameter, Euclid voids alone seem to outperform both galaxy clustering and weak lensing individually, and look to be competitive with their combination. It is worth noting, however, that this forecast is based only on the statistical constraints, obtained by rescaling a Flagship covariance up to the full Euclid volume.

We also consider a minimal one-parameter extension to the ΛCDM model known as wCDM, in which the cosmological constant Λ is replaced by a dark energy component with a constant equation of state that is not restricted to w = −1. Fig. 10 shows the forecast constraints obtained on Ω_m and w in the dark red contours. We recover $w=-1.{00}_{-0.05}^{+0.06}$ and Ω_m = 0.3183  ±  0.0028, corresponding to a relative precision of 6% on w. For comparison, the yellow contours in Fig. 10 show the current constraints on these parameters from SDSS at multiple redshifts between z = 0.15 and z = 2.33, comprising the BAO results from different galaxy and quasar samples reported by Alam et al. (2021) together with recent improved high-redshift AP measurements from the ‘full shape’ of the Lyman-α forest clustering (Cuceu et al. 2023), which give w = −0.90  ±  0.12. Due to a well-known geometrical degeneracy which is broken by low-redshift AP measurements, the constraints on this model from the CMB (including CMB lensing) are relatively weak, $w=-1.{57}_{-0.40}^{+0.50}$ (Planck Collaboration VI 2020). For completeness, Fig. 10 also shows two sets of forecasts for the Euclid void-galaxy CCF obtained by Hamaus et al. (2022) through a different method for modelling and measurement to that used in this work.

Fig. 10. Forecast constraints on the matter density parameter Ωm and the dark energy equation of state w in a wCDM cosmological model from a measurement of the void-galaxy CCF with Euclid as outlined in this work (red contours), corresponding to $w=-1.{00}_{-0.05}^{+0.06}$. We also show previous forecasts for the Euclid void-galaxy CCF – obtained using a different model and measurement method – by Hamaus et al. (2022), for the ‘independent’ (green contours) and ‘calibrated’ (blue, labelled ‘cal.’ here) cases in that paper, corresponding to whether nuisance parameters are marginalised over or fixed. For context, we show current measurements from current BAO and ‘full-shape’ measurements in the SDSS galaxy, quasar and Lyman-α datasets from the MGS, BOSS and eBOSS surveys (yellow) presented by Cuceu et al. (2023).

These results highlight the excellent promise of the void-galaxy CCF as a standalone cosmological probe. When removing the assumption of flatness and allowing more general dark energy models with additional degrees of freedom, the void measurement of D_M/D_H is no longer sufficient to provide useful model constraints on its own. However, the higher precision of the AP measurement still provides large gains in information and improves the figure of merit for dark energy in these models when voids are used in combination with other probes, as shown by Nadathur et al. (2020a). Thus this technique retains great potential to enhance Euclid science.

6. Conclusions

We have presented a forecast of the cosmological constraints that can be obtained from measurement of the anisotropic void-galaxy cross-correlation function (CCF) in the spectroscopic sample of the upcoming Euclid galaxy survey. The pattern of anisotropies in the redshift-space CCF can be used to simultaneously fit for AP distortions caused by differences between the fiducial analysis model and the true background cosmology, and for RSD caused by the peculiar motions of galaxies around voids. Similar cosmological analyses of the void-galaxy CCF have previously been applied to data from the Sloan Digital Sky Survey (Nadathur et al. 2019b, 2020b; Hamaus et al. 2020; Woodfinden et al. 2022), and have been shown to provide significantly better AP measurements than obtained from traditional galaxy-clustering analyses, including BAO. Precise AP measurements, particularly at high redshift, provide very powerful probes of cosmological models of dark energy (Nadathur et al. 2020a; Cuceu et al. 2023).

Our forecast is based on an analysis of the lightcone Flagship galaxy mock catalogue, which was calibrated to match the expected properties of the Euclid spectroscopic galaxy sample, and covers 5157 deg² of the sky between the redshifts 0.9 ≤ z < 1.8. We analysed this mock data using the same method previously used by Nadathur et al. (2019b) and Woodfinden et al. (2022). A key component of this method is the use of velocity field reconstruction to approximately remove large-scale RSD from the galaxy distribution before identifying voids. This step is necessary to remove a selection bias in the construction of the void sample that depends on the orientation of the voids to the line-of-sight direction, thus ensuring the validity of the theoretical modelling. To assess the efficiency of the particular algorithm used to perform this reconstruction, we compared the results obtained using the true real-space galaxy field – corresponding to the idealised perfect reconstruction scenario, that will not be achievable in practice with Euclid data – to those obtained with realistic reconstruction, with and without including spectroscopic redshift errors in the mock data.

We performed the analysis in four non-overlapping redshift bins, compressing the anisotropic CCF information to its first three even-order multipoles – monopole, quadrupole, and hexadecapole – which contain all of the practically measurable information. We found that using a realistic reconstruction algorithm did not significantly increase the statistical uncertainty of the recovered marginalised constraints on the cosmological parameters ϵ and fσ₈ compared to the idealised scenario, and also did not introduce significant systematic offsets at the level of statistical precision available from the Flagship mock. Realistic spectroscopic redshift errors were found to have a similarly minor effect on the results. However, we caution that the analysis here was limited by the availability of only a single simulation realisation, so we can not guarantee that systematic offsets will not be present when dealing with Euclid data. A more detailed study of systematic errors should be performed when a larger number of Euclid mocks become available.

We then extrapolated from the results obtained from the Flagship mock to the full Euclid survey, covering almost three times larger survey volume, to obtain our forecast for the achievable statistical precision. We found that void-galaxy CCF measurements with Euclid should be able to achieve a statistical precision of down to 0.3% on the measurement of the distance ratio D_M(z)/D_H(z) over the redshift range 0.9 ≤ z < 1.8, and 5–8% on the growth rate parameter fσ₈. Within the context of flat ΛCDM cosmological models, the statistical uncertainty we achieve on D_M(z)/D_H(z) translates to a relative statistical uncertainty of 0.009 on the matter density parameter Ω_m which is ∼30% better than forecasts of what is achievable from both the Euclid spectroscopic galaxy clustering and weak lensing probes individually, and is comparable to their combination (Euclid Collaboration 2020). In the wCDM model scenario, we found that voids alone can measure the dark energy equation of state to a relative precision of 6%, $w=-1.{00}_{-0.05}^{+0.06}$.

It should be noted here that the modelling method adopted for this forecast is the template-fitting approach used in several previous works. In this approach template functions describing the real-space CCF, the void velocity profile (or equivalently, matter density profile) and the velocity dispersion profile are first constructed by calibrating against simulation results and then allowed to vary with cosmological parameters, as described in Sect. 2. The recent work of Massara et al. (2022) suggests that improvements to this template-fitting approach may be required in the future, but our aim is to forecast the constraints possible with currently available techniques so we have not considered these here. We have also not accounted for possible systematic effects which may enter this analysis if the simulations from which the templates are constructed differ strongly from the true cosmological model. A quantitative study of these effects (such as conducted in the context of SDSS by Woodfinden et al. 2022) would require mock catalogues constructed at different cosmologies, but at present we only have access to a single realisation of the Euclid Flagship mock at a fixed cosmology. A large number of mock catalogues would also allow a better estimate of the covariance matrix than the jackknife resampling technique used here.

In this work, the template we used for the real-space CCF was necessarily estimated from the single realisation of the Flagship data itself. This increases the noise in this estimate, in a manner that is correlated with measurement noise in the data vector. When the final Euclid data are analysed, we anticipate that large numbers of mock realisations will be available, so that while a similar procedure could be used in the final analysis, it will not be mandatory. If indeed both real-space and redshift-space CCFs are measured from the same data, it would be necessary to account for the correlation between them for a fully self-consistent analysis. We describe how to do this in Appendix A, although we choose not to include this in the main analysis presented above. However, as shown in Appendix A, the effect of accounting for this correlation is to decrease the effective statistical uncertainties on the recovered cosmological parameters. The forecasts using the larger statistical uncertainties presented in our headline analysis are therefore expected to be conservative.

Our treatment of possible observational systematics is quite simple and ignores some effects that we expect will be present in the Euclid data. In particular, although we have incorporated the effects of a uniformly low completeness of the spectroscopic sample and a Gaussian distribution of redshift errors, we have not considered the purity of the sample and errors due to spectral line misidentification. These effects will be explored in detail in a future work. This paper is one of several companion papers forecasting the constraining power from different void observables in Euclid (Hamaus et al. 2022; Contarini et al. 2022; Bonici et al. 2023) which all use the same common baseline of assumptions, to allow a consistent comparison across methods.

Hamaus et al. (2022) have also investigated forecasts for constraints obtained from the void-galaxy CCF with Flagship, using a different approach. Unlike in our work, that paper did not use reconstruction to try to remove the orientation-dependent sample selection bias in the void catalogues, and thus obtained sharply different measurements of the CCF than those we find here. They then introduced additional nuisance parameters and modified the growth-dependent terms of the theoretical model derived from the conservation equation in Sect. 2 in order to describe the Flagship data. Other minor differences in approach include the assumption that the matter density profile around voids is simply related to the galaxy profile by the large-scale linear galaxy bias b (instead of our template method for predicting the mean galaxy velocity), and neglecting velocity dispersion around the mean.

Despite these differences, the forecasts of Hamaus et al. (2022) are broadly compatible with ours. In their baseline model (which they refer to as ‘independent’), which marginalises over nuisance parameters as we do here, they predict a 0.5% measurement of the AP parameter, somewhat weaker than our 0.3%, and a 4% measurement of the growth rate, slightly better than our value of 5%–8%. They further report that if the nuisance parameters are instead fixed to their best-fit values (referred to as the ‘calibrated’ case), these results could improve to 0.4% and 1% precision on the AP parameter and growth rate, respectively. The additional improvement in the AP measurement that is possible from our method can have a large effect on the cosmological constraints obtained, as shown in Fig. 10, where uncertainties in w and Ω_m are significantly reduced. Taken together, the results of these two papers showcase both the promise and the robustness of the void-galaxy cross-correlation method applied to Euclid data.

Acknowledgments

SR acknowledges support from the Research Council of Norway through project 325113. SN acknowledges support from an STFC Ernest Rutherford Fellowship, grant reference ST/T005009/2. AW acknowledges the support of the Natural Sciences and Engineering Research Council of Canada (NSERC) [funding reference number 547744]. Research at Perimeter Institute is supported in part by the Government of Canada through the Department of Innovation, Science and Economic Development Canada and by the Province of Ontario through the Ministry of Colleges and Universities. This work has made use of CosmoHub (Carretero et al. 2017; Tallada et al. 2020). 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 and IEEC), and was partially funded by the “Plan Estatal de Investigación Científica y Técnica y de Innovación” programme 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: Accounting for correlation between data and theory

Equation (25) defines the jackknife estimate of the covariance C^data of the redshift-space data vector ${\mathit{\xi }}_{\mathrm{data}}^{\mathrm{s}}$, which we have used for computing the likelihood in the main analysis. However, this is not the only source of uncertainty: the computation of the theory model ${\mathit{\xi }}_{\mathrm{theory}}^{\mathrm{s}}$ depends on the multipoles of the real-space CCF ${\xi }_{\mathcal{l}}^{\text{rr}}$, which are themselves determined from measurements performed on the galaxy catalogue using the estimators described in Sect. 4.2. These estimators have an associated uncertainty, which is negligible when ${\xi }_{\mathcal{l}}^{\text{rr}}$ is estimated from the mean of a large number (N ≈ 1000) of mock realisations as done, for example, by Nadathur et al. (2019b) and Woodfinden et al. (2022), but is significant in our case as the estimate is derived from a single Flagship realisation. Equally, since ${\xi }_{\mathcal{l}}^{\text{rr}}$ is estimated from the same Flagship void and galaxy catalogues which are used to measure the data vector, the uncertainties in ${\mathit{\xi }}_{\mathrm{theory}}^{\mathrm{s}}$ and ${\mathit{\xi }}_{\mathrm{data}}^{\mathrm{s}}$ will in principle be correlated with each other.

In order to self-consistently account for these model uncertainties and correlations with the data vector, we can define a combined effective data vector ${\mathit{\xi }}^{\mathrm{tot}}\equiv {\mathit{\xi }}_{\mathrm{theory}}^{\mathrm{s}}-{\mathit{\xi }}_{\mathrm{data}}^{\mathrm{s}}$ and estimate the covariance C^tot using the same jackknife procedure described in Sect. 4.3. Then Eq. (26) for the log-likelihood remains unchanged, provided that C⁻¹ is understood to refer to the inverse covariance matrix of ξ^tot and not ${\mathit{\xi }}_{\mathrm{data}}^{\mathrm{s}}$ alone. We can write this schematically as:

$$\begin{array}{cc}& \mathsf{C}(\mathrm{data}-\mathrm{theory})=\mathsf{C}\left(\mathrm{data}\right)+\mathsf{C}\left(\mathrm{theory}\right)\hfill \\ \hfill & -2\mathsf{C}(\mathrm{data},\mathrm{theory}),\hfill \end{array}$$(A.1)

where the last term represents the cross-covariance between the data vector and the theory model. The familiar scenario corresponds to a situation where the second term on the right is negligible and the last term completely vanishes. This is the case by construction in previous studies in the literature (e.g. Nadathur et al. 2020b; Woodfinden et al. 2022), since ${\xi }_{\mathcal{l}}^{\text{rr}}$ was taken as the mean from many mocks and thus data were theory are not correlated. For the method used in this work, given the limitations of only a single Flagship mock sample, this is not true. For future Euclid data analyses, whether this term is important or not will depend on the choice of how to determine ${\xi }_{\mathcal{l}}^{\text{rr}}$.

For Flagship, we found a high degree of correlation between the model and the data vectors, which changes the correlation structure of the full covariance C^tot with respect to C^data. This is shown in Fig. A.1, and can be compared to Fig. 5. Two features are immediately noticeable: off-diagonal correlations between the monopole and quadrupole components in particular are enhanced, and within each diagonal block of the covariance matrix, correlations between different radial bins are relatively suppressed. This behaviour arises because uncertainties in the measurement of the real and redshift-space monopoles ${\xi }_0^{\text{rr}}$ and ${\xi }_0^{\text{rs}}$ in each bin are naturally strongly correlated, as the real and redshift-space galaxy positions are closely correlated. In each radial bin, the monopole component ${\xi }_0^{\text{rr}}$ strongly influences the corresponding modelled multipoles of ${\mathit{\xi }}_{\mathrm{theory}}^{\mathrm{s}}$, thus creating the correlation between the elements of ${\mathit{\xi }}_{\mathrm{theory}}^{\mathrm{s}}$ and the monopole component of ${\mathit{\xi }}_{\mathrm{data}}^{\mathrm{s}}$.

Fig. A.1. Normalised full covariance matrix Ctot from the Flagship zeff = 1.0 redshift bin, estimated using jackknife resampling. Dashed lines differentiate between monopole, quadrupole and hexadecapole blocks of the data vector.

Even more importantly, the strong correlation between real- and redshift-space multipole moments causes the amplitude of the diagonal elements of C^tot to be significantly smaller than those of C^data. This is shown in Fig. A.3. At first glance, it may seem unusual that accounting for additional uncertainty in the model effectively reduces the total uncertainty C^tot, but this is simply a reflection of the correlation that means that variations in ${\mathit{\xi }}_{\mathrm{theory}}^{\mathrm{s}}-{\mathit{\xi }}_{\mathrm{data}}^{\mathrm{s}}$ are reduced since uncertainties move both components in the same direction. It is intuitive that when ${\mathit{\xi }}_{\mathrm{theory}}^{\mathrm{s}}$ and ${\mathit{\xi }}_{\mathrm{data}}^{\mathrm{s}}$ are highly correlated with each other, models that deviate from observation should be more severely penalised in the log-likelihood, and this is achieved by the smaller values of C^tot. An equivalent way to view this is that when ${\mathit{\xi }}_{\mathrm{theory}}^{\mathrm{s}}$ and ${\mathit{\xi }}_{\mathrm{data}}^{\mathrm{s}}$ are both obtained from the same realisation of the initial conditions, the cosmic variance in them cancels out.

Changing the covariance matrix used in the likelihood evaluation in this way naturally propagates through to the recovered statistical uncertainty on the cosmological parameters. To estimate this, we performed all the model fits again using the appropriate full covariance C^tot instead of C^data. The results obtained for D_M/D_H and fσ₈ are shown in Fig. A.2, where the black points are for the perfect reconstruction scenario fit using C^tot, and the blue points are for fits to the same data using C^data. It is clear from this that correctly accounting for the correlation between model and data leads to a very significant increase in the statistical precision of the fit, roughly by a factor of 3 in both D_M/D_H and fσ₈. The main results of our paper, which do not account for this effect, are therefore conservative overestimates of the statistical error that can be achieved with Euclid if the real-space cross-correlation function is measured from the data instead of being taken from a simulation.

Fig. A.2. Redshift dependence of the measured cosmological parameters obtained from Flagship using different covariance prescriptions. The blue points show results for fits to the data vector for the perfect reconstruction case using Cdata and are the same as in Fig. 8. Black points are for fits to the same data vector but using the full covariance Ctot. Light red points are for the fit to the realistic reconstruction case when using Ctot. The use of the full covariance Ctot greatly reduces the statistical uncertainty but can introduce systematic offsets when combined with the use of reconstruction.

However, the reduction in statistical uncertainty makes systematic errors more important. While in the perfect reconstruction scenario we are able to recover unbiased estimates of D_M/D_H and fσ₈ at the higher precision, this is not so for the realistic reconstruction scenario (red points in Fig. A.2), where systematic offsets from the fiducial values are apparent. The source of these offsets is likely to be the residuals in the modelling of the monopole moment at small separations, an example of which can be seen in Fig. 6. These are caused by the imperfections of the practical reconstruction technique. Figure A.3 shows that the relative reduction in covariance terms from using C^tot instead of C^data is largest for the monopole, so these residuals have a larger effect on the recovered fit. For the realistic reconstruction scenario, the full covariance matrix C^tot therefore does not give a good estimate of the total error budget, and systematic uncertainties would have to be folded in.

Fig. A.3. Diagonal elements of the full covariance matrix Ctot (black lines), as well as the three different contributions from Eq. (A.1), showing the reduced variance that results from properly accounting for correlation between data and theory. The three different panels represent the monopole, quadrupole and hexadecapole part, respectively.

As we only have one Flagship simulation available, we cannot quantitatively assess the true size of these systematic uncertainties. However, it appears that the apparent large improvement in statistical precision would not be practically realised for Euclid in the reconstruction scenarios, even if ${\xi }_{\mathcal{l}}^{\text{rr}}$ were estimated from the data and the full covariance C^tot were used. Therefore the pessimistic forecasts we have obtained for the existing methodology using C^data alone are more realistic.

All Tables

All Figures

Fig. 1. Left: distribution of void sizes in each of the four redshift bins of the Flagship mock catalogue. The vertical bars indicate the void size cuts applied to each catalogue, and dashed lines represent the voids which are discarded. Right: comoving number density n(z) of galaxies and voids as function of redshift in Flagship. For voids, n(z) is multiplied by 100 for visibility. The edges of the redshift bins are shown by dashed vertical lines. The decrease in the number of voids near the edges of each bin arises because void-finding is performed on each bin individually as discussed in the text; some void candidates near the redshift edges of the bin are then missed.

In the text

Fig. 2. Quadrupole moments of the galaxy-galaxy power spectrum in each redshift bin of the Flagship mock, measured in redshift space (yellow), real space (blue) and after applying reconstruction to remove RSD from the galaxy positions (red). Reconstruction was performed using the fiducial values of βfid listed in Table 1 and using smoothing scales Rs of 8 and 9 h−1Mpc for the first two and last two bins, respectively. The similarity of the red and blue lines and their proximity to zero is because reconstruction is successfully removing the large-scale RSD-induced anisotropy in the galaxy positions.

In the text

Fig. 3. Void-galaxy CCF in real space, measured in the Flagship simulation zeff = 1.0 redshift bin. From top to bottom, the panels show the monopole, quadrupole and hexadecapole moments of the CCF. The multipoles of ξrr, corresponding to the perfect reconstruction scenario and measured using the true real-space positions are shown in blue, while light and dark red respectively show the multipoles of ξpp obtained using reconstruction without and with redshift errors. The inset panels show the residuals with respect to the perfect reconstruction case of ${\xi }_{\mathcal{l}}^{\text{rr}}$ for each multipole. In all cases ξ2 and ξ4 are close to zero, indicating that reconstruction successfully removes orientation-dependent selection bias.

In the text

	Fig. 4. Measured radial outflow velocity profile vr(r) (top), and normalised line-of-sight velocity dispersion σv||(r) (bottom) of galaxies around voids found in the four redshift bins of Flagship. Shading represents the 68.3% (1σ) measurement uncertainty coming from Poisson noise. The dispersion amplitude σv for the four bins can be found in Table 1.
In the text

	Fig. 5. Normalised covariance matrix Cdata for the data vector from the Flagship zeff = 1.0 redshift bin, estimated using jackknife resampling. Dashed lines differentiate between monopole, quadrupole and hexadecapole blocks of the data vector.
In the text

Fig. 6. Measured multipole moments of the redshift-space void-galaxy CCF, ${\xi }_{\mathcal{l}}^{\text{ps}}$, for realistic reconstruction applied to the catalogue including redshift errors, are shown as the data points. Error bars shown are determined from the diagonal entries of the covariance matrix Cdata. The solid lines show the corresponding best-fit models, and insets show the residuals of the model with respect to the data. The rows correspond to (from top to bottom) the monopole, quadrupole and hexadecapole moments (ℓ = 0, 2, 4). Columns correspond to the different redshift bins, labelled with the effective redshift zeff.

In the text

Fig. 7. Marginalised 1D and 2D posterior constraints for parameters fσ8 and ϵ from fits to the Flagship void-galaxy CCF measurements in four redshift bins. Contours show the 68.3% (1σ) and 95.5% (2σ) confidence intervals obtained in the case of perfect reconstruction (blue), realistic reconstruction (light red), and realistic reconstruction with added redshift error (dark red). Dashed crosshairs indicate the true values for the Flagship cosmology.

In the text

Fig. 8. Redshift dependence of the measured cosmological parameters obtained from fits to Flagship simulation data. The top panel shows constraints on the ratio of comoving transverse and Hubble distances divided by redshift z, DM/(zDH), which is derived from ϵ (inset). The bottom panel shows measurements of the growth rate fσ8. Coloured data points indicate the different reconstruction scenarios considered, plotted at slightly shifted redshift locations for clarity. The grey shaded bands represent the 68.3% (1σ) and 95.5% (2σ) confidence intervals from fits to Planck Collaboration VI (2020), extrapolated down to these redshifts assuming ΛCDM, and centred on the Flagship cosmology.

In the text

Fig. 9. Redshift dependence of the forecast measurements that could be obtained from the 15 000 deg2Euclid survey (purple data points). These are generated from a random realisation of the Euclid CCF (Sect. 5.2) and so the central values can be shifted relative to the Flagship results in Fig. 8. Low redshift measurements from current SDSS data obtained using the same analysis method by Nadathur et al. (2020b) and Woodfinden et al. (2022) are included for comparison.

In the text

Fig. 10. Forecast constraints on the matter density parameter Ωm and the dark energy equation of state w in a wCDM cosmological model from a measurement of the void-galaxy CCF with Euclid as outlined in this work (red contours), corresponding to $w=-1.{00}_{-0.05}^{+0.06}$. We also show previous forecasts for the Euclid void-galaxy CCF – obtained using a different model and measurement method – by Hamaus et al. (2022), for the ‘independent’ (green contours) and ‘calibrated’ (blue, labelled ‘cal.’ here) cases in that paper, corresponding to whether nuisance parameters are marginalised over or fixed. For context, we show current measurements from current BAO and ‘full-shape’ measurements in the SDSS galaxy, quasar and Lyman-α datasets from the MGS, BOSS and eBOSS surveys (yellow) presented by Cuceu et al. (2023).

In the text

	Fig. A.1. Normalised full covariance matrix Ctot from the Flagship zeff = 1.0 redshift bin, estimated using jackknife resampling. Dashed lines differentiate between monopole, quadrupole and hexadecapole blocks of the data vector.
In the text

Fig. A.2. Redshift dependence of the measured cosmological parameters obtained from Flagship using different covariance prescriptions. The blue points show results for fits to the data vector for the perfect reconstruction case using Cdata and are the same as in Fig. 8. Black points are for fits to the same data vector but using the full covariance Ctot. Light red points are for the fit to the realistic reconstruction case when using Ctot. The use of the full covariance Ctot greatly reduces the statistical uncertainty but can introduce systematic offsets when combined with the use of reconstruction.

In the text

	Fig. A.3. Diagonal elements of the full covariance matrix Ctot (black lines), as well as the three different contributions from Eq. (A.1), showing the reduced variance that results from properly accounting for correlation between data and theory. The three different panels represent the monopole, quadrupole and hexadecapole part, respectively.
In the text