© The Authors 2023

1. Introduction

The late-time cosmic acceleration of the Universe on large scales is an established observational fact (Riess et al. 1998; Perlmutter et al. 1998). The cosmological standard model, the so-called ΛCDM model, explains the acceleration as the effect of a cosmological constant Λ, postulating the existence of a dark energy (DE) component with negative pressure, whose physical nature remains poorly understood. Minimal extensions to such a model include a variation of DE with time, parametrized by an equation of state following Chevallier-Polarski-Linder (CPL): w(z) = w₀ + z/(1 + z) w_a (Chevallier & Polarski 2001; Linder 2002). Such pressure impacts not only the background evolution of the Universe, but also the growth of large-scale structure (LSS): the greater the acceleration, the slower the structures can grow. While different DE models have been proposed (Ratra & Peebles 1988; Steinhardt et al. 1999; Armendariz-Picon et al. 2000; Copeland et al. 2006), the physical nature of DE remains unknown. Precise measurements of w(z) promise to shed light on its features.

Massive relic neutrinos can also affect the growth of structure and the expansion history of the Universe. Due to their free streaming with hot thermal velocities (Lesgourgues & Pastor 2006, 2012, 2014; Kiakotou et al. 2008), they alter the epoch of matter-radiation equality and suppress the growth of structure at mildly nonlinear and small scales. The Universe can thus be used as a laboratory to constrain the neutrino mass scale.

Upcoming galaxy surveys are designed to probe the DE equation of state w and the neutrino mass scale, with direct observations and measurements of the growth of structure and of distance scales. Cosmological probes used in these investigations include galaxy clustering (GC) and weak lensing (WL).

At the GC level, measurements of the galaxy correlation function from LSS surveys are used to constrain cosmological parameters, exploiting the sensitivity of the galaxy density fluctuations to the underlying dark matter (DM) density field. At the WL level, images of large ensembles of galaxies provide the so-called cosmic shear: tiny distortions in the shapes of galaxies due to the gravitational potential produced by intervening density perturbations, crossed by light propagating from the source to the observer. Images of galaxy shapes, complemented by an estimation of their redshifts, allow us to measure the structure growth and improve the inference on cosmological parameters.

Recently, aside from traditional GC, under-dense regions known as cosmic voids have been used to extract cosmological information (Lavaux & Wandelt 2010; Sutter et al. 2014; Hamaus et al. 2016, 2017, 2020; Nadathur et al. 2019; Pisani et al. 2019). Their sizes range from about ten to a few hundreds of Mpc (Sheth & Van De Weygaert 2004). Galaxy redshift surveys allow us to build cosmic void catalogs (Sutter et al. 2012; Micheletti et al. 2014; Clampitt & Jain 2015; Sánchez et al. 2017; Mao et al. 2017; Hawken et al. 2020; Pollina et al. 2019; Hamaus et al. 2020; Aubert et al. 2020; Nadathur et al. 2020). As for galaxies and galaxy clusters, it is possible to study the correlation function of cosmic voids and their possible cross-correlation with other cosmological probes (Granett et al. 2008, 2015; Ilić et al. 2013; Hamaus et al. 2014a). In particular, it is possible to measure the cross-correlation between WL and under-dense regions, and exploit this signal to infer cosmological information, and break possible parameter degeneracy (Krause et al. 2013; Melchior et al. 2014; Clampitt & Jain 2015; Sánchez et al. 2017; Fang et al. 2019).

The use of cosmic voids as a cosmological probe presents several advantages. Firstly, voids are less affected than DM halos by shell-crossing and virialization effects, which makes their dynamical evolution more amenable to theoretical models and easier to describe (Sheth & Van De Weygaert 2004; Hamaus et al. 2014b; Stopyra et al. 2020).

Secondly, recent works in the literature (e.g., Lavaux & Wandelt 2012; Pisani et al. 2015, 2019; Massara et al. 2015; Kreisch et al. 2019, 2021; Schuster et al. 2019; Verza et al. 2019; Bos et al. 2012; Lee & Park 2009) have shown that void formation and evolution are sensitive to the DE equation of state, and to the total neutrino mass, in a qualitatively and quantitatively different way with respect to over-dense structures, because of the different scales involved and their different nature.

For instance, the DE contribution within voids can be more dominant than in the Universe, on average, allowing such a probe to have a large sensitivity to its amount and evolution. In this work we consider DE as an effective quintessence “fluid” that does not cluster significantly and is nearly homogeneously distributed in space. Due to their hot thermal velocities, massive neutrinos have a free-streaming length that, depending on their mass and redshift, can range between one hundred to a few tens of h⁻¹ Mpc (Lesgourgues & Pastor 2006), thus matching the typical void sizes (Kreisch et al. 2019).

This paper presents the first forecasts on cosmological parameters based on the combination of WL, void-lensing cross-correlation, and angular void clustering, where voids are found in the galaxy photometric catalog of the Euclid survey, an ESA medium-class mission, currently scheduled for launch in 2023. The paper belongs to a series of companion papers investigating the scientific return that can be expected from voids from the Euclid mission (Hamaus et al. 2021; Contarini et al. 2022).

It is expected that the Euclid photometric galaxy catalog will be characterized by a surface density of n_g = 30 arcmin⁻² (Laureijs et al. 2011), and the imaging Euclid catalog will contain the shapes of about 1.5 billion galaxies, observed in the visible range with the VIS instrument (Cropper et al. 2016). The redshifts of such galaxies will be measured in photometric mode, using the Near Infrared Spectroscopic Photometric (NISP) instrument (Costille et al. 2019), complemented by ground-based observations in different bands.

The parameter forecasts presented here closely follow the recipe of the Euclid Inter Science Task-Force for Forecasts (IST:F) group (Euclid Collaboration 2020, hereafter EC20) in the case of WL and photometric galaxy clustering (GC_ph): the covariance between cosmological parameters is evaluated with the Fisher matrix approach, using angular power spectra computed within the Limber approximation (Limber 1953), and exploiting a tomography technique. Forecasts are given for spatially flat ΛCDM, νΛCDM, w₀w_aCDM, and νw₀w_aCDM model cosmologies, adopting the following set of cosmological parameters: the reduced Hubble constant h defined via H₀ = h 100 km s⁻¹ Mpc⁻¹, the baryon density parameter Ω_b at present time, the total matter density parameter Ω_m at present time, the sum of the three active neutrino masses M_ν = ∑_im_νi, the massive neutrino density parameter Ω_ν = M_ν (eV)/(93.14 h²) at present time, the cold DM density parameter Ω_c = Ω_m − Ω_b − Ω_ν at present time, where Ω_m = Ω_b + Ω_c + Ω_ν, the DE density parameter Ω_DE = 1 − Ω_m at present time, the parameters of the DE equation of state, w₀ and w_a, the rms of the matter linear density fluctuations, σ₈, inside a radius of 8 h⁻¹ Mpc, and the scalar spectral index n_s. The cosmological parameters varied in the analysis, together with their fiducial values, are summarized in Table 1.

The article is organized as follows. Section 2 illustrates the theoretical modeling of void clustering and bias. Section 3 describes the angular power spectra exploited as observables for parameter forecasts. Section 4 details the evaluation of the Fisher matrix. Section 5 reports the results and interpretation of our analysis. Finally, Sect. 6 presents our concluding remarks.

2. The clustering of cosmic voids

We evaluated the angular power spectrum, $C_{\ell }^{\mathrm{vv}}$, of cosmic voids from the void auto-power spectrum, P_vv(k, z), while the angular cross-power spectra, $C_{\ell }^{\gamma \mathrm{v}}$, between voids and WL were evaluated from the void-matter cross-spectrum P_vm(k, z). The latter are shown, together with the nonlinear matter power spectrum P_mm(k, z), in the top left panel of Fig. 1. P_vv(k, z) was obtained from the matter power spectrum, P_mm(k, z), assuming the Poisson shot noise (Hamaus et al. 2014a; Chan et al. 2014, 2019; Jamieson & Loverde 2019) and the “effective void bias” $b_{\text{v}}^{\text{eff}}(z)$ defined in Eq. (19). To this aim we adopted the following relations:

$$\begin{array}{cc}\hfill P_{\mathrm{vv}}(k,z)& =[1-{\mathrm{S}}_N(k)]b_{\mathrm{v}}^{\mathrm{eff}}{(z)}^2{P}_{\mathrm{mm}}(k,z)+1/{\overline{n}}^{\mathrm{v}}(z)\hfill \\ \hfill & \equiv b_{\mathrm{v}}^{\mathrm{eff}}{(z)}^2{\widehat{P}}_{\mathrm{mm}}(k,z)+1/{\overline{n}}^{\mathrm{v}}(z),\hfill \end{array}$$(1)

Fig. 1. Figure containing the power spectra used in this work. Top left: nonlinear matter auto-power spectrum Pmm (solid red line), the void auto-power spectrum Pvv without the shot noise (dashed blue line), computed assuming the pessimistic void bias scenario with bv eff ≈ −11.9, and the absolute value of the void-matter cross-power spectrum Pvm (dot-dashed black line), all computed at redshift z = 0.001. Top right: tomographic void angular auto-power spectra $C_{ij}^{\text{vv}}(\mathcal{l})$ for the diagonal (i = j) tomographic bins, without shot noise. Bottom left: tomographic lensing auto-power spectra $C_{ij}^{\gamma \gamma }(\mathcal{l})$ for the diagonal (i = j) tomographic bins, without shape noise. Bottom right: tomographic void-lensing angular cross-spectra $C_{ij}^{\text{v}\gamma }(\mathcal{l})$ for the diagonal (i = j) tomographic bins. All the spectra are theoretically evaluated in the reference cosmology reported in Table 1 and the bin endpoints are given in Eq. (24).

where the void bias, $b_{\text{v}}^{\text{eff}}(z)$, is assumed to be scale-independent¹, ${\overline{n}}^{\mathrm{v}}(z)$ is the void number density, and ${\widehat{P}}_{\mathrm{mm}}(k,z)$ is the nonlinear matter power spectrum when an additional wavenumber filter, S_N(k), is applied to cut small scales. The low-k pass filter, S_N(k), is necessary since in this work we focus on void clustering at large scales, as, in this case, the behavior of the void power spectrum, P_vv(k, z), is well described by linear theory via a simple multiplicative factor for the void bias. This can be understood in the context of the halo model formalism (Cooray & Sheth 2002), which has an analog formulation for voids (Voivodic et al. 2020). Such an assumption is inaccurate at small scales (excluded from the analyses in the present work), which are sensitive to the void density profile (Hamaus et al. 2014a). While the inclusion of small scales could improve the results presented in this work, this poses a serious modeling challenge, as to date in the literature there are no accurate models of the cosmological dependence of the void density profile. Furthermore, we point out that the high shot noise would likely reduce the contribution of the small scales. Therefore, we decided for a more conservative approach, considering only scales independent of the particular void profile, which we can confidently include in our analysis.

Analogously, the void-matter cross-spectrum P_vm(k, z) is defined as

$$\begin{array}{c}\hfill P_{\mathrm{vm}}(k,z)=[1-{\mathrm{S}}_N(k)]b_{\mathrm{v}}^{\mathrm{eff}}(z)P_{\mathrm{mm}}(k,z)\equiv b_{\mathrm{v}}^{\mathrm{eff}}(z){\widehat{P}}_{\mathrm{mm}}(k,z).\end{array}$$(2)

To avoid numerical instabilities, we prefer the low-k pass filter rather than a sharp cut in k. S_N(k) is the so-called smoothstep function, defined as

$$\begin{array}{c}\hfill {\mathrm{S}}_N(k)=\{\begin{array}{cc}0\hfill & \mathrm{if}\frac{k}{k_{\mathrm{ref}}}\le a\hfill \\ {\left(\frac{k}{k_{\mathrm{ref}}}\right)}^{N+1}{\sum }_{n=0}^N\left(\begin{array}{c}2N+1\\ N-n\end{array}\right)\left(\begin{array}{c}N+n\\ n\end{array}\right){(-\frac{k}{k_{\mathrm{ref}}})}^n\hfill & \mathrm{if}a\le \frac{k}{k_{\mathrm{ref}}}\le b\hfill \\ 1\hfill & \mathrm{if}\frac{k}{k_{\mathrm{ref}}}\ge b,\hfill \end{array}\end{array}$$(3)

where N measures the degree of smoothness of the function itself: the first discontinuous derivative of S_N(k) is the (N + 1)th derivative. The order used in S_N for producing P_vv(k, z) and P_vm(k, z) is N = 3. The smoothstep parameters are set to: a = 1, b = 1.8, and k_ref = 0.25 h Mpc⁻¹; the filter suppresses P_vv(k, z) and P_vm(k, z) for k > k_ref. The cutoff scale k_ref was determined combining the mean radius of voids ${\overline{r}}_{\mathrm{v}}$ in the catalog with the void exclusion principle, from which we know that the cutoff scale is given by $k\simeq 2\pi /{\overline{r}}_{\mathrm{v}}$ (Hamaus et al. 2014a). We chose this particular filter because it allowed us to better control its effect, since this function has some more parameters. Furthermore, we checked that the exact shape of the filter did not alter sensibly our results, as it was used to remove numerical artifacts.

Summarizing, we define the probe-dependent power spectrum ${\widehat{P}}_{\mathrm{AB}}$ as

$$\begin{array}{c}\hfill {\widehat{P}}_{\mathrm{AB}}(k,z)=\{\begin{array}{cc}{P}_{\mathrm{mm}}(k,z)\hfill & \mathrm{if}\mathrm{A}=\mathrm{B}=\gamma \hfill \\ [1-{\mathrm{S}}_N(k)]P_{\mathrm{mm}}(k,z)\hfill & \mathrm{else},\hfill \end{array}\end{array}$$(4)

where γ stands for cosmic shear.

As we describe below, to obtain the effective void bias, $b_{\text{v}}^{\text{eff}}(z)$, we used the peak-background split (PBS) formalism (Sheth & Van De Weygaert 2004), and weighted the void bias, b_v(z), over the void size function (Sheth & Van De Weygaert 2004; Jennings et al. 2013), that is to say, the comoving number density of voids per radius interval. The void size function was predicted from the “excursion set formalism” in the form of the Sheth & van de Weygart (SvdW) model (Sheth & Van De Weygaert 2004), according to which voids are treated as isolated objects, whose evolution is, therefore, not affected by the environment and is assumed to be spherically symmetric.

Under such assumptions, the prediction of the void abundance is completely described by two physical parameters. The former is the density threshold for structure collapse, δ_c, linearly extrapolated to the present time as δ_c(z) = δ_c(0)/G(z), where G(z) is the so-called linear growth factor normalized at z = 0, and δ_c(0) slightly depends on cosmology, here assumed to be δ_c(0) = 1.686. The latter is the density threshold for the void formation δ_v, linearly extrapolated to the present time as δ_v(z) = δ_v(0)/G(z), with δ_v(0) = − 0.9 (see, e.g., Chan et al. 2014; Ronconi & Marulli 2017; Contarini et al. 2019; Verza et al. 2019, recently showing that various values for the threshold δ_v(0) can be used to model voids more reliably). We chose a value of δ_v(0) by comparing our effective bias prediction against the measurement from the Euclid Flagship simulation (more on this in Sect. 3).

In the void abundance characterization, both δ_v and δ_c have an important role. The void formation threshold, δ_v, is the value an under-dense region needs to overcome to turn into a void. The density threshold, δ_c, is necessary to account for the “void-in-cloud problem”, where an under-dense region inside an over-dense one will not become a void, since the latter will collapse.

The SvdW model predicts the void abundance to be a function of both δ_c and δ_v, via the so-called void-in-cloud parameter, 𝒟, and the dimensionless parameter, x, defined as

$$\begin{array}{c}\hfill \mathcal{D}=\frac{\left|{\delta }_{\mathrm{v}}\right|}{{\delta }_{\mathrm{c}}+\left|{\delta }_{\mathrm{v}}\right|},x=\frac{\mathcal{D}}{\left|{\delta }_{\mathrm{v}}\right|}\sigma ,\end{array}$$(5)

where σ² is the variance of the filtered linear density field on a scale R

$$\begin{array}{c}\hfill {\sigma }^2(R)\equiv S(R)={\displaystyle \int }\frac{\mathrm{d}k}{k}\frac{k^3{P}_{\mathrm{lin}}(k)}{2{\pi }^2}{|W(kR)|}^2.\end{array}$$(6)

Here P_lin(k) is the linear matter power spectrum at z = 0, and W(kR) is the top-hat filter in Fourier space:

$$\begin{array}{c}\hfill W(kR)=\{\begin{array}{cc}1\hfill & kR\le 1\hfill \\ 0\hfill & kR>1.\hfill \end{array}\end{array}$$(7)

The abundance of voids with mass M is predicted by the SvdW model as

$$\begin{array}{c}\hfill \frac{\mathrm{d}n}{\mathrm{d}lnM}=\frac{{\rho }_{\mathrm{m}}}{M}{f}_{ln\sigma }(\sigma )\frac{\mathrm{d}ln{\sigma }^{-1}}{\mathrm{d}lnM},\end{array}$$(8)

where ρ_m is the mean background density of matter and the multiplicity function, f_lnσ, represents the fraction of voids in a unit range of lnM, and it is shown to be (Sheth & Van De Weygaert 2004)

$$\begin{array}{c}\hfill f_{ln\sigma }(\sigma )=2{\displaystyle \sum _{j=1}^{\infty }exp(-\frac{{(j\pi x)}^2}{2})j\pi x^{2}sin(j\pi \mathcal{D}).}\end{array}$$(9)

Dealing with observations, it is more convenient to work in terms of void volumes rather than void masses. In linear theory, the relation between the void volume and its mass is

$$\begin{array}{c}\hfill \frac{{\rho }_{\mathrm{m}}}{M(r_{\mathrm{L}})}=\frac{1}{V(r_{\mathrm{L}})},\end{array}$$(10)

where $V(r_{\text{L}})\equiv 4\pi r_{\text{L}}^3/3$ is the volume of a sphere with radius r_L computed in the linear regime.

Even if the SvdW approach accurately predicts the void mass abundance, the predicted volume fraction of the Universe inside voids is larger than one. In order to fix this issue, the SvdW approach has been extended to the V dn size function from Jennings et al. (2013), who noticed that such a volume fraction issue originates from the incorrect matching between linear and corresponding nonlinear quantities. Since this approach correctly predicts the volume fraction of the Universe inside voids, it is known as the volume-conserving void size function.

Jennings et al. (2013) corrected this match by requiring the void volume fraction of the Universe in the nonlinear regime to have the same value as in the linear one:

$$\begin{array}{c}\hfill V(r_{\mathrm{v}})\mathrm{d}n={V\left(r_{\mathrm{L}}\right)\mathrm{d}{n}_{\mathrm{L}}|}_{r_{\mathrm{L}}=r_{\mathrm{L}}(r_{\mathrm{v}})},\end{array}$$(11)

where the subscript L denotes quantities computed in the linear regime. The V dn void size function then reads (Jennings et al. 2013)

$$\begin{array}{c}\hfill \frac{\mathrm{d}n}{\mathrm{d}ln{r}_{\mathrm{v}}}={\frac{f_{ln\sigma }(\sigma )}{V(r_{\mathrm{v}})}\frac{\mathrm{d}ln{\sigma }^{-1}}{\mathrm{d}ln{r}_{\mathrm{L}}}\frac{\mathrm{d}ln{r}_{\mathrm{L}}}{\mathrm{d}ln{r}_{\mathrm{v}}}|}_{r_{\mathrm{L}}=r_{\mathrm{L}}(r_{\mathrm{v}})},\end{array}$$(12)

where r_v indicates the actual void radius. The conversion from r_L to r_v is given by

$$\begin{array}{c}\hfill \frac{r_{\mathrm{v}}}{r_{\mathrm{L}}}={\left(\frac{{\rho }_{\mathrm{m}}}{{\rho }_{\mathrm{v}}}\right)}^{1/3},\end{array}$$(13)

where ρ_v is the average density within the void.

However, the expression of the multiplicity function, f_lnσ, in Eq. (9) is difficult to evaluate in practice. Therefore, Jennings et al. (2013) suggested the following approximation that works at 0.2% accuracy:

$$\begin{array}{c}\hfill f_{ln\sigma }(\sigma )\approx \{\begin{array}{cc}\sqrt{\frac{2}{\pi }}\frac{\left|{\delta }_{\mathrm{v}}\right|}{\sigma }exp(-\frac{{\delta }_{\mathrm{v}}^2}{2{\sigma }^2})\hfill & x\le 0.276\hfill \\ 2{\sum }_{j=1}^{4}exp(-\frac{{(j\pi x)}^2}{2})j\pi x^{2}sin(j\pi \mathcal{D})\hfill & x>0.276.\hfill \end{array}\end{array}$$(14)

Finally, the V dn size function is given by

$$\begin{array}{c}\hfill \frac{\mathrm{d}n}{\mathrm{d}ln{r}_{\mathrm{v}}}={\frac{f_{ln\sigma }(\sigma )}{V(r_{\mathrm{v}})}\frac{\mathrm{d}ln{\sigma }^{-1}}{\mathrm{d}ln{r}_{\mathrm{L}}}|}_{r_{\mathrm{L}}=r_{\mathrm{v}}/1.7},\end{array}$$(15)

and, after some rearrangements, its final expression is

$$\begin{array}{c}\hfill \frac{\mathrm{d}n}{\mathrm{d}{r}_{\mathrm{v}}}=\frac{f_{ln\sigma }(\sigma )}{r_{\mathrm{v}}V(r_{\mathrm{v}})}{[-\frac{r_{\mathrm{L}}}{\sigma }\frac{\mathrm{d}\sigma }{\mathrm{d}{r}_{\mathrm{L}}}]}_{r_{\mathrm{L}}=r_{\mathrm{v}}/1.7},\end{array}$$(16)

where r_v = 1.7 r_L is predicted by considering a spherical expansion.

As mentioned previously, the void bias, b_v, at a fixed void radius, r_v, is computed using the PBS formalism (Sheth & Van De Weygaert 2004; Chan et al. 2014):

$$\begin{array}{c}\hfill b_{\mathrm{v}}(z,r_{\mathrm{v}})=1+\frac{{\nu }^2-1}{{\delta }_{\mathrm{v}}}+\frac{{\delta }_{\mathrm{v}}\mathcal{D}}{4{\delta }_{\mathrm{c}}^2{\nu }^2},\end{array}$$(17)

where the z dependence is encapsulated in δ_c and δ_v, and the r_v dependence in ν:

$$\begin{array}{c}\hfill \nu \equiv \frac{|{\delta }_{\mathrm{v}}|}{\sigma (r_{\mathrm{v}})}·\end{array}$$(18)

We defined the effective void bias by weighting Eq. (17) with the void size function V dn in Eq. (16):

$$\begin{array}{c}\hfill b_{\mathrm{v}}^{\mathrm{eff}}(z)=\frac{{\int }_{r_{\mathrm{v},\min }}^{r_{\mathrm{v},\max }}\mathrm{d}{r}_{\mathrm{v}}\frac{\mathrm{d}n}{\mathrm{d}{r}_{\mathrm{v}}}{b}_{\mathrm{v}}(z,r_{\mathrm{v}})}{{\int }_{r_{\mathrm{v},\min }}^{r_{\mathrm{v},\max }}\mathrm{d}{r}_{\mathrm{v}}\frac{\mathrm{d}n}{\mathrm{d}{r}_{\mathrm{v}}}},\end{array}$$(19)

exploiting in this way all void radii r ≥ r_v, minh⁻¹ Mpc. In this work, we assumed r_v, min = 25 h⁻¹ Mpc, and we provide forecasts treating the evolution of the void bias in two main scenarios. A “pessimistic void bias” scenario, where the redshift evolution of the void bias is supposed to be known from linear theory, but its absolute normalization is unknown. In this case, $b_{\text{v}}^{\text{eff}}$ is evolved with the growth factor of the fiducial cosmology. Its fiducial value today is assumed to be $b_{\text{v}}^{\text{eff}}(z=0)=-11.9$, according to Eq. (19), and it is marginalized over as a nuisance parameter. This void bias scenario will be combined with the pessimistic WL scenario as described in Sect. 3. An “optimistic void bias” scenario, where we exploit the cosmology dependence of both the void size function and void bias. In particular, such dependence is encapsulated in the growth factor evolution of the density thresholds δ_c and δ_v, which appear both in the void bias definition and in the void size function. This configuration is considered to be optimistic, since we assume that the value of b_v, entering $b_{\text{v}}^{\text{eff}}$, is known from linear theory, and therefore its dependence on the cosmological parameters is exploitable. This void bias scenario will be combined with the optimistic WL scenario as described in Sect. 3. To summarize, we present forecasts in two scenarios: a pessimistic one, where we treat $b_{\text{v}}^{\text{eff}}(z=0)$ as an unknown quantity to be marginalized over, fixing its evolution to the fiducial cosmology; and an optimistic one, where we assume we are able to model the cosmological dependence of $b_{\text{v}}^{\text{eff}}(z)$ and use it as a further source of information on cosmological parameters. Providing forecasts in both cases is, in our opinion, the fairest approach, since we present parameter constraints spanning from the worst to the best case scenario.

3. Observables: Angular power spectra

This section summarizes the observables used for cosmological parameter forecasts based on WL, voids, and their cross-correlation. Here we describe in particular the angular power spectra C(ℓ), which are functions of the multipole number ℓ, and are defined as the spherical harmonic transform of a two-point correlation function. The a_ℓm, the coefficients of the spherical harmonics decomposition of a field ℱ on the sphere (e.g., the cosmic shear or void density fields), are defined on the sky as

$$\begin{array}{c}\hfill a_{\ell m}={\displaystyle \int }\mathrm{d}\mathrm{\Omega }\mathcal{F}(\theta ,\varphi )Y_{\ell m}(\theta ,\varphi ).\end{array}$$(20)

The C(ℓ) represent the two-point correlation function of the a_ℓm and they are diagonal in ℓ and independent of m because of statistical isotropy.

In the forecasts presented in this work, we used three kinds of C(ℓ): the void-void auto-correlation C^vv(ℓ), the lensing-lensing auto-correlation C^γγ(ℓ), and the void-lensing cross-correlation C^vγ(ℓ). We computed the C(ℓ) “tomographically” (Hu 1999) in a set of redshift bins, that is to say they were evaluated over two redshift bins i and j, and were generically denoted as $C_{ij}^{\text{AB}}(\mathcal{l})$. For example, $C_{ij}^{\text{v}\gamma }(\mathcal{l})$ is the spherical harmonic transform of the correlation function between the lensing signal in the ith bin and the void signal in the jth bin.

We evaluated C^vv(ℓ) and C^vγ(ℓ) for ℓ ∈ [10, 1500], while for C^γγ(ℓ), two scenarios were considered: a “pessimistic” WL scenario, with ℓ ∈ [10, 1500], and an “optimistic” WL scenario, with ℓ ∈ [10, 5000].

The Limber approximation (Limber 1953; LoVerde & Afshordi 2008) allowed us to write a simple integral expression for the tomographic C(ℓ), where the power spectrum ${\widehat{P}}_{\mathrm{AB}}(k,z)$, defined in Eq. (4), enters the integral as

$$\begin{array}{c}\hfill C_{\mathit{ij}}^{\mathrm{AB}}(\ell )\simeq \frac{c}{H_0}{\displaystyle {\int }_{z_{\min }}^{z_{\max }}\mathrm{d}z\frac{W_i^{\mathrm{A}}(z)W_j^{\mathrm{B}}(z)}{E(z)r^2(z)}{\widehat{P}}_{\mathrm{AB}}[\frac{\ell +1/2}{r(z)},z],}\end{array}$$(21)

where $W_i^{\mathrm{A}/\mathrm{B}}(z)$ are suitable “weight” functions for the different probes A and B defined in Eqs. (27) and (30), E(z) is the dimensionless Hubble parameter,

$$\begin{array}{cc}\hfill E(z)& =\frac{H(z)}{H_0}\hfill \\ \hfill & =\sqrt{{\mathrm{\Omega }}_{\mathrm{m}}{(1+z)}^3+(1-{\mathrm{\Omega }}_{\mathrm{m}}){(1+z)}^{3(1+w_0+w_a)}{\mathrm{e}}^{-3w_{a}z/(1+z)}},\hfill \end{array}$$(22)

and r(z) is the comoving distance,

$$\begin{array}{c}\hfill r(z)=\frac{c}{H_0}{\displaystyle {\int }_0^z\frac{\mathrm{d}{z}^{\prime }}{E(z^{\prime })}·}\end{array}$$(23)

Throughout the forecast, the integral expression of $C_{ij}^{\text{AB}}(\mathcal{l})$ was computed numerically with the Simpson method (Brew & Atkinson 1979). The signals $C_{ij}^{\text{vv}}(\mathcal{l})$, $C_{ij}^{\gamma \gamma }(\mathcal{l})$, and $C_{ij}^{\text{v}\gamma }(\mathcal{l})$, evaluated using the reference cosmology in Table 1, are shown in Fig. 1.

In the Euclid photometric survey, the minimum and the maximum redshifts are z_min = 0.001 and z_max = 2.5, respectively. The C(ℓ) and the distributions of galaxies and voids are sampled in ten redshift bins i, with the following endpoints:

$$\begin{array}{cc}\hfill \{z_{\mathrm{e}}\}=\{& 0.001,0.418,0.560,0.678,0.789,0.900,1.019,1.155,\hfill \\ \hfill & 1.324,1.576,2.50\}.\hfill \end{array}$$(24)

These bins were chosen to be equi-populated in the photometric galaxy catalog. Although, as explained in Sect. 3, the mock void catalog was split in ten “equi-spaced” redshift bins, we preferred to convert the redshift void distribution to the same binning adopted for galaxies.

In order to simplify the notation of the integrand in Eq. (21), we defined the “kernel” functions as

$$\begin{array}{c}\hfill K_{\mathit{ij}}^{\mathrm{AB}}(z)=\frac{c}{H_0}\frac{W_i^{\mathrm{A}}(z)W_j^{\mathrm{B}}(z)}{E(z)r^2(z)},\end{array}$$(25)

with i, j being the tomographic indices, and A, B the considered probes.

The tomographic void-void angular power spectra $C_{ij}^{\text{vv}}(\mathcal{l})$ were computed as

$$\begin{array}{c}\hfill C_{\mathit{ij}}^{\mathrm{vv}}(\ell )={\displaystyle {\int }_{z_{\min }}^{z_{\max }}\mathrm{d}z{K}_{\mathit{ij}}^{\mathrm{vv}}(z){\widehat{P}}_{\mathrm{mm}}[\frac{\ell +1/2}{r(z)},z],}\end{array}$$(26)

where ${\widehat{P}}_{\mathrm{mm}}(k,z)$ is defined in Eq. (1) and we omitted the shot-noise additive term $1/{\overline{n}}_i^{\mathrm{v}}$ defined in Eq. (40).

The void weight function is defined as

$$\begin{array}{c}\hfill W_i^{\mathrm{v}}(z)=\frac{H(z)}{c}{n}_i^{\mathrm{v}}(z)b_{\mathrm{v},i}^{\mathrm{eff}}.\end{array}$$(27)

The void effective bias, b_veff, i, entering in the void weight function, is given by b_veff(z), as defined in Eq. (19), evaluated in the center of the ith tomographic bin. The un-normalized projected void density distribution in redshift was measured directly from the Euclid Flagship mock photometric galaxy catalog², as explained in Sect. 3.

The tomographic void-lensing angular power spectra, $C_{ij}^{\text{vv}}(\mathcal{l})$, and the lensing-lensing angular power spectra, $C_{ij}^{\gamma \gamma }(\mathcal{l})$, were computed, respectively, as

$$\begin{array}{c}\hfill C_{\mathit{ij}}^{\mathrm{v}\gamma }(\ell )={\displaystyle {\int }_{z_{\min }}^{z_{\max }}\mathrm{d}z{K}_{\mathit{ij}}^{\mathrm{v}\gamma }(z){\widehat{P}}_{\mathrm{mm}}[\frac{\ell +1/2}{r(z)},z]}\end{array}$$(28)

and

$$\begin{array}{c}\hfill C_{\mathit{ij}}^{\gamma \gamma }(\ell )={\displaystyle {\int }_{z_{\min }}^{z_{\max }}\mathrm{d}z{K}_{\mathit{ij}}^{\gamma \gamma }(z)P_{\mathrm{mm}}[\frac{\ell +1/2}{r(z)},z],}\end{array}$$(29)

where in $C_{ij}^{\gamma \gamma }(\mathcal{l})$ we omitted the shape-noise term given by Eq. (40). Following EC20, in the kernels of Eqs. (28) and (29), the lensing weight function reads

$$\begin{array}{c}\hfill W_i^{\gamma }(z)=\frac{3}{2}{\left(\frac{H_0}{c}\right)}^2{\mathrm{\Omega }}_{\mathrm{m}}(1+z)r(z){\stackrel{\sim }{W}}_i^{\gamma }(z),\end{array}$$(30)

with

$$\begin{array}{c}\hfill {\stackrel{\sim }{W}}_i^{\gamma }(z)={\displaystyle {\int }_z^{z_{max}}\mathrm{d}{z}^{\prime }{n}_i^{\mathrm{g}}\left(z^{\prime }\right)[1-\frac{r(z)}{r\left(z^{\prime }\right)}]·}\end{array}$$(31)

Here ${\stackrel{\sim }{W}}_i^{\gamma }(z)$ is the lensing efficiency, defined as an integral expression of the observed galaxy distribution $n_i^{\text{g}}(z)$.

The density distribution of the “observed” galaxies in the ith tomographic bin, $n_i^{\text{g}}(z)$, was computed as a convolution of the galaxy distribution, n^g(z), and the photometric instrument response, as described below. According to EC20, the galaxy redshift distribution, n^g(z), adopted in this work is

$$\begin{array}{c}\hfill n^{\mathrm{g}}(z)\propto {\left(\frac{z}{z_0}\right)}^{2}exp[-{\left(\frac{z}{z_0}\right)}^{3/2}],\end{array}$$(32)

where $z_0=0.9/\sqrt{2}$ (Laureijs et al. 2011). Redshift measurements are affected by experimental and reconstruction effects. We modeled the probability p_ph(z_p|z) that a true galaxy redshift z will be reconstructed with a photometric redshift z_p as

$$\begin{array}{cc}\hfill p_{\mathrm{ph}}(z_{\mathrm{p}}|z)=& \frac{1-f_{\mathrm{out}}}{\sqrt{2\pi }{\sigma }_{\mathrm{b}}(1+z)}exp\{-\frac{1}{2}{\left[\frac{z-c_{\mathrm{b}}{z}_{\mathrm{p}}-z_{\mathrm{b}}}{{\sigma }_{\mathrm{b}}(1+z)}\right]}^2\}\hfill \\ \hfill & +\frac{f_{\mathrm{out}}}{\sqrt{2\pi }{\sigma }_{\mathrm{o}}(1+z)}exp\{-\frac{1}{2}{\left[\frac{z-c_{\mathrm{o}}{z}_{\mathrm{p}}-z_{\mathrm{o}}}{{\sigma }_{\mathrm{o}}(1+z)}\right]}^2\}·\hfill \end{array}$$(33)

The modeling and the parameters entering p_ph(z_p|z) are reported in Table 2, and were taken from EC20. The convolution of n^g(z) and p_ph(z_p|z) estimates the galaxy redshift distribution as measured by NISP in the photometric mode, together with ground-based observation. The density distribution, $n_i^{\text{g}}(z)$, of the “observed” galaxies in the ith tomographic bin is obtained convolving p_ph(z_p|z) with n^g(z) in the ith bin

$$\begin{array}{c}\hfill n_i^{\mathrm{g}}(z)=\frac{{\int }_{z_i^{-}}^{z_i^{+}}\mathrm{d}{z}_{\mathrm{p}}{n}^{\mathrm{g}}(z)p_{\mathrm{ph}}(z_{\mathrm{p}}|z)}{{\int }_{z_{\min }}^{z_{\max }}\mathrm{d}z{\int }_{z_i^{-}}^{z_i^{+}}\mathrm{d}{z}_{\mathrm{p}}{n}^{\mathrm{g}}(z)p_{\mathrm{ph}}(z_{\mathrm{p}}|z)},\end{array}$$(34)

where $z_i^{+}$ and $z_i^{-}$ are the edges of the ith bin as defined in Eq. (24). The density distributions, $n_i^{\text{g}}(z)$, in the different tomographic bins are shown in the left panel of Fig. 2.

Fig. 2. Galaxy and void distributions used in this work. Left: dimensionless normalized projected galaxy distributions $n_i^{\text{g}}(z)$ over the tomographic bins i. Right: dimensionless normalized projected void distribution in equi-spaced redshift bins. The corresponding un-normalized projected galaxy and void distributions are reported in Table 3 in units of arcmin−2.

The matter power spectrum, P_mm(k, z), was computed numerically using as Boltzmann solver the “Code for Anisotropies in the Microwave Background”³ (CAMB; Lewis et al. 1999)⁴. In order to obtain the cosmological parameter forecasts presented in this work, we computed the matter power spectra, P_mm(k, z), for each set of input cosmological parameters, which were varied with respect to the fiducial cosmology in Table 1 to also compute the derivatives entering the Fisher matrix approach (Fisher 1935). CAMB evaluates P_mm(k, z) over a k − z grid. According to EC20, the range of redshifts chosen in the present analysis is z ∈ [0.001,  2.5], the range of wave numbers is k ∈ [10⁻⁵,  400] Mpc⁻¹, and at least 300 equispaced points in redshift are used. CAMB performs a logarithmic k-binning. In the forecast, at least 50 steps per decade are used. The impact of different k and z binnings on the final results is reported in Sect. 5.6.

Since both $n_i^{\text{g}}(z)$ and $n_i^{\text{v}}(z)$ are normalized quantities, the variation of the cosmological parameters in the Fisher forecast computation does not (or negligibly) affect these quantities. In fact, the main effect on number counts from changing the model cosmology is due to the associated volume change, which cancels out for normalized quantities. Instead, the shot-noise matrix in Eq. (40) is un-normalized and computed at the fiducial cosmology assumed in the Fisher forecast computation. The cosmological parameter variation is needed for the computations of the C(ℓ) derivatives, and enters the matter power spectra, P_mm(k, z), the growth factor, G(z), the Hubble parameter, H(z), and the comoving distance, r(z), where the latter two terms change the weight functions $W_i^{\text{A}}$.

In order to use realistic estimates of void bias and void distribution from the Euclid photometric galaxy sample, in the present analysis we employed the Flagship mock galaxy catalog (Euclid Collaboration, in prep.). This is based on an N-body simulation of 12 600³ DM particles in a periodic box of 3780 h⁻¹ Mpc on a side (Potter et al. 2017), with a flat ΛCDM cosmology very similar to that reported in Table 1, namely Ω_m = 0.319, Ω_b = 0.049, Ω_Λ = 0.681, σ₈ = 0.83, and n_s = 0.96, h = 0.67, as obtained by Planck (Planck Collaboration XIII 2016). Dark matter halos were identified with the ROCKSTAR halo finder (Behroozi et al. 2013), and populated with central and satellite galaxies using a halo occupation distribution (HOD) framework to reproduce the relevant observables for the Euclid main cosmological probes. The HOD algorithm (Carretero et al. 2015; Crocce et al. 2015) was calibrated exploiting several observational constraints, such as the local luminosity function for the faintest galaxies (Blanton et al. 2003, 2005) and GC statistics as a function of luminosity and color (Zehavi et al. 2011). The resulting Flagship galaxy mock lightcone spanned one octant of the sky and simulated both spectroscopic and photometric Euclid galaxy samples. In this paper we consider the latter, which extends up to redshift z = 2.3 and in which a Gaussian photometric redshift error of Δz = 0.05(1 + z) was applied to each galaxy.

To identify cosmic voids in this catalog, we applied the 2D void finder of Sánchez et al. (2017), Vielzeuf et al. (2019) to the photometric galaxy sample, computing the number of voids in ten “equi-spaced” redshift bins, as shown in the right panel of Fig. 2, where the ratio between the void number in each bin and the bin width represents the void projected density at the bin center. We find that the 2D void population traced by photometric galaxies in the Flagship catalog extends from r_v, min ∼ 25 h⁻¹ Mpc up to r_v, max ∼ 300 h⁻¹ Mpc, and its projected spatial density, n^v(z), is obtained by interpolating the void density in each bin center.

Moreover, to verify that the bias modeling, Eq. (19), used in our forecasts was representative of Euclid observations, we measured the void bias redshift evolution in the obtained Flagship photometric void catalog. To this aim, we followed the methodology presented in Hamaus et al. (2014a). We used the open-source code nbodykit (Hand et al. 2018) to compute the void auto power-spectrum in eight redshift bins (z = [0.3, 0.5, 0.75, 1.0, 1.25, 1.5, 1.75, 2.]) with a bin size of Δz = 0.2. Then, from Eq. (1), we inferred the void bias as

$$\begin{array}{c}\hfill b_{\mathrm{v}}^{\mathrm{eff}}(z)=\sqrt{\frac{P_{\mathrm{vv}}(k,z)-1/{\overline{n}}^{\mathrm{v}}}{P_{\mathrm{mm}}(k,z)}},\end{array}$$(35)

where the void power spectrum was measured in the Flagship photometric void catalog. Then, we averaged the measured bias in the range 0.05 < k(h⁻¹ Mpc) < 0.1, that it to say, in the regime in which the void bias is constant as a function of the scale (see Hamaus et al. 2014a).

The evolution of the void bias as a function of redshift in the Flagship photometric void catalog (green error bars) is shown in Fig. 3. The measurement errors were estimated using a jackknife resampling. The red stars in Fig. 3 show the galaxy bias measured from the Flagship mock catalog using Eq. (35). Finally the solid blue line shows the theoretical bias modeling, Eq. (19), obtained via the PBS formalism and adopted for our forecast computation. We note that the measured negative void bias shows a redshift evolution close enough to the theoretical one. Therefore, we consider the PBS prescription of the void bias evolution acceptable for the purposes of the forecast analysis presented in this work. We accounted for any theoretical uncertainty on the void threshold and size function by marginalizing over the normalization of the effective bias in the pessimistic case. We verified that the agreement between the theoretical void size function and the one measured from the photometric void sample in the Flagship simulation was good enough for the forecasting purposes of this work, and we postpone to a future work a detailed calibration against simulated data, when more realistic mocks, accounting for detailed Euclid survey specification and systematics, will be available.

Fig. 3. Effective void bias (green error bars) and galaxy bias (red error bars) as measured in the Flagship simulation versus the theoretical effective void bias (blue line) of Eq. (19). Void bias measurements have a bigger error bar, compared to galaxy bias measurements, since there are fewer voids than galaxies.

4. The void-lensing Fisher matrix

The Fisher matrix formalism (Fisher 1935) is used to predict uncertainties of cosmological parameter measurements from different probes. This section summarizes the approach and provides details on its application to the cases considered in this work.

The Fisher matrix is defined as the expectation value of the second derivatives of the logarithm of the likelihood L:

$$\begin{array}{c}\hfill F_{\alpha \beta }=-\langle \frac{{\partial }^{2}lnL}{\partial \alpha \partial \beta }\rangle ,\end{array}$$(36)

where α and β are the parameters of interest. The expected “error covariance matrix” is the inverse of the Fisher matrix:

$$\begin{array}{c}\hfill {\mathcal{C}}_{\alpha \beta }={(F^{-1})}_{\alpha \beta }.\end{array}$$(37)

The diagonal elements of the error covariance matrix are the squares of the marginalized 1-σ errors on the parameters:

$$\begin{array}{c}\hfill {\sigma }_{\alpha }=\sqrt{{\mathcal{C}}_{\alpha \alpha }}.\end{array}$$(38)

As we deal with angular power spectra, it is convenient to introduce the matrix ${\mathrm{\Sigma }}_{\mathit{ij}}^{\mathrm{AB}}$ associated with a given C(ℓ):

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\mathit{ij}}^{\mathrm{AB}}(\ell )=\sqrt{\frac{2}{(2\ell +1)\mathrm{\Delta }\ell f_{\mathrm{sky}}}}[C_{\mathit{ij}}^{\mathrm{AB}}(\ell )+N_{\mathit{ij}}^{\mathrm{AB}}(\ell )],\end{array}$$(39)

where Δℓ is the multipole bin width, and f_sky the sky fraction covered by the survey. The “shot-noise matrix”, $N_{\mathit{ij}}^{\mathrm{AB}}(\ell )$, depends on the particular probe combination. For the C(ℓ) of the void clustering, WL, and void-lensing, assuming Poisson statistics, they are respectively

$$\begin{array}{c}\hfill N_{\mathit{ij}}^{\mathrm{vv}}=\frac{1}{{\overline{n}}_i^{\mathrm{v}}}{\delta }_{\mathit{ij}},N_{\mathit{ij}}^{\gamma \gamma }=\frac{{\sigma }_{ϵ}^2}{{\overline{n}}_i^{\mathrm{g}}}{\delta }_{\mathit{ij}},N_{\mathit{ij}}^{\gamma \mathrm{v}}=N_{\mathit{ij}}^{\mathrm{v}\gamma }=0,\end{array}$$(40)

where δ_ij is the Kronecker delta, ${\sigma }_{ϵ}^2$ is the galaxy shape noise, and ${\overline{n}}_i^{\mathrm{g}}$ and ${\overline{n}}_i^{\mathrm{v}}$ are the un-normalized average galaxy and void surface densities in the ith tomographic bin computed in the fiducial cosmology, respectively. The survey specifications used to compute the covariances are shown in Table 3.

In the case of a single probe (A = B), the covariance matrix of the a_ℓm is simply given by Eq. (39). When two or more probes are combined together, one needs to construct the full covariance matrix, Σ^XC, composed of matrix blocks defined in Eq. (39). For the probes γ and v considered in this work, Σ^XC is given by

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}^{\mathrm{XC}}(\ell )=\left(\begin{array}{cc}{\mathrm{\Sigma }}^{\gamma \gamma }(\ell )& {\mathrm{\Sigma }}^{\gamma \mathrm{v}}(\ell )\\ {\mathrm{\Sigma }}^{\mathrm{v}\gamma }(\ell )& {\mathrm{\Sigma }}^{\mathrm{vv}}(\ell )\end{array}\right).\end{array}$$(41)

The covariance and the C(ℓ) matrix share the same structure, and therefore the block matrix reads

$$\begin{array}{c}\hfill {\mathbf{C}}^{\mathrm{XC}}(\ell )=\left(\begin{array}{cc}{\mathbf{C}}^{\gamma \gamma }(\ell )& {\mathbf{C}}^{\gamma \mathrm{v}}(\ell )\\ {\mathbf{C}}^{\mathrm{v}\gamma }(\ell )& {\mathbf{C}}^{\mathrm{vv}}(\ell )\end{array}\right).\end{array}$$(42)

Assuming that the a_ℓm, defined in Eq. (20), follow a multivariate Gaussian distribution, an analytical expression for the Fisher matrix elements is given by⁵

$$\begin{array}{c}\hfill F_{\alpha \beta }={\displaystyle \sum _{\ell ={\ell }_{\min }}^{{\ell }_{\max }}\text{Tr}\left\{{[\mathrm{\Sigma }(\ell )]}^{-1}\frac{\partial \mathbf{C}(\ell )}{\partial \alpha }{[\mathrm{\Sigma }(\ell )]}^{-1}\frac{\partial \mathbf{C}(\ell )}{\partial \beta }\right\}.}\end{array}$$(43)

The formula above applies both to the single- and two-probe correlation case, provided that the C(ℓ) and the covariance matrices are chosen accordingly. The Fisher matrix computation involves the derivatives of the C(ℓ) with respect to cosmological and nuisance parameters. In order to ensure reliable results, numerical derivatives were computed with⁶ a numerical approach based on the SteM fitting procedure (Camera et al. 2017), which is based on a iterative linear regression, and a semi-analytical approach, in which both analytical and numerical derivatives are used.

5. Results: Correlations, errors, and the inclusion of the void bias evolution

In this section we present the results obtained implementing the analysis described in previous sections. We discuss the forecasts of cosmological parameters, expressed as Fisher matrix marginalized contours, obtained combining the expected Euclid WL and angular void clustering.

We report parameter forecasts in two different scenarios, namely the pessimistic and optimistic setups, described in Sect. 2 for the void bias, and in Sect. 3 for the WL, respectively. They are also summarized in Table 4. The 1-σ errors, together with the DE figure of merit⁷ (FoM), in both scenarios, and different cosmological models, are reported in Table 5. The corresponding contour plots are shown in Figs. 4–8. Focusing on the DE equation of state and on the total neutrino mass, the contour plots for w₀, w_a, M_ν, evaluated respectively for the pessimistic and optimistic scenarios, are also shown in the top panels of Fig. 9.

Fig. 4. Fisher matrix marginalized contours for the (baseline) ΛCDM model, in the pessimistic (top) and in the optimistic (bottom) scenarios.

Fig. 5. Fisher matrix marginalized contours for the νΛCDM model, in the pessimistic (top) and in the optimistic (bottom) scenarios.

Fig. 6. Fisher matrix marginalized contours for the w0waCDM model, in the pessimistic (top) and in the optimistic (bottom) scenarios.

Fig. 7. Full contour plots when GCsp is combined with WL + V + XC for the w0waCDM model, in the pessimistic (top) and in the optimistic (bottom) scenarios.

Fig. 8. Fisher matrix marginalized contours for the νw0waCDM model, in the pessimistic (top) and in the optimistic (bottom) scenarios.

5.1. The standard ΛCDM cosmology

The contour plots for the ΛCDM cosmology are reported in Fig. 4, the pessimistic case in the top panel and the optimistic case in the bottom panel. In both the pessimistic and optimistic scenarios, we can observe that WL has a larger constraining power than photometric void clustering, except for h and Ω_b in the optimistic case, where for h void clustering provides better constraints than WL by a factor of ∼3, while for Ω_b the constraints from the two single probes look comparable. This is due to the form of the WL kernel and the integration along the line of sight, and to the presence of baryon acoustic oscillations (BAO) in the void angular auto power spectrum, C^vv(ℓ), and the void-lensing cross spectrum, C^vγ(ℓ), (see Fig. 1), which are instead completely washed out in the WL angular power spectrum, C^γγ(ℓ). We also remark that, in general, both in the pessimistic and optimistic scenarios, the parameters for which the constraints improve most, when combining void clustering and WL, are not only h and Ω_b, but also n_s; similar conclusions (without WL) were found by Kreisch et al. (2021).

In the case of n_s, this happens since the WL and void clustering ellipses in the n_s − Ω_b plane happen to be orthogonal. We verified that, even if n_s and Ω_b produce an increment in the same direction for C^γγ(ℓ) and C^vv(ℓ), which would imply that the two parameters are negatively correlated, this orthogonality comes from projecting the large parameter space onto the n_s − Ω_b 2D space; in other words, it is due to the Fisher matrix inversion.

Finally, when moving from the pessimistic to the optimistic scenario, all marginalized errors decrease: in particular the error on n_s decreases by a factor of ∼2, the constraints on σ₈ and Ω_m get ∼40% tighter, and the uncertainty on h decreases by a factor of ∼3. The increase in the constraining power on σ₈ is somehow expected: the void bias affects the overall amplitude of the C(ℓ), so it is degenerate with σ₈. On the one hand, in the pessimistic scenario, both the value of the bias parameters and σ₈ were assumed to be measured from data, and this increases the uncertainties on σ₈. On the other hand, in the optimistic scenario, we assumed the cosmological dependence of the effective void bias to be known and we exploited it in the forecasts of the cosmological parameters.

5.2. The massive neutrino cosmology

Here we present the results in the case where the total neutrino mass, M_ν, is not kept fixed to its fiducial value, but it is a quantity whose value has to be determined from the data. Strictly speaking, in the case of the neutrino mass, the Fisher matrix approach goes beyond its range of validity, as the associated likelihood is non-Gaussian being truncated at M_ν = 0. However, as shown in Fig. 1 of Brinckmann & Lesgourgues (2019), the posterior obtained from the MCMC is in good agreement with the one obtained from the truncated Fisher matrix. Furthermore, we notice that, in the case of WL, the neutrino constraints presented in this work do not represent the final Euclid-WL constraints, as intrinsic alignments are not taken into account in this analysis. Moreover, a varying M_ν is not considered when the spectroscopic galaxy clustering, GC_sp, is included in Sect. 5.4. The contour plots for this cosmology are reported in Fig. 5, with the pessimistic case on the left and the optimistic one on the right. On the one hand, adding M_ν as a free parameter mainly impacts the measurement of σ₈, weakening its constraints by a factor of ∼2 with respect to the baseline ΛCDM case. This is expected since the effects of σ₈ and massive neutrinos on the matter power spectrum are similar: the former regulates the normalization of P_mm, and the latter suppresses P_mm in a scale-dependent way, due to neutrinos free streaming (Lesgourgues & Pastor 2006). On the other hand, such a suppression allows void clustering to impact positively on the neutrino mass measurements: the error on M_ν decreases with respect to the WL case alone, by ∼5% and ∼15% in the pessimistic and optimistic scenarios, respectively. In addition to considering a fiducial total mass of M_ν = 0.06 eV in the normal hierarchy scenario, we also computed parameter forecasts assuming a neutrino mass degenerate scenario with fiducial total mass M_ν = 0.15 eV⁸, from the best-fit value in Pellejero-Ibanez et al. (2017). With this choice, the marginalized errors on M_ν decrease by ∼5%. This is expected: the higher the neutrino mass, the more the matter power spectrum is sensitive to its effects (Lesgourgues & Pastor 2006).

5.3. The dynamical dark energy cosmology

In this cosmological scenario, we analyzed the w₀w_aCDM model, considering a time-dependent DE equation of state. The contour plots for this cosmology are reported in Fig. 6, both in the pessimistic and optimistic scenarios. Again, adding w₀ and w_a mainly affects h, σ₈, and Ω_m, increasing their errors. The impact on h and Ω_m is explained considering that w₀ and w_a enter the Hubble parameter H(z) (22), where the dependence from the DE equation of state has been exploited. Moreover, w₀ and w_a also enter the linear growth factor (Linder & Jenkins 2003), and therefore impact Ω_m and the normalization σ₈ again. When moving from the pessimistic to the optimistic scenario, the constraints on w₀ and w_a improve, and consequently the FoM is enhanced by a factor of ∼3. This can be explained since, in the optimistic scenario, we exploited both the cosmology dependence of the growth factor in the void bias evolution, b_v, as well as the cosmology dependence of the void size function used to compute the effective void bias in Eq. (19). This is confirmed by previous works (Pisani et al. 2015; Verza et al. 2019), although these focused on upcoming spectroscopic data and found that the void size function is sensitive to the DE equation of state. We also stress the increase in the constraining power when the two probes are combined together: comparing from Table 5 the constraints produced by WL alone against the ones obtained from its combination with the void clustering and void-lensing cross-correlation, we can observe that the FoM is enhanced by ∼10% in the pessimistic scenario, and by a factor of ∼2 in the optimistic one.

Finally, in order to measure the impact on parameter forecasts of the void-lensing cross-correlation signal, we also evaluated the marginalized errors when WL and void clustering are treated as independent probes, that is, when the Fisher matrices of the single probes are directly summed up. Then, we compared the FoM against the case when the void-lensing cross-correlation is included in the analysis. We found that, in the latter case, the DE FoM is enhanced by 5% in the pessimistic scenario, and by 20% in the optimistic one, as reported in Table 5.

5.4. Including the spectroscopic galaxy clustering

In this section we combine the spectroscopic galaxy clustering, GC_sp, to the WL, void clustering, and void-lensing cross-correlation probes. To this purpose, we considered GC_sp as a probe independent of the others, and made use of the GC_sp Fisher matrix as provided by EC20 for the observed anisotropic galaxy power spectrum, adding it to the ones computed in this work. Here our goal was to evaluate if and by how much the void clustering and the void-lensing correlation can still improve the Euclid performance even when its two primary probes are both accounted for. The corresponding results are reported in Fig. 7 and summarized in Table 6.

In the pessimistic scenario, constraints given by GC_sp + V are weaker than the ones given by GC_sp + WL; for instance, the FoM has a value of 15 in the former case, and of 112 in the latter one. We find that the constraints given by GC_sp + WL + V are only slightly improved with respect to GC_sp + WL. Nevertheless, the constraints given by GC_sp + WL + V + XC are tighter than in the GC_sp + WL case. In particular, the FoM increases by ∼5% when considering GC_sp + WL + V + XC against GC_sp + WL.

In the optimistic scenario, constraints given by GC_sp + V are comparable to the ones given by GC_sp + WL. Moreover, the FoM increases from 192 for GC_sp + WL, up to 702 for GC_sp + WL + V, in which case constraints on the other cosmological parameters become tighter: by ∼10% for h, by ∼30% for n_s, by ∼50 − 55% for Ω_b and σ₈, and by ∼70% for Ω_m.

Finally, when also including the void-lensing cross-correlation, while in the pessimistic case there is at most a 10% improvement in the FoM, in the optimistic setup we find that the most improved parameters are Ω_m and σ₈, with their constraints improved by a factor of 2 and 3, respectively; the total FoM increases by ∼10%, that is, from 702 for GC_sp + WL + V, up to 791 for GC_sp + WL + V + XC.

5.5. Combining massive neutrinos with dynamical dark energy

In this section we consider M_ν, w₀, and w_a as free parameters. The contour plots for this model cosmology are reported in Fig. 8, both for the pessimistic and optimistic cases. As expected, in this scenario the constraint on σ₈ gets even weaker with respect to the baseline ΛCDM case, due to the combined variation of M_ν, w₀, and w_a, entering the matter power spectrum, the Hubble parameter, and the growth factor. Comparing the pessimistic and optimistic scenarios, we observe improvements on the constraints in the latter case, but the orientation of the ellipses remains unchanged: this is reassuring as it means that the results do not change qualitatively when the two different setups are considered. In addition, we also performed a forecast for this model with a fiducial neutrino mass of M_ν = 0.33 eV with a degenerate mass spectrum, using the best-fit value in Pellejero-Ibanez et al. (2017): with this choice, the marginalized errors on M_ν decrease by ∼10% with respect to the fiducial case with M_ν = 0.06 eV. Again, this can be explained since, by increasing the value of M_ν, we increase the impact of neutrinos on the matter power spectrum, which suffers a larger suppression with respect to the fiducial neutrino mass case, becoming more sensitive to the presence of neutrinos.

From Table 5, it is possible to observe that, when we consider the case including the void-lensing cross-correlation signal compared to the one where the two probes are joined independently, the error on M_ν gets reduced by ∼5% and ∼10%, and the FoM is enhanced by 10% and 25%, in the pessimistic and optimistic scenarios, respectively.

As a final consideration, we comment on the mutual impact of M_ν and DE. Including w₀ and w_a increases the error on M_ν and vice versa. The increase in the w₀ and w_a errors, due to the variation of M_ν, is at most ∼15%, while the error on M_ν increases at most by 70% due to the variation of w₀ and w_a. However, especially for w₀ and w_a, the constraints are less impacted when the void-lensing cross-correlation is included in the analysis, compared to the constraints given by the single probes alone, as shown in the top panels of Fig. 9.

Fig. 9. Forecast contours for the νw0waCDM model, in the pessimistic (left) and in the optimistic (right) scenarios.

5.6. Systematics checks

The analysis presented in this work relies on various choices for the computation of the Fisher matrices, such as the number of bins in multipoles ℓ, redshift z, and wave-number k, as well as on differentiation methods. In this section we present the stability of the final results against different computation choices.

We performed stability tests where the multipole range was divided into 80, 100, and 120 logarithmically equispaced bins, where the C(ℓ) were evaluated in the linear center of the bin. The impact of this ℓ-binning was negligible: the marginalized 1-σ errors stayed unchanged up to the third digit. We also performed a forecast with a reduced multipole range, from ℓ = 20 to ℓ = 1350. Both for the pessimistic and optimistic void bias scenarios, the impact of the multipole range reduction on the FoM and M_ν 1-σ uncertainty was about 10%.

The impact of the k-binning was negligible too. CAMB evaluates the matter power spectrum P_mm(k, z) with a logarithmic k-binning grid. For the reported forecasts, P_mm(k, z) was evaluated with 60 k values per decade. Stability tests were performed also using 50 and 90 k values per decade; the marginalized 1-σ errors remained unchanged up to the third digit.

The impact of the z-binning was minor. CAMB evaluates the matter power spectrum P_mm(k, z) with a linear z-binning grid. For the reported forecast, P_mm(k, z) was evaluated for 450 z values. Stability tests were also performed using 300 and 600 z-bins: the differences in the marginalized 1-σ errors were always smaller than 3%.

The impact of the differentiation method was also negligible. In the reported forecast, the SteM derivative method was used. A forecast was also performed with the semi-analytical derivative method; the marginalized 1-σ errors stayed unchanged up to the fourth digit.

We also performed a forecast using equi-populated redshift bins for voids, rather than the equi-populated redshift bins for galaxies. This modification led to slightly (about 10 − 15%) worse constraints. This can be understood looking at the void distribution in Fig. 2, which increases with the redshift z. Thus the requirement of equi-populated redshift bins translates into broader (and so sparser) bins at low z, and narrower (and so denser) bins at high z, with respect to the equi-populated galaxy configuration. Since the effective void bias decreases in absolute value with z, this leads to more bins with a lower void bias, resulting in a decrease in the C(ℓ) in those bins, and hence to weaker constraints, with, as expected, a larger impact in the optimistic case than in the pessimistic case, where we marginalize over the void bias. Overall we conclude that our results are robust against various computation choices.

6. Conclusion and outlook

In this work we present the first forecast on cosmological parameter inference obtained combining the void-lensing cross-correlation with the WL and void angular two-point correlations, as will be measured from the Euclid photometric galaxy catalogs. In order to reach this goal we considered the following approach:

• we measured the projected void density distribution in redshift, n^v(z), directly from the Euclid Flagship photometric mock galaxy catalog, using the 2D void finder of Sánchez et al. (2017);

• we evaluated the C(ℓ) numerical derivatives, which enter the Fisher matrix expression, using the SteM technique (Camera et al. 2017). We tested the stability of results against different differentiation techniques, and found them to be robust;

• we used the void bias as obtained from the PBS and excursion set formalisms (Sheth & Van De Weygaert 2004), and considered the volume-conserving V dn void size function model (Jennings et al. 2013).

We present our results for two different choices of the effective void bias $b_{\text{v}}^{\text{eff}}(z)$: an “optimistic scenario”, where we exploited the cosmology dependence of both the void size function and void bias, and a “pessimistic scenario”, where we assumed the void bias evolution to be given by the growth factor in the reference cosmology, but its absolute normalization was supposed to be unknown, so that $b_{\text{v}}^{\text{eff}}(z=0)$ was marginalized over as a nuisance parameter.

We present parameter forecasts for different cosmological models: starting from a flat ΛCDM model, we first separately added as free parameters the total neutrino mass and the CPL parametrization of the DE equation of state (Chevallier & Polarski 2001; Linder 2002), then we add both in combination. Our main findings are presented in Table 5 and Figs. 4–8, and can be summarized as follows:

• ΛCDM cosmology: the WL angular correlation function is able to constrain all the cosmological parameters better than the photometric void angular spectrum; however, void constraints on h and Ω_b are competitive in the optimistic scenario. For the Hubble constant, h, this is due to the form of the weak lensing kernel and the integration along the line of sight. For Ω_b this is due to the presence of the BAO features in the void angular spectrum, which on the other hand are washed out in the lensing angular spectrum (see, e.g., Fig. 1). In general, these, together with n_s, are the parameters for which the constraints improve most when void clustering is combined with WL.

• νΛCDM cosmology: adding to the analysis the total neutrino mass, M_ν, as a free parameter mainly impacts the constraints on σ₈, as expected, since σ₈ and M_ν both affect the amplitude of the matter power spectrum (even if in the M_ν-case this is a scale-dependent effect). The void clustering, together with the void-lensing cross-correlation, improves the determination of the neutrino mass scale with respect to WL alone: the constraint on M_ν becomes stronger by ∼5% and ∼15%, respectively, in the pessimistic and optimistic scenarios.

• w₀w_aCDM cosmology: when the CPL parametrization of the DE equation of state is considered, we find that the most affected parameters are h, Ω_m, and σ₈. The expression of the Hubble parameter Eq. (22) explains the effect on h and Ω_m, while the impact of w₀ and w_a on the linear growth factor (Linder & Jenkins 2003) accounts for the effect on σ₈ and, again, Ω_m. Also, in this scenario void clustering helps to tighten the constraints on the cosmological parameters. In particular, the DE FoM increases by ∼10% in the pessimistic scenario and by a factor of ∼2 in the optimistic scenario.

• νw₀w_aCDM cosmology: in this model we evaluated, in particular, the mutual impact among w₀, w_a, and M_ν. On the one hand, we find that, when adding M_ν, the marginalized errors of w₀ and w_a increase at most by ∼15%. On the other hand, when adding w₀, w_a as free parameters, the impact on the M_ν constraints is higher, as its error increases even by 70%. However, especially for w₀ and w_a, the constraints are less impacted when the void-lensing cross-correlation is included, compared to the constraints provided by the two probes combined independently.

• FoM and the void-lensing cross-signal: in order to evaluate the improvement in the constraining power provided by adding the void-lensing cross-correlation signal to void clustering and weak lensing, we also considered the case where these two probes were assumed to be independent. When including the void-lensing correlation, the error on M_ν was reduced by 5% and 10%, in the pessimistic and optimistic scenarios, respectively. The FoM was enhanced by 10% and 25%, in the pessimistic and optimistic scenarios, respectively.

• Adding GC_sp: finally, we combined the galaxy lensing, void clustering, and void-lensing cross-correlation probes with the spectroscopic galaxy clustering, GC_sp. In this case, we exploited the Fisher matrices as obtained by EC20 in the so-called “pessimistic” and “optimistic” settings, and considered the flat w₀w_aCDM scenario, since GC_sp forecasts are computed keeping M_ν = 0.06 eV fixed. When combining WL + V + XC with GC_sp, both assumed in the pessimistic configuration, we find that the FoM increases from 23 for WL + V + XC alone up to 117 for WL + V + XC + GC_sp, and from 105 up to 791 in the optimistic case. In particular in the latter case, results are promising and competitive with other kinds of probe combination, as, for example, when galaxy lensing is combined with photometric galaxy clustering, GC_ph, and their cross-correlation (Tutusaus et al. 2020). As intrinsic alignments are not included in this analysis, we cannot make a proper comparison with the results reported in Table 14 in EC20, where the authors show a total FoM of 377 and 1257 in the pessimistic and optimistic settings, respectively. The results are, nevertheless, extremely encouraging: we show that including void clustering and the void-lensing cross-correlation in the total likelihood analysis will allow us to considerably improve Euclid’s performance.

The forecasts presented in this work show that photometric void clustering and its cross-correlation with WL deserve to be exploited in the data analysis of the Euclid galaxy survey, as they could be able to improve constraints on several cosmological parameters, such as, in particular, the total neutrino mass (Kreisch et al. 2021) and the DE equation of state. In this respect, it is worth noting that, for a full comparison with the Euclid performance from primary probes, in our analysis we should have included also the information from photometric galaxy clustering. However, we decided not to add this probe to the analysis in order to avoid double counting the information from the photometric sample. A way to include GC_ph, without incurring in this issue, is offered by a full-field inference approach, which allows us to self-consistently analyze GC and voids in the DM distribution (Leclercq et al. 2015). It would be important to also include GC_ph since, as shown in this work, the combination of void-clustering, void-lensing, and GC_sp does not allow us to reach the same Euclid performance as the combination of Euclid primary probes (EC20). However, this does not imply that the information coming from the photometric void sample should be neglected, as voids and galaxies are tracers with a very different bias, and hence probe different regimes of LSS (Wang & Zhao 2020), therefore providing therefore complementary information about our Universe. Moreover, in this paper we only consider only some observational and astrophysical uncertainties, neglecting baryonic physics and intrinsic alignment of galaxies.

While the aim of this work was to show that the inclusion of the void-void auto-correlation and void-lensing cross-correlation will improve the Euclid survey performance, and we performed some checks against the Flagship simulation to ensure the reliability of our theoretical modeling, additional work needs to be carried out to prepare the pipelines for the analyses of forthcoming real data. First and foremost, the theoretical model should be further improved. As already mentioned, we did not include nonlinear scales, since this would require a modeling of the void density profile and its dependence on cosmological parameters. Second, we need to rely on more realistic mock data, which will account for more effects: Euclid survey specification and systematics and possibly mispecifications induced by the HOD, which could induce some errors in the identified void catalogs. Finally, the framework presented here can also be easily extended to investigate non-flat scenarios, more specific DE models, theories of modified gravity, and primordial non-Gaussianities (Chan et al. 2019).

Acknowledgments

M.B. acknowledges financial support from the ASI agreement n. I/023/12/0 “Euclid attivitá relativa alla fase B2/C”. N.H. is supported by the Excellence Cluster ORIGINS, which is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC-2094 – 390783311. A.P. is supported by NASA ROSES grant 12-EUCLID12-0004, and NASA grant 15-WFIRST15-0008 to the Nancy Grace Roman Space Telescope Science Investigation Team “Cosmology with the High Latitude Survey”. A.K. has been supported by a Juan de la Cierva Incorporación fellowship with project number IJC2018-037730-I, and funding for this project was also available in part through SEV-2015-0548 and AYA2017-89891-P. S.N. acknowledges support from an STFC Ernest Rutherford Fellowship, grant reference ST/T005009/1. 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’Etudes 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 Economia y Competitividad, 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 website (http://www.euclid-ec.org).

References

Appendix A: Fisher matrix expression proof

This paper presents a forecast for cosmological parameter measurements using the Fisher matrix technique. This appendix is devoted to the derivation of the Fisher matrix expression considering as observables either the a_ℓm (field perspective) or the C(ℓ) (estimator perspective).

A.1. The field perspective

We consider the spherical harmonics’ expansion coefficients $a_{\ell m}^{\mathrm{A}i}$, the coefficients of the spherical harmonics decomposition of the 2D field A in the i-th tomographic bin, which we assume to follow a multivariate gaussian distribution with zero mean and covariance matrix S

$$\begin{array}{c}\hfill L({\mathit{a}}_{\ell }|\mathit{\theta })=\frac{exp(-\frac{1}{2}{\mathit{a}}_{\ell }^{\mathrm{T}}{\mathbf{S}}^{-1}{\mathit{a}}_{\ell })}{\sqrt{{(2\pi )}^n|\mathbf{S}|}},\end{array}$$(A.1)

where θ is the vector of the parameters and a_ℓ is a vector collecting the $a_{\ell m}^{\mathrm{A}i}$. The logarithm of this probability distribution function is

$$\begin{array}{cc}\hfill lnL({\mathit{a}}_{\ell }|\mathit{\theta })& =-\frac{1}{2}[nln(2\pi )+ln|\mathbf{S}|+{\mathit{a}}_{\ell }^{\mathrm{T}}{\mathbf{S}}^{-1}{\mathit{a}}_{\ell }]\hfill \\ \hfill & =-\frac{1}{2}[nln(2\pi )+ln|\mathbf{S}|+\text{Tr}({\mathbf{S}}^{-1}\mathbf{A})],\hfill \end{array}$$(A.2)

with $\mathbf{A}\equiv {\mathit{a}}_{\mathcal{l}}{\mathit{a}}_{\mathcal{l}}^T$. The computation of the second derivative, after taking the expectation value, yields

$$\begin{array}{c}\hfill F_{\alpha \beta }(\ell )=-\langle lnL{({\mathit{a}}_{\ell }|\mathit{\theta })}_{,\alpha \beta }\rangle =\frac{1}{2}\text{Tr}\left({\mathbf{S}}^{-1}{\mathbf{S}}_{,\alpha }{\mathbf{S}}^{-1}{\mathbf{S}}_{,\beta }\right),\end{array}$$(A.3)

where S = ⟨A⟩ and we used the following matrix identities:

$$\begin{array}{c}\hfill {(ln|\mathbf{S}|)}_{,\alpha }=\text{Tr}({\mathbf{S}}^{-1}{\mathbf{S}}_{,\alpha }),{({\mathbf{S}}^{-1})}_{,\alpha }=-{\mathbf{S}}^{-1}{\mathbf{S}}_{,\alpha }{\mathbf{S}}^{-1}.\end{array}$$(A.4)

Now, we specify how the $a_{\ell m}^{\mathrm{A}i}$ were collected into the vector a_ℓ. First of all, we recall that

$$\begin{array}{c}\hfill \langle a_{\ell m}^{\mathrm{A}i}{a}_{{\ell }^{\prime }{m}^{\prime }}^{\mathrm{B}j}\rangle ={\mathrm{\Sigma }}_{\mathit{ij}}^{\mathrm{AB}}(\ell ){\delta }_{\ell {\ell }^{\prime }}{\delta }_{m{m}^{\prime }},{\mathrm{\Sigma }}_{\mathit{ij}}^{\mathrm{AB}}(\ell )=C_{\mathit{ij}}^{\mathrm{AB}}(\ell )+N_{\mathit{ij}}^{\mathrm{AB}}(\ell ),\end{array}$$(A.5)

where $N_{\mathit{ij}}^{\mathrm{AB}}$ is the Poisson shot noise for probes combination AB in tomographic bins ij. Working at fixed multipole ℓ, the vector a_ℓ has a multi-index I = (m, i, A). Each of these indices runs in a different range, and the range of the tomographic index i may depend on the probe A considered. To summarize:

• m varies from −ℓ to ℓ

• i varies from 1 to N_A

• A varies from 1 to 𝒩,

with 𝒩 being the number of probes and N_A the number of tomographic bins for probe A. We chose to order the array a_ℓ varying the three indices with a significance increasing from left to right. With this choice the matrix S has the following block form:

$$\begin{array}{c}\hfill \mathbf{S}=\mathrm{\Sigma }(\ell )\otimes {\mathbb{1}}_{2\ell +1},\mathrm{\Sigma }(\ell )=\left(\begin{array}{cccc}{\mathrm{\Sigma }}^{11}(\ell )& {\mathrm{\Sigma }}^{12}(\ell )& \cdots & {\mathrm{\Sigma }}^{1\mathcal{N}}(\ell )\\ {\mathrm{\Sigma }}^{21}(\ell )& {\mathrm{\Sigma }}^{22}(\ell )& \cdots & {\mathrm{\Sigma }}^{2\mathcal{N}}(\ell )\\ ⋮& ⋮& \ddots & ⋮\\ {\mathrm{\Sigma }}^{\mathcal{N}1}(\ell )& {\mathrm{\Sigma }}^{\mathcal{N}2}(\ell )& \cdots & {\mathrm{\Sigma }}^{\mathcal{N}\mathcal{N}}(\ell )\end{array}\right).\end{array}$$(A.6)

Here Σ^AB(ℓ) is the tomographic covariance matrix, that is to say, ${[{\mathrm{\Sigma }}^{\mathrm{AB}}(\ell )]}_{\mathit{ij}}={\mathrm{\Sigma }}_{\mathit{ij}}^{\mathrm{AB}}(\ell )$, and ⊗ is the Kronecker product. Equation (A.6) can be understood as follows: the matrix S is made up of diagonal blocks of size ℓ, each of them proportional to the identity 𝟙_2ℓ+1 with a different factor ${\mathrm{\Sigma }}_{\mathit{ij}}^{\mathrm{AB}}(\ell )$. Here the identity matrix 𝟙_2ℓ+1 is exactly the Kronecker delta δ_mm′, which appears in Eq. (A.5). In order to compute the trace in Eq. (A.3) we make use of some properties of the Kronecker product (Steeb 1997):

$$\begin{array}{cc}& (\mathbf{A}\otimes \mathbf{B})(\mathbf{C}\otimes \mathbf{D})=\mathbf{A}\mathbf{C}\otimes \mathbf{B}\mathbf{D},\hfill \end{array}$$(A.7)

$$\begin{array}{cc}\hfill & {(\mathbf{A}\otimes \mathbf{B})}^{-1}={\mathbf{A}}^{-1}\otimes {\mathbf{B}}^{-1},\hfill \end{array}$$(A.8)

$$\begin{array}{cc}\hfill & \text{Tr}(\mathbf{A}\otimes \mathbf{B})=\text{Tr}(\mathbf{A})\text{Tr}(\mathbf{B}).\hfill \end{array}$$(A.9)

Using these properties, we have

$$\begin{array}{cc}\hfill \text{Tr}\left[{\mathbf{S}}^{-1}{\mathbf{S}}_{,\alpha }{\mathbf{S}}^{-1}{\mathbf{S}}_{,\beta }\right]& =\text{Tr}[\left(\mathrm{\Sigma }{(\ell )}^{-1}\mathrm{\Sigma }{(\ell )}_{,\alpha }\mathrm{\Sigma }{(\ell )}^{-1}\mathrm{\Sigma }{(\ell )}_{,\beta }\right)\otimes {\mathbb{1}}_{2\ell +1}]\hfill \\ \hfill & =\text{Tr}\left(\mathrm{\Sigma }{(\ell )}^{-1}\mathrm{\Sigma }{(\ell )}_{,\alpha }\mathrm{\Sigma }{(\ell )}^{-1}\mathrm{\Sigma }{(\ell )}_{,\beta }\right)\text{Tr}\left({\mathbb{1}}_{2\ell +1}\right)\hfill \\ \hfill & =(2\ell +1)\text{Tr}\left(\mathrm{\Sigma }{(\ell )}^{-1}\mathrm{\Sigma }{(\ell )}_{,\alpha }\mathrm{\Sigma }{(\ell )}^{-1}\mathrm{\Sigma }{(\ell )}_{,\beta }\right)\hfill \\ \hfill & =(2\ell +1)\text{Tr}\left(\mathrm{\Sigma }{(\ell )}^{-1}\mathbf{C}{(\ell )}_{,\alpha }\mathrm{\Sigma }{(\ell )}^{-1}\mathbf{C}{(\ell )}_{,\beta }\right),\hfill \end{array}$$

where in the last line we used the fact that the shot noise is independent of the cosmological parameters, and the C(ℓ)’s block matrix C(ℓ) has the same form of Σ(ℓ) of Eq. (A.6), with entries $C_{ij}^{\text{AB}}(\mathcal{l})$ instead of $\mathrm{\Sigma }{(\ell )}_{\mathit{ij}}^{\mathrm{AB}}$. Then from Eq. (A.3), we can finally write

$$\begin{array}{c}\hfill F_{\alpha \beta }(\ell )=\frac{2\ell +1}{2}\text{Tr}\left(\mathrm{\Sigma }{(\ell )}^{-1}\mathbf{C}{(\ell )}_{,\alpha }\mathrm{\Sigma }{(\ell )}^{-1}\mathbf{C}{(\ell )}_{,\beta }\right).\end{array}$$

Now, we sum over the multipoles ℓ, considering them as independent

$$\begin{array}{c}\hfill F_{\alpha \beta }={\displaystyle \sum _{\ell }\frac{2\ell +1}{2}\text{Tr}\left[\mathrm{\Sigma }{(\ell )}^{-1}\mathbf{C}{(\ell )}_{,\alpha }\mathrm{\Sigma }{(\ell )}^{-1}\mathbf{C}{(\ell )}_{,\beta }\right].}\end{array}$$(A.10)

In order to obtain the expression Eq. (43), it sufficient to redefine the covariance as

$$\begin{array}{c}\hfill \mathrm{\Sigma }(\ell )⟶\sqrt{\frac{2}{(2\ell +1)\mathrm{\Delta }\ell f_{\mathrm{sky}}}}\mathrm{\Sigma }(\ell ).\end{array}$$(A.11)