© J. Kuruvilla 2022

1. Introduction

One of the outstanding questions in modern physics is related to the determination of the mass of neutrinos, which has fundamental implications for both particle physics and cosmology. Neutrino oscillation experiments have established that neutrinos ought to have a non-zero mass (e.g., Forero et al. 2014; Gonzalez-Garcia et al. 2016; Capozzi et al. 2017; de Salas et al. 2017). However, these oscillation experiments are only sensitive to the mass splittings between the neutrino mass eigenstates and in order to measure the absolute scale of the neutrino mass, other experiments are required. Recently, the Karlsruhe Tritium Neutrino (KATRIN) experiment reported the first direct detection of sub-eV neutrino mass, with an upper limit on the ‘effective neutrino mass’ of 0.8 eV (Aker et al. 2021). This is based on the kinematic measurements through the observation of the energy spectrum of tritium β-decay in a way that is model independent.

All the while, stronger constraints on the summed neutrino mass (M_ν) can be obtained by combining various cosmological probes since the massive neutrinos leave an imprint on various cosmological observables (e.g., Wong 2011; Lesgourgues & Pastor 2012). Yet these constraints have an additional model dependence. Assuming the standard ‘lambda cold dark matter’ (ΛCDM) model, one of the strongest constraints on the summed neutrino mass has been obtained by combining the cosmic microwave background (CMB, Planck Collaboration VI 2020), baryonic acoustic oscillation, and redshift-space galaxy clustering (eBOSS Collaboration 2021) to obtain an upper limit of M_ν < 0.102 eV. However, by considering extensions of the standard cosmological model, the upper limit becomes less stringent (e.g., Vagnozzi et al. 2018; Choudhury & Hannestad 2020). The community can expect the summed neutrino mass to be measured with increased precision from cosmological probes in the foreseeable future with the advent of next generation of CMB surveys (e.g., the Simons Observatory¹ (SO, Ade et al. 2019), CMB-S4² (Abazajian et al. 2016)) as well as the stage IV galaxy redshift surveys (e.g., the ‘Dark Energy Spectroscopic Instrument’³ (DESI, DESI Collaboration 2016), Euclid⁴ (Laureijs et al. 2011) and Nancy Grace Roman space telescope⁵ (Spergel et al. 2015)).

Currently, the M_ν constraints from galaxy clustering are mainly obtained using two-point statistics, namely, the power spectrum in the Fourier space or the two-point correlation function in the configuration space. The impact of massive neutrinos on the two-point clustering statistics has been studied quite extensively using N-body simulations both in the real (e.g., Saito et al. 2008; Wong 2008; Castorina et al. 2015) and the redshift space (e.g., Villaescusa-Navarro et al. 2018; García-Farieta et al. 2019). However, these statistics are affected by the M_ν − σ₈ degeneracy, thus acting as a limitation in measuring the summed neutrino mass. The three-point clustering statistics in the Fourier space (i.e., the bispectrum) has been shown to break this degeneracy (Hahn et al. 2020; Hahn & Villaescusa-Navarro 2021). In addition, it has been shown that the three-point cluster statistics contain additional cosmological information compared to its two-point counterpart and it is thus able to obtain substantial improvements with regard to constraints on other cosmological parameters as well (e.g., Yankelevich & Porciani 2019; Chudaykin & Ivanov 2019; Gualdi & Verde 2020; Agarwal et al. 2021; Samushia et al. 2021). Currently there are efforts to understand the possibility of constraining summed neutrino mass using various summary statistics, among others, such as the one-point probability distribution function of the matter density (e.g., Uhlemann et al. 2020) and the void size function (e.g., Bayer et al. 2021).

Another plausible avenue is to use velocity statistics such as the mean pairwise velocity, which provides a complementary view to the clustering information, either through the peculiar velocity surveys or the kinetic Sunyaev–Zeldovich (kSZ) effect. Mueller et al. (2015) showed that the mean pairwise velocity can be utilised to constrain the summed neutrino mass, while Kuruvilla et al. (2020) studied its interplay between the baryonic feedback and the summed neutrino mass effects at nonlinear separation scales. Furthermore, the three-point mean relative velocity statistics is able to obtain stronger constraints on the summed neutrino mass when compared to the mean pairwise velocity (Kuruvilla & Aghanim 2021). However, the growth rate measurement from the kSZ effect of the CMB are degenerate with the optical depth (e.g., Keisler & Schmidt 2013; Battaglia 2016; Flender et al. 2017) and this degeneracy acts as a limitation in measuring the cosmological parameters (e.g., Smith et al. 2018), which is commonly referred to as the optical depth degeneracy. It has been suggested that the use of fast radio bursts can be used to break this degeneracy (Madhavacheril et al. 2019). In this paper, we develop a new statistic that is independent of the optical depth and based on the use of the first moment of the three-point relative velocities, namely, the mean relative velocities between pairs in a triplet. Thus, we are able to circumvent the problem of optical depth as a limitation factor in kSZ experiments – thereby formulating one of the main goals of this paper.

The remainder of this work is structured as follows. In Sect. 2, we describe the newly introduced summary statistic based on three-point mean relative velocities. The Quijote suite of simulation used in this work is introduced briefly in Sect. 3.1. The information content in the velocity statistics is studied using the Fisher-matrix formalism, which is briefly summarised in Sect. 3.2. Our results are presented in Sect. 4. We present our conclusions in Sect. 5.

2. Cosmology using the kSZ effect

2.1. Kinetic Sunyaev–Zeldovich effect

As the CMB photons interact with the free electrons of hot ionised gas along the line of sight (LOS), the apparent CMB temperature changes. This is due to the fact that there is a transfer of energy from electrons to the resulting scattered photons as the electrons have a significantly higher kinetic energy than the photons. In this work, we focus on the secondary effect that is known as the kinetic Sunyaev–Zeldovich (kSZ; Sunyaev & Zeldovich 1972, 1980). It arises when the scattering medium is moving relative to the Hubble flow. The fractional temperature fluctuation caused due to kSZ is

$$\begin{array}{cc}\hfill {\frac{\mathrm{\Delta }T(\widehat{n})}{T_{\mathrm{cmb}}}|}_{\mathrm{kSZ}}& =-{\displaystyle \int }\mathrm{d}l{\sigma }_{\mathrm{T}}\left(\frac{{\mathit{v}}_{\mathrm{e}}·\widehat{n}}{c}\right)n_{\mathrm{e}},\hfill \\ \hfill & =-\tau \left(\frac{{\mathit{v}}_{\mathrm{e}}·\widehat{n}}{c}\right),\hfill \end{array}$$(1)

where σ_T is the Thomson scattering cross-section, T_cmb is the CMB temperature, c is the speed of light, v_e is the peculiar velocity of free electrons, and n_e is the physical free electron number density. The integral ∫dl is computed along the LOS which is given by $\widehat{n}$. The optical depth is defined as τ = ∫dlσ_Tn_e, that is, the integrated electron density.

The kSZ signal detection is challenging because of its small amplitude and the fact that its spectrum is identical to that of primary CMB temperature fluctuations. One of the approaches used to detect the kSZ signal is to employ the pairwise statistic (e.g., Hand et al. 2012; Hernández-Monteagudo et al. 2015; Planck Collaboration Int. XXXVII 2016; Schaan et al. 2016; Soergel et al. 2016; De Bernardis et al. 2017; Li et al. 2018; Calafut et al. 2021; Chen et al. 2022). There is existing evidence of the kSZ signal using other techniques as well (e.g., Hernández-Monteagudo et al. 2015; Schaan et al. 2016, 2021; Nguyen et al. 2020; Tanimura et al. 2021; Chaves-Montero et al. 2021).

In the case of the kSZ pairwise signal, the temperature acts as a proxy for the peculiar velocity, and, as such, it probes the optical depth weighted pairwise velocity (e.g., Hand et al. 2012; Soergel et al. 2018)

$$\begin{array}{c}\hfill \frac{⟨\mathrm{\Delta }{T}^{\mathrm{kSZ}}(r_{12})⟩}{T_{\mathrm{cmb}}}\simeq -\tau \frac{\overline{w}(r_{12})}{c},\end{array}$$(2)

where $⟨\mathrm{\Delta }{T}_{12}^{\mathrm{kSZ}}⟩$ is the mean temperature difference between the objects ‘1’ and ‘2’, and $\overline{w}(r_{12})$ is the mean radial component of the pairwise velocity that can be defined in the single streaming regime as follows:

$$\begin{array}{cc}\hfill {⟨{\mathit{w}}_{12}|{\mathit{r}}_{12}⟩}_{\mathrm{p}}=& {\displaystyle \frac{⟨(1+{\delta }_1)(1+{\delta }_2)({\mathit{v}}_2-{\mathit{v}}_1)⟩}{⟨(1+{\delta }_1)(1+{\delta }_2)⟩},}\hfill \end{array}$$(3)

where δ_i ≡ δ(x_i) is the density contrast, v_i ≡ v(x_i)≡u(x_i)/aH is the normalised peculiar velocity, a is the scale factor, and H is the Hubble constant. Using perturbation theory at leading order (LO), the mean radial matter pairwise velocity can be written as (e.g., Fisher 1995; Juszkiewicz et al. 1998; Reid & White 2011)

$$\begin{array}{c}\hfill {\displaystyle {⟨{\mathit{w}}_{12}|{\mathit{r}}_{12}⟩}_{\mathrm{p}}=\overline{w}(r_{12}){\widehat{\mathit{r}}}_{12}\simeq -\frac{f}{{\pi }^2}{\widehat{\mathit{r}}}_{12}{{\displaystyle \int }}_0^{\infty }k{j}_1(k{r}_{12})P(k)\mathrm{d}k,}\end{array}$$(4)

where ${\widehat{\mathit{r}}}_{12}$ is the unit vector along the pair ‘12’, the subscript p implies that the averages are computed over all pairs with separation r₁₂, P(k) denotes the linear matter power spectrum, and j₁(x) = sin(x)/x² − cos(x)/x. It should be noted that Eq. (2) assumes that there is no correlation between optical depth and velocity field. Following Eqs. (2) and (4), we can see that

$$\begin{array}{c}\hfill \mathrm{\Delta }{T}^{\mathrm{kSZ}}\propto \tau f{\sigma }_8^2,\end{array}$$(5)

thus implying that the growth rate measurement from the pairwise kSZ is perfectly degenerate with optical depth (Keisler & Schmidt 2013). Here, we presented the argument for the matter component, however, in the observations there is an additional bias dependence entering in the above equation, which comes about as a result of the galaxy bias and while assuming no velocity bias.

2.2. New statistics based on mean relative velocity between pairs in a triplet

In the previous section, we mentioned about the mean relative velocity between two tracers (i.e., between a pair) or the mean pairwise velocity. However, this can be generalised to the case of three tracers, in which we can consider two mean relative velocity between pairs in a triplet with separations △₁₂₃ = (r₁₂, r₂₃, r₃₁): (i) ⟨w₁₂|△₁₂₃⟩_t, and (ii) ⟨w₂₃|△₁₂₃⟩_t. The subscript t here implies that the averages are computed over all triplets with separations (r₁₂, r₂₃, r₃₁). Similarly to Eq. (3), in the single stream fluid approximation, the mean relative velocity between pair 12 in a triplet can be written as (Kuruvilla & Porciani 2020):

$$\begin{array}{cc}\hfill {⟨{\mathit{w}}_{12}|{▵}_{123}⟩}_{\mathrm{t}}& {\displaystyle =\frac{⟨(1+{\delta }_1)(1+{\delta }_2)(1+{\delta }_3)({\mathit{v}}_2-{\mathit{v}}_1)⟩}{⟨(1+{\delta }_1)(1+{\delta }_2)(1+{\delta }_3)⟩}}\hfill \\ \hfill & \simeq ⟨{\delta }_1{\mathit{v}}_2⟩-⟨{\delta }_2{\mathit{v}}_1⟩+⟨{\delta }_3{\mathit{v}}_2⟩-⟨{\delta }_3{\mathit{v}}_1⟩\hfill \\ \hfill & =\overline{w}(r_{12}){\widehat{\mathit{r}}}_{12}-\frac{1}{2}[\overline{w}(r_{23}){\widehat{\mathit{r}}}_{23}+\overline{w}(r_{31}){\widehat{\mathit{r}}}_{31}].\hfill \end{array}$$(6)

The three-point mean relative velocity statistics can be composed into both its radial (R_ij) and transverse (T_ij) component in the plane of the triangle defined by the particles. This is in contrast to the mean pairwise velocity for which the transverse component is zero. In the case of ⟨w₁₂|△₁₂₃⟩_t, it follows as:

$$\begin{array}{cc}\hfill {⟨{\mathit{w}}_{12}|{▵}_{123}⟩}_{\mathrm{t}}& ={⟨{\mathit{w}}_{12}·{\widehat{\mathit{r}}}_{12}|{▵}_{123}⟩}_{\mathrm{t}}{\widehat{\mathit{r}}}_{12}+{⟨{\mathit{w}}_{12}·\widehat{\mathit{t}}|{▵}_{123}⟩}_{\mathrm{t}}\widehat{\mathit{t}}\hfill \\ \hfill & =R_{12}({▵}_{123}){\widehat{\mathit{r}}}_{12}+T_{12}({▵}_{123})\widehat{\mathit{t}},\hfill \end{array}$$(7)

where $\widehat{\mathit{t}}=({\widehat{\mathit{r}}}_{23}-cos\chi {\widehat{\mathit{r}}}_{12})/sin\chi$, ${\widehat{\mathit{r}}}_{23}$ is the unit vector along the pair ‘23’, and $\chi =arccos({\widehat{\mathit{r}}}_{12}·{\widehat{\mathit{r}}}_{23})$. In this work, we make use of only the radial component and for the pair marked 12 for the triplet, it can be written as

$$\begin{array}{cc}\hfill R_{12}({▵}_{123})& =\overline{w}(r_{12})-\frac{1}{2}[\overline{w}(r_{23})cos\chi \hfill \\ \hfill & -\overline{w}(r_{31})\frac{r_{12}+r_{23}cos\chi }{\sqrt{r_{12}^2+r_{23}^2+2r_{12}{r}_{23}cos\chi }}].\hfill \end{array}$$(8)

Similarly, the mean radial relative velocity between the pair 23 in △₁₂₃ can be written as

$$\begin{array}{cc}\hfill R_{23}({▵}_{123})& =\overline{w}(r_{23})-\frac{1}{2}[\overline{w}(r_{12})cos\chi \hfill \\ \hfill & -\overline{w}(r_{31})\frac{r_{23}+r_{12}cos\chi }{\sqrt{r_{12}^2+r_{23}^2+2r_{12}{r}_{23}cos\chi }}].\hfill \end{array}$$(9)

Similarly to Eq. (2), the three-point mean relative temperature difference from kSZ can be expressed as:

$$\begin{array}{c}\hfill \frac{\mathrm{\Delta }{T}_{\mathit{ij}}^{\mathrm{kSZ}}({▵}_{123})}{T_{\mathrm{cmb}}}\simeq -\tau \frac{R_{\mathit{ij}}^{\mathrm{h}}({▵}_{123})}{c},\end{array}$$(10)

where $R_{\mathit{ij}}^{\mathrm{h}}({▵}_{123})$ respresents the three-point mean relative velocity statistics for haloes (biased tracers), and to first approximation it can be written down as linear bias term times R_ij(△₁₂₃) (Kuruvilla & Aghanim 2021). Based on the radial mean relative velocities between pairs in a triplet, we can introduce a new ratio statistic, namely:

$$\begin{array}{cc}\hfill \mathcal{R}({▵}_{123})& =\frac{{⟨{\mathit{w}}_{12}·{\widehat{\mathit{r}}}_{12}|{▵}_{123}⟩}_{\mathrm{t}}}{{⟨{\mathit{w}}_{23}·{\widehat{\mathit{r}}}_{23}|{▵}_{123}⟩}_{\mathrm{t}}},\hfill \end{array}$$(11)

which tells us how quickly the average infall velocity of pair 12 performs in comparison to the average infall velocity of pair 23 for a specific triangular configuration △₁₂₃. On linear scales, using perturbation theory at LO, ℛ(△₁₂₃) can be written as the ratio between Eqs. (8) and (9)

$$\begin{array}{c}\hfill \mathcal{R}({▵}_{123})=\frac{R_{12}({▵}_{123})}{R_{23}({▵}_{123})}=\frac{R_{12}^{\mathrm{h}}({▵}_{123})}{R_{23}^{\mathrm{h}}({▵}_{123})}\equiv \frac{\mathrm{\Delta }{T}_{12}^{\mathrm{kSZ}}({▵}_{123})}{\mathrm{\Delta }{T}_{23}^{\mathrm{kSZ}}({▵}_{123})}.\end{array}$$(12)

The above statistic is thus independent of optical depth, σ₈, and linear bias. In the following sections, we take a detailed look at whether the Ansatz of σ₈ and linear bias independence holds up. Additionally, we study the cosmological information content in ℛ(△₁₂₃).

3. Data and analysis

3.1. Quijote simulation suite

In this work, we make use of the Quijote⁶ (Villaescusa-Navarro et al. 2020) suite of simulations, which was run using the tree-PM code GADGET-3 (Springel 2005). Spanning more than a few thousand cosmological models, it contains 44 100 N-body simulations. These simulations have a box length of 1 h⁻¹ Gpc, tracking the evolution of 512³ cold dark matter (CDM) particles. The initial conditions (ICs) were generated at redshift z = 127 using the second-order Lagrangian perturbation theory. The fiducial cosmological parameters (assuming zero-summed neutrino mass) for the simulation is as follows: the total matter density: Ω_m = 0.3175; baryonic matter density: Ω_b = 0.049; primordial spectral index of the density perturbations: n_s = 0.9624; amplitude of the linear power spectrum on the scale of 8 h⁻¹ Mpc: σ₈ = 0.834; and the present-day value of the Hubble constant: H₀ ≡ H(z = 0) = 100 h km s⁻¹ Mpc⁻¹, with h = 0.6711. This is broadly consistent with the Planck 2018 result (Planck Collaboration VI 2020). The suite consists of 15 000 random realisations for the fiducial cosmology. For the purpose of calculating derivatives, Quijote provides a set of 500 random realisations wherein only one parameter is varied with respect to the fiducial cosmology. The variations are as follows: $\{{\mathrm{\Omega }}_{\text{m}}^{+},{\mathrm{\Omega }}_{\text{m}}^{-},{\mathrm{\Omega }}_{\text{b}}^{+},{\mathrm{\Omega }}_{\text{b}}^{-},n_{\text{s}}^{+},n_{\text{s}}^{-},{\sigma }_8^{+},{\sigma }_8^{-}\}=\{0.3275,0.3075,0.051,0.047,0.9824,0.9424,0.849,0.819\}$ and {h⁺, h⁻}={0.6911, 0.6511}.

In addition, the suite also provides 500 realisations for three massive neutrino cosmology, where the summed neutrino masses are 0.1, 0.2, and 0.4 eV, respectively. The initial conditions for these simulations were produced using the Zeldovich approximation (ZA), with 512³ neutrino particles in addition to the CDM particles. To compute the numerical derivatives with respect to massive neutrinos, the Quijote suite provides an addition 500 random realisations for the fiducial cosmology, in which the ICs were also generated using ZA.

In this work, we use the halo catalog data from 22 000 N-body simulations of the Quijote suite. These halos were identified using a friends-of-friends algorithm. We selected halos that have a halo mass of M_h > 5 × 10¹³h⁻¹ M_⊙ (corresponding to groups and clusters of galaxies) at z = 0, which gives a mean number density of $\overline{n}\sim 0.92\times {10}^{-4}{h}^3{\mathrm{Mpc}}^{-3}$ for the reference simulations. Additionally, in the case of (i) fiducial cosmology, and (ii) for variations in σ₈ (both ${\sigma }_8^{+}$ and ${\sigma }_8^{-}$), we used 30 realisations of the particle data (randomly down-sampled to 100³ particles) to compute ℛ(△₁₂₃).

3.2. Fisher-matrix formalism

To quantify the error estimates on the cosmological parameters, we used the Fisher-matrix formalism which can be defined as (e.g., Tegmark et al. 1997; Heavens 2009; Verde 2010):

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

where θ_α and θ_β are two of the cosmological model parameters and ℒ is the likelihood of the data given a model. Assuming a Gaussian likelihood, we can express the Fisher information matrix as

$$\begin{array}{c}\hfill F_{\alpha \beta }=\frac{\partial \mathcal{R}}{\partial {\theta }_{\alpha }}·{\widehat{\mathbf{C}}}^{-1}·\frac{\partial {\mathcal{R}}^{\mathsf{T}}}{\partial {\theta }_{\beta }},\end{array}$$(14)

where ℛ represents the data vector for the ratio statistic we introduced in Eq. (11) and ${\widehat{\mathbf{C}}}^{-1}$ is the precision matrix (i.e., the inverse covariance matrix). It should be noted that in the definition of F_αβ, we neglected a term that appears due to the cosmology dependence of the covariance matrix. However, the correction has been shown to have a negligible effect (Kodwani et al. 2019). We compute the covariance matrix of ℛ directly from the simulations as follows:

$$\begin{array}{c}\hfill \stackrel{\sim }{\mathbf{C}}=\frac{1}{N_{\mathrm{sims}}-1}{\displaystyle \sum _{i=1}^{N_{\mathrm{sims}}}({\mathcal{R}}_i-\overline{\mathcal{R}}){({\mathcal{R}}_i-\overline{\mathcal{R}})}^{\mathsf{T}},}\end{array}$$(15)

where $\overline{\mathcal{R}}=N_{\mathrm{sims}}^{-1}\sum _{i=1}^{N_{\mathrm{sims}}}{\mathcal{R}}_i$, and N_sims denotes the total number of simulations used to compute the covariance matrix (in this work N_sims = 15 000). While Eq. (15) gives an unbiased estimate of the covariance matrix, its inversion leads to a biased estimate of the precision matrix. This, however, can be statistically corrected by applying a multiplicative correction factor to the precision matrix (Kaufmann 1967; Anderson 2003; Hartlap et al. 2007):

$$\begin{array}{c}\hfill {\widehat{\mathbf{C}}}^{-1}=\frac{N_{\mathrm{sims}}-N_{\mathrm{bins}}-2}{N_{\mathrm{sims}}-1}{\stackrel{\sim }{\mathbf{C}}}^{-1},\end{array}$$(16)

where N_bins is the number of bins in ℛ.

We numerically computed the derivatives required to construct the Fisher information matrix using the Quijote suite of simulations, which provides 500 realisations where only one cosmological parameter is varied while the rest are fixed at its fiducial value. Thus, in the case when the model parameters are one of the following: θ ≡ {Ω_m, Ω_b, h, n_s, σ₈}, we can make use of the central difference approximation to compute the derivative numerically, namely:

$$\begin{array}{c}\hfill \frac{\partial \mathcal{R}}{\partial \theta }\simeq \frac{\mathcal{R}(\theta +\mathrm{d}\theta )-\mathcal{R}(\theta -\mathrm{d}\theta )}{2\mathrm{d}\theta }·\end{array}$$(17)

In the case of the neutrino mass, the fiducial value is 0.0 eV and it cannot have negative values, hence, we obtain the partial derivative using

$$\begin{array}{c}\hfill \frac{\partial \mathcal{R}}{\partial M_{\nu }}\simeq \frac{-\mathcal{R}(M_{\nu }=0.4)+4\mathcal{R}(M_{\nu }=0.2)-\mathcal{R}(M_{\nu }=0)}{0.4}·\end{array}$$(18)

Thus, we utilise two sets of massive neutrino simulations from Quijote, with M_ν = 0.2 eV and M_ν = 0.4 eV for the Fisher information matrix. However, the initial condition of the simulations with the massive neutrinos were generated using ZA. To maintain consistency in computing the partial derivative, we made use of another 500 realisations of the fiducial cosmology (with M_ν = 0 eV) simulation in which the initial conditions were also generated using ZA.

4. Results

In Fig. 1, we show the direct measurements of ℛ(△₁₂₃) from the 15 000 reference halo catalogs (solid blue line) and compare it with the LO predictions (dashed orange line). We consider all triangular configurations with r_min ∈ (40, 45) and r_max ∈ (115, 120), and such that r₁₂ ≥ r₂₃ ≥ r₃₁. All the separation scales have a bin width of 5 h⁻¹ Mpc. It thus corresponds to a total of 766 triangular configurations, spanning from configuration ‘0’ being the smallest (i.e., △₁₂₃ ∈ {(40, 45),(40, 45),(40, 45)} h⁻¹ Mpc) to configuration ‘765’ being the largest (△₁₂₃ ∈ {(115, 120), (115, 120),(115, 120)} h⁻¹ Mpc). In Appendix A, we show the comparison between the theoretical predictions and the direct measurements when r₁₂ is fixed to (115, 120) h⁻¹ Mpc for the triangular configurations. We can see from Eqs. (8) and (9) that the LO prediction for the mean three-point relative velocities (R₁₂ and R₂₃) will be equal to each other when r₁₂ = r₂₃, irrespective of the length of the third side. This is directly visible in Fig. 1, where ℛ^theory = 1 when this condition is met. When comparing the theoretical predictions with the direct measurements from the halo catalogs, we see that it is generally accurate within 4–5% for configurations with all separation lengths greater than 55 h⁻¹ Mpc. As is expected when the separation length decreases, the fidelity of the LO prediction also decreases with the maximum deviation at about 27% for the triangular configuration {(100, 105),(50, 55),(50, 55)} h⁻¹Mpc. This motivates us to directly measure ℛ from the simulations to compute the derivatives for the Fisher information matrix.

Fig. 1. Comparison of theoretical prediction (orange dashed line) for ℛ(△123) using perturbation theory at LO against the direct measurement (blue solid line) from the halo catalogs of the Quijote suite of simulation is shown in the top panel. While the residual showing the deviation of the theoretical prediction from the direct measurement from the simulations is shown in the bottom panel. The blue shaded region denotes the 5% region. The triangle configuration go from the smallest being {(40, 45),(40, 45),(40, 45)} h−1 Mpc to the largest which corresponds to {(115, 120),(115, 120),(115, 120)} h−1 Mpc. The vertical dashed lines correspond to a change (increase) in the value of r12 bin in △123.

In addition, ℛ(△₁₂₃) is unaffected by variation in σ₈ as mentioned earlier in Sect. 2.2. To demonstrate this, we computed the ratio of ℛ(△₁₂₃) for ${\sigma }_8^{+}$ = 0.849 and ${\sigma }_8^{-}$ = 0.819 using dark matter particles from 30 realisations and showcase this result in Fig. 2. The (blue) dot represents the mean of the measurement, while the scatter is shown using the (orange) error bars. Thus we can conclude that ℛ is independent of σ₈. Similarly in Fig. 3, we take a look at the bias dependence of ℛ(△₁₂₃) (black solid line), and R₁₂ (blue dashed line), where we show the ratio of each between the halo and matter component for each of the summary statistics. The bias term for the mean relative velocity between pair ‘23’ in a triplet is similar to $R_{12}^{\mathrm{h}}$, and is thus not shown in the figure. As reported in Kuruvilla & Aghanim (2021), for these triangular configurations (assuming a scale independent bias), it yields a bias factor around 1.85 for R₁₂ and R₂₃. For the purpose of computing these ratios in the figure, we used the mean relative velocity information for matter from 30 realisations of the dark-matter only simulations; and for the halo, we utilised the 15 000 catalogues. The shaded regions represent the 1σ errors from the propagation of uncertainties of the mean relative velocity statistics for matter and the halo. We can see that on large separation scales, the newly introduced statistic (black solid line) is bias-independent, while for the smallest triangle configuration, there is a very weak dependence of bias when considering the newly introduced statistic. This thus supports the Ansatz presented in Eq. (12), where the LO in perturbation theory renders ℛ to be bias independent on linear scales. For all triangular configurations considered in this work, the bias is found to be equal to one within 1–2%; hence, for the purposes of Fisher matrix formalism, we consider ℛ being independent of a (constant) linear bias term.

Fig. 2. Ratio of ℛ at ${\sigma }_8^{+}$ = 0.849 to ℛ at ${\sigma }_8^{-}$ = 0.819. The (blue) solid dot, and the (orange) error bar represents the mean and the relative error on the mean, respectively.

Fig. 3. Dashed (blue) line: bias for R12, i.e. the mean radial relative velocity between pairs 1 and 2 in a triplet. The solid (black) line shows the (weak) bias dependence of the ratio statistics ℛ, which is equal to one within 1–2% for all triangular configurations considered here.

Since ℛ is found to be independent of σ₈ at the scales we are probing (i.e., r_min ≥ 40 h⁻¹ Mpc), this sets the statistic in the unique position of being unaffected by the degeneracy in the M_ν − σ₈ parameter plane. We checked the impact of summed neutrino mass on ℛ utilising three non-zero neutrino mass in Fig. 4. The solid (blue) line shows the impact of M_ν = 0.1 eV on ℛ when compared to zero neutrino-mass cosmology. Similarly, the dashed (orange) and dash-dotted (green) lines showcases the impact of M_ν = 0.2 and M_ν = 0.4 eV, respectively. As can be seen, when the neutrino mass increases there is a decrease in the infall velocity between a pair in most of the triangular configurations. This is related to the free-streaming of neutrinos as a result of them having large thermal velocities. As a result, below the free-streaming scale, neutrinos do not cluster, which further slows down the collapse of the matter overall. This leads to an overall reduction in the growth of overall density perturbations at scales below free-streaming scale, and thus causes a suppression of power on large Fourier modes when looking at the matter power spectrum (e.g., Wong 2011; Lesgourgues & Pastor 2012). When looking at ℛ for all the configurations we measured, the maximal effect of suppression is seen in the case when M_ν = 0.4 eV, and for the triangular configuration {(100, 105),(50, 55),(50, 55)} h⁻¹ Mpc when compared to the zero neutrino mass cosmology.

Fig. 4. Effect of summed neutrino mass on ℛ(△123), as measured directly from the simulations, when compared to zero neutrino mass fiducial cosmology. The summed neutrino mass considered here are denoted in the legend, and its units are given in eV.

4.1. Cosmological parameters

We now turn our attention to the information content in the new ratio statistic, ℛ, and we can see its viability in constraining the cosmological model. As discussed in Sect. 3.2, we achieve this using the Fisher information matrix, and the ingredients are the partial derivatives of ℛ with respect to the cosmological model parameters and the covariance matrix (C). We showcase the correlation matrix in Fig. 5, which is given as ${\mathrm{C}}_{\mathit{ij}}/\sqrt{{\mathrm{C}}_{\mathit{ii}}{\mathrm{C}}_{\mathit{jj}}}$, wherein the covariance matrix is directly measured from the simulations using 15 000 realisations. We notice the presence of non-diagonal terms within it being positively correlated at similar triangular configurations, while also tending to be negatively correlated, as the configurations differ substantially.

Fig. 5. Correlation matrix (i.e., the covariance matrix of ℛ normalised by its diagonal elements) computed using 15 000 realisations of the Quijote simulations. The triangle configurations are same as in Fig. 1.

We now take a look at the information content in ℛ using the Fisher information matrix formalism, as defined in Sect. 3.2. As mentioned, we computed both the elements of it directly from the simulations. Since the bias dependence of ℛ was shown to being very weak even at small scales (∼40–50 h⁻¹ Mpc) and independent at large scales (≥80 h⁻¹ Mpc), we did not consider the bias parameter in the Fisher-matrix formalism. Thus, the model parameters are Ω_m, Ω_b, h, n_s, and M_ν. We show the results of our Fisher forecast obtained from ℛ (orange colour) in Fig. 6, where the contour denote the 68.3% joint credible region for all possible model parameters. The 1σ marginalised error for any model parameter θ_α is given by $\sqrt{F_{\alpha \alpha }^{-1}}$, and they are as follows for the parameters we considered: {σ_Ωm, σ_Ωb, σ_h, σ_ns, σ_Mν} ≡ {0.0158, 0.0041, 0.0391, 0.0394, 0.1175} (where σ_Mν is given in units of eV).

Fig. 6. Joint 68.3% credible region for all the pairs of cosmological model parameters at z = 0. Different colours indicate the different summary statistics: $R_{12}^{\mathrm{h}}+R_{23}^{\mathrm{h}}$ (blue contour, taken from Kuruvilla & Aghanim 2021) and ℛ (orange contour). The fiducial values are indicated by the black dashed lines, and Mν is given in units of eV.

We compare these constraints with those obtained from the mean pairwise velocity, and the mean relative velocity between pairs in a triplet as reported in Kuruvilla & Aghanim (2021), shown in Fig. 6 using blue contour. As a fair comparison, we use the constraints obtained from them using r_min = 40 h⁻¹ Mpc. It presents a factor of improvement of {6.2, 7.6, 9.8, 12.9, 8.87} for {Ω_m, Ω_b, h, n_s, M_ν}, respectively, over the mean pairwise velocity. However, when comparing against the constraints from $R_{12}^{\mathrm{h}}({▵}_{123})+R_{23}^{\mathrm{h}}({▵}_{123})$ (i.e., the mean three-point relative velocities), ℛ has its constraining power reduced by a factor of 1.34–1.44 for all the model parameters. This shrinkage in constraining power was also seen when considering $R_{12}^{\mathrm{h}}$ and $R_{23}^{\mathrm{h}}$ separately as compared to its combination in Kuruvilla & Aghanim (2021).

To study the information content in ℛ(△₁₂₃), so far, we made use of all possible triangular configurations. However, we could also study how it changes if we were to omit certain triangular configurations. We first checked this assumption by omitting the equilateral and isosceles (where sides 12 and 23 are equal) triangles. This selection corresponds to a total of 630 triangular configuration, which equals to omitting 136 configurations from the total possible set. We find that for this selection, the information content of ℛ is reduced by 10% for all the cosmological parameters considered. We can make another selection criteria, namely: to omit configurations that have more than a 10% deviation from the theoretical predictions, as shown in Fig. 1. In this case, the selection leads to a total configuration of 754 (which implies that only 12 configurations were omitted from the total possible set). The reduction of the information content as such for this selection is very minimal compared to the total information content, whereby the constraints for all the cosmological parameters have been reduced by less than 1%.

4.2. Discussions

It is useful to consider how the constraints from ℛ fares against those obtained from clustering statistics. In order to answer that question, we compare the constraints obtained in this work with those obtained from the redshift space halo power spectrum and the halo bispectrum in Hahn et al. (2020). Compared with the constraints from the redshift-space power multipoles for k_max = 0.2 h Mpc⁻¹ (which is closest to the r_min considered in this work), ℛ obtains a factor of improvement of {2.3, 3.6, 4.5, 5.4, 5.7} for {Ω_m, Ω_b, h, n_s, M_ν}, respectively. However it is slightly reduced to {1.4, 2.9, 3.1, 3.3, 2.5} when comparing with the constraints from power spectrum when k_max = 0.5 h Mpc⁻¹. This improvement over the power spectrum (a two-point summary statistics) is not surprising as ℛ is based on the first moment of the three-point relative velocity statistics. Hence, it is interesting to compare the constraints against those obtained from the redshift-space bispectrum monopole, and when using all triangular configurations for k_max = 0.2 h Mpc⁻¹, ℛ still enjoys an improvement by a factor of {1.8, 2.9, 3.2, 3.1, 1.8} for {Ω_m, Ω_b, h, n_s, M_ν}, respectively. However, when considering a larger set of triangular configurations for the bispectrum with k_max = 0.5 h Mpc⁻¹, there is less constraining power for ℛ with a factor of {0.7, 1.0, 1.0, 0.9, 0.4} for {Ω_m, Ω_b, h, n_s, M_ν}, respectively. This is not surprising, as the bispectrum monopole in this case probes further into the nonlinear scales, whereas ℛ was analysed for triangular configurations with separation scales of 40 h⁻¹ Mpc and above.

As it currently stands, one of the limitations in applying the new statistic, ℛ(△₁₂₃), directly to any observational data is the lack of an estimator to measure the three-point mean radial relative velocities using the LOS velocities (as it will be the LOS velocity that can be measured from either the peculiar velocity surveys or kSZ experiments). In the case of the mean pairwise velocity, Ferreira et al. (1999) demonstrated how such an estimator can be constructed. Similarly for the case of R_ij(△₁₂₃), and ℛ(△₁₂₃), we will be constructing such an estimator in a future work. The application of a ratio statistic to observational data has its own complications. This could arise due to the noise in the data; additionally, in a survey we might find that R₂₃(△₁₂₃) to be zero for one of the triangular configuration bin, which could lead to a catastrophic failure of our ratio statistic. However, this would also be dealt with in the future when creating the estimator for ℛ.

In the case of the mean pairwise velocities, an alternative estimator exists, based on each tracers’ transverse velocity component (Yasini et al. 2019). On the other hand, the three-point mean relative velocity consists of a non-vanishing mean transverse component in the plane of the triangle (unlike in the case of the pairwise velocity, which has its transverse component equal to zero). Thus, we would construct an estimator for the non-vanishing three-point mean transverse relative velocity in a future work. Furthermore, the analysis we present here takes only the radial component into consideration and, hence, a combination of both radial and transverse components of the three-point mean relative velocity could further improve the chances of constraining the cosmological model with greater accuracy.

Another caveat that we do not discuss in this work is the mass dependence of the optical depth parameter, which has been shown to increase as the halo mass increases (e.g., Battaglia 2016). We considered it as an averaged quantity, as shown in Eqs. (2) and (10). However, we do not envision the mass dependence to affect ℛ(△₁₂₃), as long as the ratio is measured using the same mass bin. On the other hand, in assuming a fixed cosmology, we could consider a scenario where measuring $R_{23}^{\mathrm{h}}$ is fixed to a high halo mass bin while measuring $R_{12}^{\mathrm{h}}$ for various mass bins in Eq. (11). Thus, it could lead to making potential measurements of the (scaled) mass dependence of optical depth (and degenerate with the bias factor) from direct kSZ experiment directly.

5. Conclusions

The determination of neutrino mass using cosmological observables has become one of the main goals set by forthcoming cosmological surveys. However, two-point statistics are affected overall by the M_ν − σ₈ degeneracy, whether we are using clustering or relative velocity statistics (which limits the potential of constraining neutrino mass from cosmology). With regard to relative velocities, Kuruvilla & Porciani (2020) introduced the three-point mean relative statistics (i.e., the mean relative velocity between pairs in a triplet) and, subsequently, in Kuruvilla & Aghanim (2021), they quantified the cosmological information content within. It was found to offer substantial information gain when compared to two-point statistics (both power spectrum and mean pairwise velocity), while being competitive with the constraints from the bispectrum.

In this paper, we extend the application of the mean three-point relative velocity statistics and we introduce a new ratio statistic, ℛ, in Eq. (11), which is unaffected by σ₈. This enables us to constrain the neutrino mass, in addition to other cosmological parameters, independently of σ₈. Moreover, in the context of kSZ experiments, this statistic is independent of optical depth, hence it circumvents the optical depth degeneracy which currently acts a limiting factor in the determination of cosmological parameters from the kSZ experiments. Furthermore, the leading order perturbation theory prediction suggests that ℛ will be bias-independent on linear scales. We verified this assumption by measuring ℛ values for both halos and matter, finding that the bias is consistent with one at 1–2% for all triangular configurations we probed in this work (r_min = 40 h⁻¹ Mpc and r_max = 120 h⁻¹ Mpc).

We also studied the effect of summed neutrino mass on ℛ and found that as the neutrino mass increases the amplitude of ℛ decreases. This can be understood by the fact that due to the free streaming of neutrinos, the collapse of matter slows down; also, by the virtue that ℛ acts as a proxy to the mean infall velocity between pairs in a triplet, ℛ decreases as M_ν increases.

We used the Fisher-matrix formalism to quantify the information content in ℛ, where the necessary derivatives and the covariance matrices were directly measured from the Quijote suite of simulations. We utilised 15 000 realisations of the reference cosmology to compute the covariance matrix, and the partial derivatives were also computed directly from the simulations. We find that constraints obtained from ℛ show a factor of 6.2–12.9 in improvement when compared to the constraints obtained from the mean pairwise velocity. When compared to the power spectrum and bispectrum, it still achieves an improvement in the constraints by a factor of 2.3–5.7 and 1.8–3.2, respectively.

In summary, we have introduced a new statistic based on the mean radial relative velocity between pairs in a triplet and we have shown that it can act as robust cosmological observable that could lead to sizeable information gain in comparison to the mean radial pairwise velocity. One of the limitations of the kSZ experiments is the optical depth degeneracy and, thus, breaking this degeneracy requires some form of an external data set (for e.g., using fast radio bursts as suggested in Madhavacheril et al. 2019). This new statistic thus provides a way forward in which the cosmological parameters can be constrained using data from future kinetic Sunyaev–Zeldovich experiments alone, without being affected by the optical depth parameter.

Acknowledgments

We would like to thank the referee for their comments, Nabila Aghanim and Francisco Villaescusa-Navarro for useful discussions. JK acknowledges funding for the ByoPiC project from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program grant agreement ERC-2015-AdG 695561. We are thankful to the community for developing and maintaining open-source software packages extensively used in our work, namely CYTHON (Behnel et al. 2011), MATPLOTLIB (Hunter 2007) and NUMPY (Harris et al. 2020).

References

Appendix A: Shape dependence of ℛ

We compare the direct measurement of ℛ from the halo catalogs (left panel) and the theoretical predictions (right panel) for a small subset of triangular configurations in Fig. A.1. In this case, we fix one side of the triangle, namely, r₁₂, and vary the other two sides. The top right bin in both panels (r₂₃/r₁₂ = r₃₁/r₁₂ = 1) corresponds to the equilateral configuration. In Fig. A.2, we show the ratio between them and see that the largest discrepancy can be seen at r₃₁/r₁₂ ∼ 0.36 and r₂₃/r₁₂ ∼ 0.68. This is expected as the separation decreases the fidelity of the LO prediction also decreases.

Fig. A.1. Comparison of ℛ measured directly from the halo catalogs (left panel), and the LO analytical prediction (right panel). The triangular configurations corresponds to the set where one of the leg of the triangle (r12) is fixed to (115, 120) h−1 Mpc.

Fig. A.2. Ratio between ℛ measured directly from the halo catalogs, and the LO analytical prediction. The triangular configurations are the same as in Fig. A.1.

All Figures

Fig. 1. Comparison of theoretical prediction (orange dashed line) for ℛ(△123) using perturbation theory at LO against the direct measurement (blue solid line) from the halo catalogs of the Quijote suite of simulation is shown in the top panel. While the residual showing the deviation of the theoretical prediction from the direct measurement from the simulations is shown in the bottom panel. The blue shaded region denotes the 5% region. The triangle configuration go from the smallest being {(40, 45),(40, 45),(40, 45)} h−1 Mpc to the largest which corresponds to {(115, 120),(115, 120),(115, 120)} h−1 Mpc. The vertical dashed lines correspond to a change (increase) in the value of r12 bin in △123.

In the text

	Fig. 2. Ratio of ℛ at ${\sigma }_8^{+}$ = 0.849 to ℛ at ${\sigma }_8^{-}$ = 0.819. The (blue) solid dot, and the (orange) error bar represents the mean and the relative error on the mean, respectively.
In the text

	Fig. 3. Dashed (blue) line: bias for R12, i.e. the mean radial relative velocity between pairs 1 and 2 in a triplet. The solid (black) line shows the (weak) bias dependence of the ratio statistics ℛ, which is equal to one within 1–2% for all triangular configurations considered here.
In the text

	Fig. 4. Effect of summed neutrino mass on ℛ(△123), as measured directly from the simulations, when compared to zero neutrino mass fiducial cosmology. The summed neutrino mass considered here are denoted in the legend, and its units are given in eV.
In the text

	Fig. 5. Correlation matrix (i.e., the covariance matrix of ℛ normalised by its diagonal elements) computed using 15 000 realisations of the Quijote simulations. The triangle configurations are same as in Fig. 1.
In the text

	Fig. 6. Joint 68.3% credible region for all the pairs of cosmological model parameters at z = 0. Different colours indicate the different summary statistics: $R_{12}^{\mathrm{h}}+R_{23}^{\mathrm{h}}$ (blue contour, taken from Kuruvilla & Aghanim 2021) and ℛ (orange contour). The fiducial values are indicated by the black dashed lines, and Mν is given in units of eV.
In the text

	Fig. A.1. Comparison of ℛ measured directly from the halo catalogs (left panel), and the LO analytical prediction (right panel). The triangular configurations corresponds to the set where one of the leg of the triangle (r12) is fixed to (115, 120) h−1 Mpc.
In the text

	Fig. A.2. Ratio between ℛ measured directly from the halo catalogs, and the LO analytical prediction. The triangular configurations are the same as in Fig. A.1.
In the text