Issue 
A&A
Volume 660, April 2022



Article Number  A113  
Number of page(s)  9  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/202142325  
Published online  22 April 2022 
Cosmology with the kinetic Sunyaev–Zeldovich effect: Independent of the optical depth and σ_{8}
Université ParisSaclay, CNRS, Institut d’Astrophysique Spatiale, 91405 Orsay, France
email: joseph.kuruvilla@universiteparissaclay.fr
Received:
28
September
2021
Accepted:
3
February
2022
Cosmological constraints obtained by the kinetic Sunyaev–Zeldovich experiments are degenerate with the optical depth measurement – an effect that is commonly known as the opticaldepth degeneracy. In this work, we introduce a new statistic based on the first moment of relative velocity between pairs in a triplet, which is capable of constraining cosmological parameters independently of the optical depth and of σ_{8}. Using 22 000 Nbody simulations from the Quijote suite, we quantified the information content in the new statistic using Fisher matrix forecast. We find that it is able to obtain strong constraints on the cosmological parameters, particularly on the summed neutrino mass. The constraints bring an improvement on all cosmological model parameters by a factor of 6.2–12.9 and 2.3–5.7 when compared to those obtained from the mean pairwise velocity and from the redshiftspace halo power spectrum, respectively. Thus, this new statistic paves a way forward in constraining cosmological parameters independent of the optical depth and σ_{8} when using data from future kinetic Sunyaev–Zeldovich experiments alone.
Key words: largescale structure of Universe / cosmology: theory / cosmic background radiation
© J. Kuruvilla 2022
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
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 nonzero mass (e.g., Forero et al. 2014; GonzalezGarcia 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 subeV 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 redshiftspace 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^{1} (SO, Ade et al. 2019), CMBS4^{2} (Abazajian et al. 2016)) as well as the stage IV galaxy redshift surveys (e.g., the ‘Dark Energy Spectroscopic Instrument’^{3} (DESI, DESI Collaboration 2016), Euclid^{4} (Laureijs et al. 2011) and Nancy Grace Roman space telescope^{5} (Spergel et al. 2015)).
Currently, the M_{ν} constraints from galaxy clustering are mainly obtained using twopoint statistics, namely, the power spectrum in the Fourier space or the twopoint correlation function in the configuration space. The impact of massive neutrinos on the twopoint clustering statistics has been studied quite extensively using Nbody simulations both in the real (e.g., Saito et al. 2008; Wong 2008; Castorina et al. 2015) and the redshift space (e.g., VillaescusaNavarro et al. 2018; GarcíaFarieta et al. 2019). However, these statistics are affected by the M_{ν} − σ_{8} degeneracy, thus acting as a limitation in measuring the summed neutrino mass. The threepoint clustering statistics in the Fourier space (i.e., the bispectrum) has been shown to break this degeneracy (Hahn et al. 2020; Hahn & VillaescusaNavarro 2021). In addition, it has been shown that the threepoint cluster statistics contain additional cosmological information compared to its twopoint 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 onepoint 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 threepoint 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 threepoint 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 threepoint 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 Fishermatrix 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
where σ_{T} is the Thomson scattering crosssection, 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 . The optical depth is defined as τ = ∫dlσ_{T}n_{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ándezMonteagudo 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ándezMonteagudo et al. 2015; Schaan et al. 2016, 2021; Nguyen et al. 2020; Tanimura et al. 2021; ChavesMontero 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)
where is the mean temperature difference between the objects ‘1’ and ‘2’, and is the mean radial component of the pairwise velocity that can be defined in the single streaming regime as follows:
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)
where is the unit vector along the pair ‘12’, the subscript p implies that the averages are computed over all pairs with separation r_{12}, P(k) denotes the linear matter power spectrum, and j_{1}(x) = sin(x)/x^{2} − 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
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 △_{123} = (r_{12}, r_{23}, r_{31}): (i) ⟨w_{12}△_{123}⟩_{t}, and (ii) ⟨w_{23}△_{123}⟩_{t}. The subscript t here implies that the averages are computed over all triplets with separations (r_{12}, r_{23}, r_{31}). 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):
The threepoint 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_{12}△_{123}⟩_{t}, it follows as:
where , is the unit vector along the pair ‘23’, and . 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
Similarly, the mean radial relative velocity between the pair 23 in △_{123} can be written as
Similarly to Eq. (2), the threepoint mean relative temperature difference from kSZ can be expressed as:
where respresents the threepoint mean relative velocity statistics for haloes (biased tracers), and to first approximation it can be written down as linear bias term times R_{ij}(△_{123}) (Kuruvilla & Aghanim 2021). Based on the radial mean relative velocities between pairs in a triplet, we can introduce a new ratio statistic, namely:
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 △_{123}. On linear scales, using perturbation theory at LO, ℛ(△_{123}) can be written as the ratio between Eqs. (8) and (9)
The above statistic is thus independent of optical depth, σ_{8}, and linear bias. In the following sections, we take a detailed look at whether the Ansatz of σ_{8} and linear bias independence holds up. Additionally, we study the cosmological information content in ℛ(△_{123}).
3. Data and analysis
3.1. Quijote simulation suite
In this work, we make use of the Quijote^{6} (VillaescusaNavarro et al. 2020) suite of simulations, which was run using the treePM code GADGET3 (Springel 2005). Spanning more than a few thousand cosmological models, it contains 44 100 Nbody simulations. These simulations have a box length of 1 h^{−1} Gpc, tracking the evolution of 512^{3} cold dark matter (CDM) particles. The initial conditions (ICs) were generated at redshift z = 127 using the secondorder Lagrangian perturbation theory. The fiducial cosmological parameters (assuming zerosummed 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^{−1} Mpc: σ_{8} = 0.834; and the presentday value of the Hubble constant: H_{0} ≡ H(z = 0) = 100 h km s^{−1} Mpc^{−1}, 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: 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^{3} 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 Nbody simulations of the Quijote suite. These halos were identified using a friendsoffriends algorithm. We selected halos that have a halo mass of M_{h} > 5 × 10^{13}h^{−1} M_{⊙} (corresponding to groups and clusters of galaxies) at z = 0, which gives a mean number density of for the reference simulations. Additionally, in the case of (i) fiducial cosmology, and (ii) for variations in σ_{8} (both and ), we used 30 realisations of the particle data (randomly downsampled to 100^{3} particles) to compute ℛ(△_{123}).
3.2. Fishermatrix formalism
To quantify the error estimates on the cosmological parameters, we used the Fishermatrix formalism which can be defined as (e.g., Tegmark et al. 1997; Heavens 2009; Verde 2010):
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
where ℛ represents the data vector for the ratio statistic we introduced in Eq. (11) and 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:
where , 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):
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}, σ_{8}}, we can make use of the central difference approximation to compute the derivative numerically, namely:
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
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 ℛ(△_{123}) 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_{12} ≥ r_{23} ≥ r_{31}. All the separation scales have a bin width of 5 h^{−1} Mpc. It thus corresponds to a total of 766 triangular configurations, spanning from configuration ‘0’ being the smallest (i.e., △_{123} ∈ {(40, 45),(40, 45),(40, 45)} h^{−1} Mpc) to configuration ‘765’ being the largest (△_{123} ∈ {(115, 120), (115, 120),(115, 120)} h^{−1} Mpc). In Appendix A, we show the comparison between the theoretical predictions and the direct measurements when r_{12} is fixed to (115, 120) h^{−1} Mpc for the triangular configurations. We can see from Eqs. (8) and (9) that the LO prediction for the mean threepoint relative velocities (R_{12} and R_{23}) will be equal to each other when r_{12} = r_{23}, 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^{−1} 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^{−1}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 r_{12} bin in △_{123}. 
In addition, ℛ(△_{123}) is unaffected by variation in σ_{8} as mentioned earlier in Sect. 2.2. To demonstrate this, we computed the ratio of ℛ(△_{123}) for = 0.849 and = 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 σ_{8}. Similarly in Fig. 3, we take a look at the bias dependence of ℛ(△_{123}) (black solid line), and R_{12} (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 , 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_{12} and R_{23}. For the purpose of computing these ratios in the figure, we used the mean relative velocity information for matter from 30 realisations of the darkmatter 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 biasindependent, 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 = 0.849 to ℛ at = 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 R_{12}, 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 σ_{8} at the scales we are probing (i.e., r_{min} ≥ 40 h^{−1} Mpc), this sets the statistic in the unique position of being unaffected by the degeneracy in the M_{ν} − σ_{8} parameter plane. We checked the impact of summed neutrino mass on ℛ utilising three nonzero neutrino mass in Fig. 4. The solid (blue) line shows the impact of M_{ν} = 0.1 eV on ℛ when compared to zero neutrinomass cosmology. Similarly, the dashed (orange) and dashdotted (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 freestreaming of neutrinos as a result of them having large thermal velocities. As a result, below the freestreaming 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 freestreaming 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^{−1} 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 , wherein the covariance matrix is directly measured from the simulations using 15 000 realisations. We notice the presence of nondiagonal 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^{−1} Mpc) and independent at large scales (≥80 h^{−1} Mpc), we did not consider the bias parameter in the Fishermatrix 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 , 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: (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^{−1} 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 (i.e., the mean threepoint 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 and separately as compared to its combination in Kuruvilla & Aghanim (2021).
To study the information content in ℛ(△_{123}), 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 redshiftspace power multipoles for k_{max} = 0.2 h Mpc^{−1} (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^{−1}. This improvement over the power spectrum (a twopoint summary statistics) is not surprising as ℛ is based on the first moment of the threepoint relative velocity statistics. Hence, it is interesting to compare the constraints against those obtained from the redshiftspace bispectrum monopole, and when using all triangular configurations for k_{max} = 0.2 h Mpc^{−1}, ℛ 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^{−1}, 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^{−1} Mpc and above.
As it currently stands, one of the limitations in applying the new statistic, ℛ(△_{123}), directly to any observational data is the lack of an estimator to measure the threepoint 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}(△_{123}), and ℛ(△_{123}), 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_{23}(△_{123}) 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 threepoint mean relative velocity consists of a nonvanishing 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 nonvanishing threepoint 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 threepoint 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 ℛ(△_{123}), 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 is fixed to a high halo mass bin while measuring 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, twopoint statistics are affected overall by the M_{ν} − σ_{8} 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 threepoint 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 twopoint 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 threepoint relative velocity statistics and we introduce a new ratio statistic, ℛ, in Eq. (11), which is unaffected by σ_{8}. This enables us to constrain the neutrino mass, in addition to other cosmological parameters, independently of σ_{8}. 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 biasindependent 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^{−1} Mpc and r_{max} = 120 h^{−1} 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 Fishermatrix 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 VillaescusaNavarro 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 ERC2015AdG 695561. We are thankful to the community for developing and maintaining opensource software packages extensively used in our work, namely CYTHON (Behnel et al. 2011), MATPLOTLIB (Hunter 2007) and NUMPY (Harris et al. 2020).
References
 Abazajian, K. N., Adshead, P., Ahmed, Z., et al. 2016, ArXiv eprints [arXiv:1610.02743] [Google Scholar]
 Ade, P., Aguirre, J., Ahmed, Z., et al. 2019, J. Cosmol. Astropart. Phys., 2019, 056 [Google Scholar]
 Agarwal, N., Desjacques, V., Jeong, D., & Schmidt, F. 2021, J. Cosmol. Astropart. Phys., 2021, 021 [CrossRef] [Google Scholar]
 Aker, M., Beglarian, A., Behrens, J., et al. 2021, ArXiv eprints [arXiv:2105.08533] [Google Scholar]
 Alam, S., Aubert, M., Avila, S., et al. 2021, Phys. Rev. D, 103, 083533 [NASA ADS] [CrossRef] [Google Scholar]
 Anderson, T. 2003, An Introduction to Multivariate Statistical Analysis, Wiley Series in Probability and Statistics (New York: Wiley) [Google Scholar]
 Battaglia, N. 2016, J. Cosmol. Astropart. Phys., 2016, 058 [CrossRef] [Google Scholar]
 Bayer, A. E., VillaescusaNavarro, F., Massara, E., et al. 2021, ApJ, 919, 24 [NASA ADS] [CrossRef] [Google Scholar]
 Behnel, S., Bradshaw, R., Citro, C., et al. 2011, Comput. Sci. Eng., 13, 13 [Google Scholar]
 Calafut, V., Gallardo, P. A., Vavagiakis, E. M., et al. 2021, Phys. Rev. D, 104, 043502 [NASA ADS] [CrossRef] [Google Scholar]
 Capozzi, F., Di Valentino, E., Lisi, E., et al. 2017, Phys. Rev. D, 95, 096014 [NASA ADS] [CrossRef] [Google Scholar]
 Castorina, E., Carbone, C., Bel, J., Sefusatti, E., & Dolag, K. 2015, J. Cosmol. Astropart. Phys., 2015, 043 [CrossRef] [Google Scholar]
 ChavesMontero, J., HernándezMonteagudo, C., Angulo, R. E., & Emberson, J. D. 2021, MNRAS, 503, 1798 [NASA ADS] [CrossRef] [Google Scholar]
 Chen, Z., Zhang, P., Yang, X., & Zheng, Y. 2022, MNRAS, 510, 5916 [NASA ADS] [CrossRef] [Google Scholar]
 Choudhury, S. R., & Hannestad, S. 2020, J. Cosmol. Astropart. Phys., 2020, 037 [CrossRef] [Google Scholar]
 Chudaykin, A., & Ivanov, M. M. 2019, J. Cosmol. Astropart. Phys., 2019, 034 [CrossRef] [Google Scholar]
 De Bernardis, F., Aiola, S., Vavagiakis, E. M., et al. 2017, J. Cosmol. Astropart. Phys., 3, 008 [Google Scholar]
 de Salas, P. F., Forero, D. V., Ternes, C. A., Tortola, M., & Valle, J. W. F. 2017, ArXiv eprints [arXiv:1708.01186] [Google Scholar]
 DESI Collaboration (Aghamousa, A., et al.) 2016, ArXiv eprints [arXiv:1611.00036] [Google Scholar]
 Ferreira, P. G., Juszkiewicz, R., Feldman, H. A., Davis, M., & Jaffe, A. H. 1999, ApJ, 515, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Fisher, K. B. 1995, ApJ, 448, 494 [NASA ADS] [CrossRef] [Google Scholar]
 Flender, S., Nagai, D., & McDonald, M. 2017, ApJ, 837, 124 [NASA ADS] [CrossRef] [Google Scholar]
 Forero, D. V., Tórtola, M., & Valle, J. W. F. 2014, Phys. Rev. D, 90, 093006 [NASA ADS] [CrossRef] [Google Scholar]
 GarcíaFarieta, J. E., Marulli, F., Veropalumbo, A., et al. 2019, MNRAS, 488, 1987 [Google Scholar]
 GonzalezGarcia, M. C., Maltoni, M., & Schwetz, T. 2016, Nucl. Phys. B, 908, 199 [NASA ADS] [CrossRef] [Google Scholar]
 Gualdi, D., & Verde, L. 2020, J. Cosmol. Astropart. Phys., 2020, 041 [CrossRef] [Google Scholar]
 Hahn, C., & VillaescusaNavarro, F. 2021, J. Cosmol. Astropart. Phys., 2021, 029 [CrossRef] [Google Scholar]
 Hahn, C., VillaescusaNavarro, F., Castorina, E., & Scoccimarro, R. 2020, J. Cosmology Astropart. Phys., 2020, 040 [CrossRef] [Google Scholar]
 Hand, N., Addison, G. E., Aubourg, E., et al. 2012, Phys. Rev. Lett., 109, 041101 [NASA ADS] [CrossRef] [Google Scholar]
 Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 585 [Google Scholar]
 Hartlap, J., Simon, P., & Schneider, P. 2007, A&A, 464, 399 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Heavens, A. 2009, ArXiv eprints [arXiv:0906.0664] [Google Scholar]
 HernándezMonteagudo, C., Ma, Y.Z., Kitaura, F. S., et al. 2015, Phys. Rev. Lett., 115, 191301 [Google Scholar]
 Hunter, J. D. 2007, Comput. Sci. Eng., 9, 9 [Google Scholar]
 Juszkiewicz, R., Fisher, K. B., & Szapudi, I. 1998, ApJ, 504, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Kaufmann, G. M. 1967, Some Bayesian Moment Formulae, Report No. 6710. Centre for Operations Research and Econometrics (Heverlee: Catholic University of Louvain) [Google Scholar]
 Keisler, R., & Schmidt, F. 2013, ApJ, 765, L32 [NASA ADS] [CrossRef] [Google Scholar]
 Kodwani, D., Alonso, D., & Ferreira, P. 2019, Open J. Astrophys., 2, 3 [NASA ADS] [CrossRef] [Google Scholar]
 Kuruvilla, J., & Aghanim, N. 2021, A&A, 653, A130 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kuruvilla, J., & Porciani, C. 2020, J. Cosmol. Astropart. Phys., 2020, 043 [CrossRef] [Google Scholar]
 Kuruvilla, J., Aghanim, N., & McCarthy, I. G. 2020, A&A, 644, A170 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, ArXiv eprints [arXiv:1110.3193] [Google Scholar]
 Lesgourgues, J., & Pastor, S. 2012, Adv. High Energy Phys., 2012, 608515 [CrossRef] [Google Scholar]
 Li, Y.C., Ma, Y.Z., Remazeilles, M., & Moodley, K. 2018, Phys. Rev. D, 97, 023514 [NASA ADS] [CrossRef] [Google Scholar]
 Madhavacheril, M. S., Battaglia, N., Smith, K. M., & Sievers, J. L. 2019, Phys. Rev. D, 100, 103532 [NASA ADS] [CrossRef] [Google Scholar]
 Mueller, E.M., de Bernardis, F., Bean, R., & Niemack, M. D. 2015, Phys. Rev. D, 92, 063501 [NASA ADS] [CrossRef] [Google Scholar]
 Nguyen, N.M., Jasche, J., Lavaux, G., & Schmidt, F. 2020, J. Cosmol. Astropart. Phys., 2020, 011 [CrossRef] [Google Scholar]
 Planck Collaboration VI. 2020, A&A, 641, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XXXVII. 2016, A&A, 586, A140 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Reid, B. A., & White, M. 2011, MNRAS, 417, 1913 [NASA ADS] [CrossRef] [Google Scholar]
 Saito, S., Takada, M., & Taruya, A. 2008, Phys. Rev. Lett., 100, 191301 [NASA ADS] [CrossRef] [Google Scholar]
 Samushia, L., Slepian, Z., & VillaescusaNavarro, F. 2021, MNRAS, 505, 628 [NASA ADS] [CrossRef] [Google Scholar]
 Schaan, E., Ferraro, S., VargasMagaña, M., et al. 2016, Phys. Rev. D, 93, 082002 [NASA ADS] [CrossRef] [Google Scholar]
 Schaan, E., Ferraro, S., Amodeo, S., et al. 2021, Phys. Rev. D, 103, 063513 [NASA ADS] [CrossRef] [Google Scholar]
 Smith, K. M., Madhavacheril, M. S., Münchmeyer, M., et al. 2018, ArXiv eprints [arXiv:1810.13423] [Google Scholar]
 Soergel, B., Flender, S., Story, K. T., et al. 2016, MNRAS, 461, 3172 [Google Scholar]
 Soergel, B., Saro, A., Giannantonio, T., Efstathiou, G., & Dolag, K. 2018, MNRAS, 478, 5320 [NASA ADS] [CrossRef] [Google Scholar]
 Spergel, D., Gehrels, N., Baltay, C., et al. 2015, ArXiv eprints [arXiv:1503.03757] [Google Scholar]
 Springel, V. 2005, MNRAS, 364, 1105 [Google Scholar]
 Sunyaev, R. A., & Zeldovich, I. B. 1980, MNRAS, 190, 413 [NASA ADS] [CrossRef] [Google Scholar]
 Sunyaev, R. A., & Zeldovich, Y. B. 1972, Comments Astrophys. Space Phys., 4, 173 [NASA ADS] [EDP Sciences] [Google Scholar]
 Tanimura, H., Zaroubi, S., & Aghanim, N. 2021, A&A, 645, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tegmark, M., Taylor, A. N., & Heavens, A. F. 1997, ApJ, 480, 22 [NASA ADS] [CrossRef] [Google Scholar]
 Uhlemann, C., Friedrich, O., VillaescusaNavarro, F., Banerjee, A., & Codis, S. 2020, MNRAS, 495, 4006 [NASA ADS] [CrossRef] [Google Scholar]
 Vagnozzi, S., Dhawan, S., Gerbino, M., et al. 2018, Phys. Rev. D, 98, 083501 [NASA ADS] [CrossRef] [Google Scholar]
 Verde, L. 2010, Statistical Methods in Cosmology (Berlin: Springer Verlag), 800, 147 [NASA ADS] [Google Scholar]
 VillaescusaNavarro, F., Banerjee, A., Dalal, N., et al. 2018, ApJ, 861, 53 [NASA ADS] [CrossRef] [Google Scholar]
 VillaescusaNavarro, F., Hahn, C., Massara, E., et al. 2020, ApJS, 250, 2 [CrossRef] [Google Scholar]
 Wong, Y. Y. Y. 2008, J. Cosmol. Astropart. Phys., 2008, 035 [CrossRef] [Google Scholar]
 Wong, Y. Y. Y. 2011, Ann. Rev. Nucl. Part. Sci., 61, 69 [NASA ADS] [CrossRef] [Google Scholar]
 Yankelevich, V., & Porciani, C. 2019, MNRAS, 483, 2078 [NASA ADS] [CrossRef] [Google Scholar]
 Yasini, S., Mirzatuny, N., & Pierpaoli, E. 2019, ApJ, 873, L23 [NASA ADS] [CrossRef] [Google Scholar]
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_{12}, and vary the other two sides. The top right bin in both panels (r_{23}/r_{12} = r_{31}/r_{12} = 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_{31}/r_{12} ∼ 0.36 and r_{23}/r_{12} ∼ 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 (r_{12}) 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 r_{12} bin in △_{123}. 

In the text 
Fig. 2. Ratio of ℛ at = 0.849 to ℛ at = 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 R_{12}, 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: (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 (r_{12}) 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.