Information content in mean pairwise velocity and mean relative velocity between pairs in a triplet

Velocity field provides a complementary avenue to constrain cosmological information, either through the peculiar velocity surveys or the kinetic Sunyaev Zel'dovich effect. One of the commonly used statistics is the mean radial pairwise velocity. Here, we consider the three-point mean relative velocity, i.e. the mean relative velocities between pairs in a triplet. Using halo catalogs from the Quijote suite of N-body simulations, we first showcase how the analytical prediction for the mean relative velocities between pairs in a triplet achieve better than 4-5% accuracy using standard perturbation theory at leading order for triangular configurations with a minimum separation of $r \geq 50\ h^{-1}$Mpc. Furthermore, we present the three-point relative velocity as a novel probe of neutrino mass estimation. We explore the full cosmological information content of the halo mean pairwise velocities, and the mean relative velocities between halo pairs in a triplet. We undertake this through the Fisher-matrix formalism using 22,000 simulations from the Quijote suite, and considering all triangular configurations with a minimum and a maximum separation of $20\ h^{-1}$Mpc and $120\ h^{-1}$Mpc, respectively. We find that the mean relative velocities in a triplet allows a 1$\sigma$ neutrino mass ($M_\nu$) constraint of 0.065 eV, that is roughly 13 times better than the mean pairwise velocity constraint (0.877 eV). This information gain is not limited only to neutrino mass, but extends to other cosmological parameters: $\Omega_{\mathrm{m}}$, $\Omega_{\mathrm{b}}$, $h$, $n_{\mathrm{s}}$ and $\sigma_{8}$ achieving a gain of 8.9, 11.8, 15.5, 20.9 and 10.9 times respectively. These results illustrate the possibility of exploiting the mean three-point relative velocities for constraining the cosmological parameters accurately from future cosmic microwave background experiments and peculiar velocity surveys.


Introduction
The last decades have witnessed a tremendous increase in our understanding of the underlying cosmological model. This has been mainly facilitated by cosmic microwave background (CMB) experiments (e.g. Planck Collaboration 2020), galaxy redshift surveys (e.g. eBOSS Collaboration 2020), and gravitational lensing surveys (e.g. Heymans et al. 2021). However, there are several key areas in which our understanding is incomplete or lacking. The questions regarding the nature of dark matter and dark energy and the summed neutrino mass (M ν ) are the two main scientific drivers for forthcoming redshift surveys such as Euclid 1 and future CMB experiments such as the Simons Observatory 2 (SO, Ade et al. 2019) and CMB-S4 3 (Abazajian et al. 2016).
Cosmological information from galaxy redshift surveys is often extracted using the two-point correlation function in the configuration space or its Fourier analogue, the power spectrum. However, as the complex galaxy distribution represents a non-Gaussian field, higher order statistics such as the bispectrum and the three-point correlation function contain additional information. Thus, by considering both two-point and three-point clustering information, future redshift surveys will be able to obtain a substantial improvement with regard to the cosmological parameter constraints (e.g. Yankelevich Agarwal et al. 2021;Samushia et al. 2021). The determination of the summed mass of neutrinos will be one of the primary goals of future redshift surveys. This will be possible as the free-streaming of neutrinos imprint unique signatures on galaxy clustering information in both real and redshift space (e.g. Saito et al. 2008;Wong 2008;Castorina et al. 2015;Villaescusa-Navarro et al. 2018;García-Farieta et al. 2019;Kamalinejad & Slepian 2020).
Complementary to clustering statistics based on the density field is the study of summary statistics based on the peculiar velocity field. We can observe the velocity fields via: (i) peculiar velocity surveys, and (ii) the kinetic Sunyaev-Zel'dovich (kSZ) effect (Sunyaev & Zeldovich 1972, 1980. Peculiar velocity can be directly measured using empirical relations like the Tully-Fisher (Tully & Fisher 1977) relation and the fundamental plane (Djorgovski & Davis 1987;Dressler et al. 1987) relation. These surveys offer an opportunity to probe the peculiar velocities in the nearby Universe. In particular, future peculiar velocity surveys like the Taipan galaxy survey 4 (da Cunha et al. 2017), and the Widefield ASKAP L-band Legacy All-sky Blind Survey 5 (WALLABY, Koribalski et al. 2020) have been shown to constrain the growth rate of structure rather competitively with the current galaxy redshift surveys at low redshift (Koda et al. 2014;Howlett et al. 2017). Recently, Dupuy et al. (2019) measured the growth rate of the structure by applying the mean pairwise velocity estimator on the Cosmicflows-3 dataset (Tully et al. 2016).
The kSZ effect is a secondary anisotropy where CMB photons are scattered off free electrons in motion due to the bulk motion of the ionised medium; e.g. cluster peculiar motion. It results in a Doppler shift, which preserves the blackbody spectrum of the CMB. The first significant detection of the kSZ effect was achieved through the mean pairwise velocity statistics (Hand et al. 2012) using an estimator developed in Ferreira et al. (1999). Further detection of the kSZ effect using the mean pairwise velocities were presented in Hernández-Monteagudo et al. 2015, Planck Collaboration 2016, Soergel et al. 2016, De Bernardis et al. 2017, Li et al. 2018, and Calafut et al. 2021. The mean radial pairwise velocity has been shown to be capable of constraining alternative theories of gravity and dark energy (Bhattacharya & Kosowsky 2007, 2008Kosowsky & Bhattacharya 2009;Mueller et al. 2015a). There have been other methods via which the kSZ effect has been detected; for example, (i) by correlating CMB maps with reconstructed velocity field (e.g. Hernández-Monteagudo et al. 2015;Schaan et al. 2016;Nguyen et al. 2020;Tanimura et al. 2021), and (ii) by cross correlating CMB maps with angular redshift fluctuation maps (Chaves-Montero et al. 2021). The forthcoming CMB surveys including SO, CMB-S4, and CMB-HD (Sehgal et al. 2019a,b) will offer a plethora of data and allow us to measure kSZ effect accurately.
In addition to their imprint on galaxy clustering, massive neutrinos also leave their imprint on various velocity statistics. For the pairwise velocity statistics, the impact of massive neutrinos and its interplay with baryonic physics at non-linear scales was studied in  using the BAHAMAS suite of simulations (McCarthy et al. 2018). It was found that the mean matter pairwise velocity decreases with increasing neutrino mass at those pair separations. For analytical predictions, Aviles & Banerjee (2020) studied the effect of massive neutrinos on pairwise velocity using Lagrangian perturbation theory for biased tracers at quasi-linear scales and above. The pairwise velocity measurement from kSZ experiments has been shown to be a novel probe to constrain the summed neutrino mass (Mueller et al. 2015b), and provide complementary constraints with respect to galaxy clustering and CMB experiments. In line with this, we want to explore how much information the higher-order relative velocity statistics provide on the cosmological parameters, in particular with regard to summed neutrino mass. The three-point generalisation of the pairwise velocity, i.e. the relative velocities between pairs in a triplet, was recently introduced in Kuruvilla & Porciani (2020). However, the impact of cosmological parameters (especially massive neutrinos) on it remains unexplored. Thus, in this paper, we study how the different cosmological parameters (including M ν ) affect the relative velocities between pairs in a triplet and this constitutes one of the main aims of this work.
We used the Fisher-matrix formalism to quantify the viability of using the mean relative velocities between pairs in a triplet to constrain cosmological parameters, and we compared it with the cosmological information obtained from the mean pairwise velocity statistics. For this purpose, we measured the necessary derivatives and the covariance matrix directly from the Quijote suite of simulations ). Thus, this work forms a parallel line of study to Hahn et al. (2020) and Hahn & Villaescusa-Navarro (2021), which quantified the information content in the redshift-space bispectrum of halos and galaxies, respectively, using the Fisher-matrix formalism; moreover, they computed the necessary quantities directly from the Quijote simulations. It was shown that the redshift-space bispectrum was able to achieve sizeable gains in the cosmological parameter constraints over those obtained from the power spectrum; notably, it was especially able to constrain summed neutrino mass to high precision. We aim to study if the mean relative velocity statistics can provide competitive constraints when compared to those from the clustering analyses.
The paper is structured as follows. The Quijote suite of simulation, which we used in this work, is introduced briefly in Sect. 2.1. The information content of the velocity statistics was studied using the Fisher-matrix formalism, which is briefly summarised in Sect. 2.2. In Sect. 3, we describe the summary statistics based on velocity that we employed in this work. Our results are presented in Sect. 4 and finally concluded in Sect. 5.

Quijote simulation suite
The Quijote 6 (Villaescusa-Navarro et al. 2020) simulation suite contains 44,100 N-body simulations spanning more than a few thousand cosmological models, run using the tree-PM code gadget-3 (Springel 2005). The initial conditions (ICs) of the simulation, at redshift z = 127, were generated using the second-order Lagrangian perturbation theory (for runs with zero mass neutrino cosmology) and the Zel'dovich approximation (for simulations with massive neutrinos). Each of these simulations spans a volume of 1 h −3 Gpc 3 and follows the cosmological evolution of 512 3 cold dark matter (CDM) particles (and, additionally, in the case of massive neutrino simulations, 512 3 neutrino particles). The fiducial cosmological parameters for the simulation are consistent with the Planck Collaboration 2020, and these are as follows: (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 (present-day value of the Hubble constant) H 0 ≡ H(z = 0) = 100 h km s −1 Mpc −1 with h = 0.6711. However, it assumes neutrinos to be massless, i.e. M ν = 0.0 eV. The reference run with the fiducial cosmological model consists of 15,000 random realisations. The Quijote suite of simulations is thus ideal for constructing an accurate covariance matrix of any cosmological observable. Similarly, it provides a set of 500 random realisations where only one cosmological parameter is varied from the fiducial cosmology. The variation in the cosmological parameters includes an increment and decrement in the cosmological parameter in consideration. The step size for this change utilised in the Quijote suite of simulations are as follows: dΩ m = 0.010, dΩ b = 0.002, dh = 0.020, dn s = 0.020, and dσ 8 = 0.015. Furthermore, the suite provides 500 realisations for the three massive neutrino cosmologies, i.e. M ν = 0.1, 0.2 and 0.4 eV. The Quijote suite provides an additional 500 random realisations for the reference cosmology, in which the ICs were generated using the Zel'dovich approximation. This allowed us to compute the numerical derivatives with respect to massive neutrinos.
In this work, we used halo catalogue data from 22,000 Nbody simulations of the Quijote suite. These halos were identified using a friends-of-friends algorithm. We selected halos that have 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 n ∼ 0.92 × 10 −4 h 3 Mpc −3 for the reference simulations. For further details regarding the Quijote suite of simulations, we refer the reader to Villaescusa-Navarro et al. (2020).

Fisher-matrix formalism
In this work, we used the Fisher-matrix formalism to quantify the error estimates on the cosmological parameters. The Fisher information matrix can be defined as (e.g. Tegmark et al. 1997;Heavens 2009;Verde 2010) where L is the likelihood of the data given a model, and θ α and θ β are two of the model parameters. Under the assumption that the probability distribution of the data is a Gaussian (and hence has a Gaussian likelihood), we can write the Fisher information matrix as where S represents the data vector for the summary statistics we are interested in. In our case, it is the mean pairwise velocity and the mean relative velocities between pairs in a triplet. The precision matrix (i.e. the inverse covariance matrix) is given bŷ C −1 , and to obtain it, we computed the covariance matrix of the desired statistics directly from the simulations as follows: The total number of simulations used to compute the covariance matrix is denoted by N sims , and in our case N sims = 15, 000. An unbiased estimate of the covariance matrix is obtained using Eq. 3. However, its inversion leads to a biased estimate of the precision matrix. It can be corrected statistically by a multiplicative correction factor (Kaufmann 1967;Anderson 2003;Hartlap et al. 2007): where N bins is the number of bins in the data vector. It should be noted that by using Eq. 2, we assumed that the covariance matrix is independent of cosmology. However, the correction due to the cosmology dependence of the covariance matrix is expected to be negligible (Kodwani et al. 2019).
The derivatives required to construct the Fisher information matrix were numerically computed using the Quijote suite of simulations. In the case when the model parameter θ ≡ {Ω m , Ω b , h, n s , σ 8 }, we made use of the central difference approximation: As mentioned in Sect. 2.1, for each cosmological parameter considered the Quijote suite provides 500 realisations where only one parameter is varied, while the rest are fixed at the fiducial cosmological model. This enables the derivatives to be calculated numerically.
In the case of the neutrino mass, the fiducial value is 0.0 eV, and it cannot have negative values; hence, we obtained the partial derivative using Thus, we utilised two sets of massive neutrino simulations from Quijote, with M ν = 0.2 eV and M ν = 0.4 eV for the Fisher information matrix. These simulations had the initial conditions generated using Zel'dovich approximation. For the case of M ν = 0 eV for the above derivative, to be consistent, we hence made use of another 500 realisations of reference cosmology simulation in which the initial conditions were also similarly generated using the Zel'dovich approximation.

Mean relative velocity statistics
In the single stream fluid approximation, the mean matter pairwise velocity can be defined as where δ i ≡ δ(x i ) is the density contrast, and u i ≡ u(x i ) ≡ u(x i )/aH is the normalised peculiar velocity (where a is the scale factor). At leading order in standard perturbation theory, we have w 12 |r 12 p δ 1 u 2 − δ 2 u 1 . Thus, the mean pairwise velocity can be written as (e.g. Fisher 1995;Juszkiewicz et al. 1998 where j 1 (x) = sin(x)/x 2 − cos(x)/x, P(k) denotes the linear matter power spectrum, the subscript p implies that the averages are computed over all particle pairs with separation r 12 , and the symbol indicates that the mathematical expression was truncated to leading order. It should be noted that due to gravity, on average the particles in a pair approach each other, i.e.,w(r 12 ) < 0. The fidelity of this analytical prediction was tested against simulations (for both dark matter particles and halos), and we found that it reproduces the measurements well above pair separations of 40-50 h −1 Mpc (e.g. Reid & White 2011;Kuruvilla & Porciani 2018). Improvements have been made on the leading order prediction from standard perturbation theory by incorporating different flavours of perturbation theory, such as convolutional Lagrangian perturbation theory (CLPT, Carlson et al. 2013;Wang et al. 2014) and convolutional Lagrangian effective field theory (CLEFT, Vlah et al. 2016). Accurate perturbation theory prediction accounting for massive neutrinos was introduced in Aviles & Banerjee (2020) and in the case of modified gravity in Valogiannis et al. (2020). In parallel, advances in semi-analytical approaches to modelling mean pairwise velocity has been shown to achieve 10-20% precision at low redshift for fairly non-linear scales of 5 h −1 Mpc and above (Shirasaki et al. 2021). The relative velocity between pairs can be generalised to triplets with separations 123 = (r 12 , r 23 , r 31 ) and is valid for all tracers (e.g. dark matter particles, halos, and galaxies). This was first done for dark matter in Kuruvilla & Porciani (2020). In such a case, we consider two mean relative velocities: w 12 | 123 t and w 23 | 123 t (where the subscript t implies that the averages are computed over all particle triplets with separations 123 ). The relative velocity between pair 12 in a triplet, in the single stream fluid approximation, can be written as Article number, page 3 of 13 A&A proofs: manuscript no. main Contrary to the case of mean pairwise velocities, the mean relative velocity between a particle pair in a triplet is not purely radial, it also has a transverse component in the plane of the triangle defined by the particles. The transverse component is generated by the gravitational influence of the third particle on the pair. We can decompose the mean relative velocity into its radial and transverse components: w 12 | 123 t = w 12 ·r 12 | 123 tr12 + w 12 ·t| 123 tt = R 12 ( 123 )r 12 + T 12 ( 123 )t , wheret = (r 23 − cos χr 12 )/ sin χ and χ = arccos(r 12 ·r 23 ). Thus, we obtain the radial (matter) component as r 12 + r 23 cos χ r 2 12 + r 2 23 + 2r 12 r 23 cos χ . (11) Similarly, the mean radial relative velocity can be written for the pair 23 in 123 : R 23 ( 123 ) =w(r 23 ) − 1 2 w(r 12 ) cos χ −w(r 31 ) r 23 + r 12 cos χ r 2 12 + r 2 23 + 2r 12 r 23 cos χ . (12) For detailed derivations of these mean relative velocities between pairs in a triplet, we refer the reader to Kuruvilla & Porciani (2020), in which we also showcase how well the standard perturbation theory prescription at leading order works for an unbiased tracer (i.e. for dark matter species). Utilising a massive number of simulations, we quantified the information content of mean relative velocity between halo pairs in a triplet. To compare the measurement against the theoretical prediction, we employed a linear (scale-independent and mass-dependent) bias b(M h ), and it can be written as and similarly where the superscript 'h' represents the mean statistics for the halos, and in our case M h implies all halos with mass greater than 5 × 10 13 h −1 M . To quantify the information content using the Fisher-matrix formalism, we measured the three-point mean relative velocity statistics for all triangular configurations with separation lengths r 12 ≥ r 23 ≥ r 31 ≥ 20 h −1 Mpc and separation lengths less than 120 h −1 Mpc. All separations use a bin width of 5 h −1 Mpc and thus have 1168 total triangular configurations. We directly computed R h i j from the simulations by measuring the average of the radial components of the relative velocities, i.e. (u j − u i ) ·r i j , and binning them to the corresponding triangles. For each halo catalogue of the Quijote simulation, this corresponds to computing the relative velocities for about a total number of 2 billion triangles at the scales we considered. In Fig. 1, we show the theoretical prediction alongside the direct measurements from the simulation of these relative velocities between halo pairs in a triplet. However, to plot  the comparison with the theoretical prediction, we show the direct measurements of R h 12 and R h 23 for triangle configurations with r 12 ≥ r 23 ≥ r 31 ≥ 50 h −1 Mpc and separation lengths of less than 120 h −1 Mpc. We restricted the plotting to these triangular configurations for two reasons: (i) we expected the analytical prediction to work well at these quasi-linear scales and above, and (ii) we reduced the number of total triangular configurations (from 1168 to 552 configurations) considered in the plot to highlight the configurations that are most interesting for comparison purposes. As stated before for all separations, we used 5 h −1 Mpc -wide bins. Thus, the smallest triangle (triangle configuration '0') considered in Fig. 1 corresponds to r 12 , r 23 , r 31 ∈ {(50, 55), (50, 55), (50, 55)} h −1 Mpc, while the largest-scale triangle (triangle configuration '552') has r 12 , r 23 , r 31 ∈ {(115, 120), (115, 120), (115, 120)} h −1 Mpc. In the top and bottom rows, we plot the relative velocity between pairs 12 and 23 in 123, respectively. The bias factor as given in Eq. 13 was determined by fitting the prediction using leastsquares method, and we obtain b(M h > 5 × 10 13 ) = 1.856. The direct measurement of relative velocities from the 15,000 refer-ence halo catalogues is denoted by the solid (blue) line, while the theoretical prediction using perturbation theory at leading order multiplied by the bias factor is given by a dashed (orange) line. The negative values of R 12 and R 23 imply that the pairs in the triplet are infalling towards each other, on average, due to the gravitational force. The rate of infall increases as the three separation lengths of the triangle decrease. We focused on a subset of configurations with the vertical blue shaded region to highlight these effects, where we fixed one leg of the triangle (100 < r 12 < 105 h −1 Mpc) and let r 23 and r 31 vary. It clearly shows how, contrary to the mean pairwise velocity, the addition of a third halo (or any other tracer of the matter distribution) imprints quite distinct features onto the velocity statistics. The leftmost data point in the vertical shaded region has its second and third legs of the triangle fixed as r 23 = r 31 ∈ (50, 55) h −1 Mpc. Thus, as the halo '3' is closer to the pair 12, the mean velocity R 12 experiences a greater infall speed. On the other hand, in the case of R 23 for the same triangular configuration, the effect is different as the presence of halo '1' causes a reduction in the net infall velocity. The relative velocity between pairs R 12 and R 23 will be equal to each other when r 12 = r 23 , irrespective of the length of the third side. When comparing the theoretical predictions, we see that it is, overall, accurate within 4-5%. As can be seen from the residuals of both R h 12 and R h 23 , the analytical predictions seem somewhat sensitive to the triangle shapes. This is mainly due to the fact that even for the mean halo pairwise velocity at these scales, the linear theory prediction is not perfect, and it is accurate at roughly below 2-4% for certain scales of pair separation between 50 and 120 h −1 Mpc. This has been shown in the case of mean halo pairwise velocity in Fig. 6 of Uhlemann et al. (2015), wherein the analytical linear theory prescription over-predicts by 2-4% at around 90-100 h −1 Mpc, and similarly under-predicts by 3-5% at around 50 h −1 Mpc. As the leading order prediction for R h 12 and R h 23 is built from the sum of different pairwise components, these small inaccuracies propagate, and we see it in the residuals. If we focus on the triangular configuration '165', we find that both R h 12 and R h 23 over-predict the measurement by 4.6%, and this configuration corresponds to the equilateral triangle with a separation length between 90 and 95 h −1 Mpc. However, the overall degree of agreement of the analytical prediction with the direct measurement from simulation for all these configurations is quite encouraging, even with a simple biasing scheme. Hence, we believe if we were to restrict our present analysis to triangular configurations above 50 h −1 Mpc, and assuming a massless neutrino cosmology, the leading order prediction would be sufficient for an analysis such as the one we undertook here. To further push to smaller non-linear scales, one would need to improve the analytical predictions for R h i j by incorporating different types of perturbation theory like CLPT or CLEFT, as has been done in the case of mean pairwise velocities. This provides us with the motivation to utilise direct measurements from the simulations to quantify the information content in these velocity statistics as we can extend the Fisher analysis to non-linear scales where the leading order theoretical predictions would otherwise not be accurate enough.

Impact of cosmology
In this section, we show how different cosmological parameters affect the relative velocities between halo pairs in a triplet. We focus on the shape dependence of lar configurations in which the first leg of the triangle is fixed to r 12 ∈ (115, 120) h −1 Mpc. In each panel, the top right part where r 31 /r 12 = r 23 /r 12 = 1 corresponds to equilateral configuration. Bottom bins (i.e. r 31 /r 12 = r 23 /r 12 ∼ 0.5) contain the folded configurations. Since our analysis has a minimum separation length of 20 h −1 Mpc, we do not explore squeezed configurations (r 12 ∼ r 23 r 31 ). The left, middle, and right panels show the effect when the summed mass of neutrinos is 0.1, 0.2, and 0.4 eV, respectively. For the fiducial cosmology simulations, the 500 realisations of which were generated using the Zel'dovich approximation (labelled 'fidZA') were utilised. However, it will not make any difference in the ratio for these triangular configurations whether the simulations generated using second-order Lagrangian perturbation theory or Zel'dovich approximation are used, as we verified that both are numerically consistent with each other at ∼ 0.05% at these scales. In the top panel, we see that by increasing the neutrino mass, the infall rate of R h 12 increases. The maximal difference as before is for the equilateral triangle, and minimal when r 31 r 23 < r 12 . In the case of R h 23 , as seen from the bottom panel, the trend is similar when r 31 /r 12 > 0.4. However, for triangles that are close to squeezed configurations, the trend reverses, wherein on increasing the neutrino mass the infall rate of R h 23 decreases. For the particular triangular configurations as shown in relative velocity statistics. We see that increasing Ω m decreases R h 12 and vice-versa. This overall behaviour is seen with variations in the Hubble parameter also, albeit with a different amplitude. However, for Ω b and σ 8 , the behaviour is reversed, and the amplitude of the change of R h 12 with respect to its fiducial counterpart is within 6% for these configurations. For the same triangular configurations, at r 31 /r 12 0.4, the variation in cosmological parameters yields a similar response for R h 23 as R h 12 , albeit with a different amplitude. However, this changes in the case of r 31 /r 12 < 0.4, where in the behaviour reverses. The reverse is seen least in the case of variation in σ 8 . R h 12 and R h 23 show the maximal difference in response to cosmological parameters in the case of equilateral triangles and at r 31 /r 12 < 0.4, respectively. We refer the reader to Appendix B to see the figures showing how these different cosmological parameters affect the mean three-point relative velocity statistics.
When increasing the neutrino mass (as in the case of M ν = 0.4 eV), there is a more pronounced shape dependence compared to σ 8 . This should facilitate breaking the M ν -σ 8 degeneracy in two-point clustering information. To further understand this, we take a look at the derivatives for R h 12 (left panel) and R h 23 (right panel) in Fig. 3 for all triangular configurations with r 12 ≥ r 23 ≥ r 31 ≥ 20 h −1 Mpc and a maximum separation of 120 h −1 Mpc. Each row corresponds to the derivative with respect to a cosmological parameter (labelled in each row on the right panel). For the very large-scale triangles, the shape for ∂R h 23 /∂σ 8 and ∂R h 23 /∂M ν looks similar. However, for the small-scale triangles, this is not the case. This points to the possibility that R h 23 + R h 23 can be used to break the M ν -σ 8 degeneracy.

Cosmological parameters
As can be seen from Eq. 2, one of the quantities entering the Fisher-matrix formalism is the inverse of the covariance matrix. It can be determined theoretically, in principle. For the mean pairwise velocities, for example, this was done in Bhattacharya & Kosowsky (2008) and Mueller et al. (2015a). As the mean relative velocity between pairs in a triplet was introduced recently, no analytical prediction for its covariance matrix has been developed. We aim to tackle this in a future work. However in this work, we took the alternative route and measured the covariance matrix of the necessary summary statistics directly from the simulations using 15,000 realisations. Thus, in the case of R h 12 and R h 23 , this allowed us to compute the covariance matrix reliably. We show the cross-correlation coefficient (ρ i j = C i j / C ii C j j ) for both the mean halo pairwise velocity (left panel) and the mean relative velocity between halo pairs in a triplet (right panel) in Fig. 4. We see that in both panels, it consists of prominent nondiagonal contributions. We especially notice that R h 12 and R h 23 are highly correlated; this is not surprising, as in the leading order Cross -correlation coefficient prediction we can see that they are a linear combination of different pieces of pairwise velocities of the three legs of the triangle (as shown in Eqs. 11 and 12). The other quantities that enter Eq. 2 are the derivatives of the model predictions with respect to the model parameters. The most common method to compute these derivatives is to use the theoretical predictions, such as those from the leading order in standard perturbation theory (Eq. 13). We noted that these predictions are accurate at a few percent for triangular configurations above 50 h −1 Mpc. At non-linear scales, the fidelity of the analytical predictions decreases. Hence, to circumvent this effect, we also used the derivatives measured directly from the simulations. We used 500 realisations, which were evaluated at thirteen different cosmologies, and each cosmological model with a small variation as already mentioned in Sect. 2.1.
Combining both the necessary derivatives and the covariance matrix for the Fisher-matrix formalism, we present the constraints on the cosmological parameters in Fig. 5. The 68.3 and 95.4 percent credible regions are given by dark and light shaded contour regions, respectively, for all possible pairs of model parameters obtained after marginalising over all the remaining parameters. The blue and orange ellipses show the constraints from the mean halo pairwise velocity and the mean relative velocity between halo pairs in a triplet, respectively. We can immediately note the benefits of constraining the cosmological parameters using the three-point relative velocities. We used the nuisance parameter b as a scaling factor for w h 12 p , R h 12 , and R h 23 . For r min = 20 h −1 Mpc, we find that R h 12 + R h 23 has a 1σ constraints for {Ω m , Ω b , h, n s , σ 8 , M ν } = {0.0091, 0.0024, 0.0226, 0.0221, 0.0156, 0.0655 eV}, which presents an improvement over constraints from w h 12 p by a factor of {8.91, 11.81, 15.46, 20.97, 10.93, 13.40}. The three-point mean halo relative velocity statistics thus offer a substantial improvement on the model parameter constraints over the mean halo pairwise velocity. To understand the role that the non-linear scales play in this improvement, we constrained the model parameters using Fisher matrix formalism where we restricted the triangular configurations at different r min (keeping r max fixed). We show the results in Table 1. For all the parameters considered, we see how the constraints improve when going to smaller scales, as expected. The improvement obtained from R h 12 + R h 23 over those obtained from mean pairwise velocity increases as the minimum separation scale considered decreases. For the mean pairwise velocity, one can see that the spectral index n s is highly degenerate with Ω m , Ω b , and h. This degeneracy is also observed in the case of constraints from halo ) and galaxy (Hahn & Villaescusa-Navarro 2021) power spectrum. However, by considering the mean relative velocity between pairs in a triplet, the n s -Ω b degeneracy is broken. This is also supported by the fact that the impact of n s on R h 12 and R h 23 is different from those of Ω b as supported by Fig. 3  (also Figs. B.1 and B.2).
In this work, we made use of Fisher matrix formalism to study the informational content, where the derivatives and the covariance matrix was directly computed numerically from the Quijote suite of simulations. The numerical stability of our results stemming from the number of simulations used to compute the covariance matrix and derivatives is discussed in depth in Appendix A. We see that the constraints on the cosmological parameters vary less well below 1% up on using 10,000 or more simulations for computing the covariance matrix. While in the case of the derivatives, the cosmological parameters are within around 5% using 450 realisations when compared to the estimates using the full 500 realisations. The constraint on the nuisance parameter, b, varies less than 2.5% and 5% using 400 and 450 realisations, respectively. Hence, we conclude that our results are stable and devoid of any numerical inconsistencies.
So far, we considered the information content from combining R h 12 and R h 23 , but it is also worth exploring the constraining power of the two quantities when considered separately. We followed the same procedure as before to calculate the Fisher information matrix, and we computed the necessary derivatives and the covariance matrix directly from the simulations. We considered all triangular configurations with a minimum and a maximum separation of 20 h −1 Mpc and 120 h −1 Mpc. In one case, we considered the constraints from R h 12 and R h 23 separately, and in the second case we consider the joint R h 12 + R h 23 as has been done in the main text. We showcase the results in Fig. 6. The blue contours show the joint credible regions for all possible pairs of parameters using R h 12 alone, while the orange contour shows that from R h 12 +R h 23 . In the figure, we have not shown the joint credible regions from R h 23 because it is the same as R h 12 . We can see that on all parameters R h 12 + R h 23 is able to obtain tighter constraints. Quantitatively, they are 1.34 -1.46 times better than those considering only R h 12 or R h 23 . In our main analysis, we did not incorporate any observation effects except for the bias term, and thus our results are provided with the ideal constraining power from these relative velocity statistics. In the case of kSZ experiments, one would have to additionally consider the optical depth as an additional free parameter, which depends on the halo mass (Battaglia 2016;Flender et al. 2017). We tried to (naively) mimic this halo mass dependence of optical depth by computing the derivative by changing the minimum halo mass considered in addition to the bias term and included it in the Fisher matrix analysis. We find that the constraining power on all six cosmological parameters remains unchanged at less than 0.1%. However, the 1-σ constraint on b degrades from 0.019 to 0.022 (decrease by 12.7%) after adding M h,min as an additional parameter in the Fisher analysis. Accounting for additional observational effects related to the actual measurement of the kSZ signal such as beam convolution, contamination by other astrophysical signals, foreground cleaning, etc. would require an actual analysis based on kSZ-like mocks. This represents a quite different analysis that goes beyond the goals of the present work.
We can compare our constraints with the ones obtained using multipoles of the power spectrum and bispectrum to see how well the three-point mean relative velocity statistics perform. The constraints from Hahn et al. (2020) provide us with Notes. Different columns show the constraints obtained by changing the minimum separation length considered, while fixing the maximum separation to 120 h −1 Mpc. The various rows display the results for different parameters, and each of these rows is further divided into three rows to focus on the constraint coming from w 12 |r 12 p , R h 12 + R h 23 , and the factor of improvement of the constraints from R h 12 + R h 23 over w 12 |r 12 p denoted by F.I., respectively. an avenue to obtain a fair comparison of our summary statistics based on relative velocity with the clustering statistics, as both studies utilise the same halo catalogues. A caveat being the number density is slightly lower in our case as we chose a higher minimum mass of halos (M h > 5 × 10 13 h −1 M , while in Hahn et al. (2020) they considered M h > 3.2 × 10 13 h −1 M ). Using the monopole and quadrupole of the redshift-space power spectrum for k max = 0.5 h Mpc −1 , Hahn et al. (2020) were able to obtain 1-σ constraint of 0.294 eV for the summed neutrino mass. Thus, our constraints on M ν (i.e. 0.065 eV) from R h 12 + R h 23 (for r min = 20 h −1 Mpc) give an improvement factor of 4.5 over those obtained from the multipoles of the power spectrum. This improvement also crosses over to other cosmological parameters that both studies considered, with an improvement factor of (2.6, 4.8, 4.9, 5.7, 2.3) for (Ω m , Ω b , h, n s , σ 8 ), respectively. Hahn et al. (2020) also provided the cosmological constraints from using the redshift-space bispectrum monopole (the threepoint clustering statistics in Fourier space). Our constraint on M ν from the three-point relative velocity statistics is slightly less than the 1-σ 8 constraints from the bispectrum monopole (0.054 eV). However, for the other parameters such as (Ω m , Ω b , h, n s ), we notice a small increase in the constraining power of threepoint mean relative velocity when compared to those from bispectrum monopole with a factor of (1.2, 1.7, 1.7, 1.5), respectively. Albeit σ 8 produces a slight decrease in the constraining power with a factor of 0.9. This shows the potential of threepoint mean relative velocity statistics to be competitive with respect to clustering information in the future. Kuruvilla & Porciani (2020) introduced the mean relative velocity between pairs in a triplet, and derived the analytical prediction for it at leading order using standard perturbation theory. In this work, we quantified the information content in this three-point relative velocity statistics. Furthermore, we showed the fidelity of analytical the prediction for halos. We accounted for the biasing by introducing a simple (scale-independent and mass-dependent) bias term. For triangular configurations with r min = 50 h −1 Mpc and r max = 120 h −1 Mpc, we find that these predictions are accurate within a few percent.

Conclusion
We checked the effect of cosmological parameters on the mean relative velocities between halo pairs in a triplet. We find distinctive shape dependence on R h 12 and R h 23 as a result of varying cosmological parameters for the triangular configurations we considered. With regard to the summed mass of neutrinos, the shape dependence becomes more prominent as M ν increases.
To quantify the information content in the two-and threepoint mean relative velocity statistics, we used the Fisher-matrix formalism, where the necessary derivatives and the covariance