Issue 
A&A
Volume 653, September 2021



Article Number  A130  
Number of page(s)  12  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/202140552  
Published online  23 September 2021 
Information content in mean pairwise velocity and mean relative velocity between pairs in a triplet
Université ParisSaclay, CNRS, Institut d’Astrophysique Spatiale, 91405 Orsay, France
email: joseph.kuruvilla@universiteparissaclay.fr
Received:
12
February
2021
Accepted:
11
July
2021
Velocity fields provide 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 threepoint mean relative velocity (i.e. the mean relative velocities between pairs in a triplet). Using halo catalogs from the Quijote suite of Nbody 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 ≥ 50 h^{−1} Mpc. Furthermore, we present the mean relative velocity between pairs in a triplet as a novel probe of neutrino mass estimation. We explored the full cosmological information content of the halo mean pairwise velocities and the mean relative velocities between halo pairs in a triplet. We did this through the Fishermatrix formalism using 22 000 simulations from the Quijote suite and by 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σ neutrino mass (M_{ν}) constraint of 0.065 eV, which is roughly 13 times better than the mean pairwise velocity constraint (0.877 eV). This information gain is not limited to neutrino mass, but it extends to other cosmological parameters: Ω_{m}, Ω_{b}, h, n_{s}, and σ_{8}, achieving an information gain of 8.9, 11.8, 15.5, 20.9, and 10.9 times, respectively. These results illustrate the possibility of exploiting the mean threepoint relative velocities to constrain the cosmological parameters accurately from future cosmic microwave background experiments and peculiar velocity surveys.
Key words: largescale structure of Universe / cosmology: theory
© J. Kuruvilla and N. Aghanim 2021
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
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 VI 2020), galaxy redshift surveys (e.g. eBOSS Collaboration 2021), 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 CMBS4^{3} (Abazajian et al. 2016).
Cosmological information from galaxy redshift surveys is often extracted using the twopoint correlation function in the configuration space or its Fourier analogue, the power spectrum. However, as the complex galaxy distribution represents a nonGaussian field, higher order statistics such as the bispectrum and the threepoint correlation function contain additional information. Thus, by considering both twopoint and threepoint clustering information, future redshift surveys will be able to obtain a substantial improvement with regard to the cosmological parameter constraints (e.g. Yankelevich & Porciani 2019; Chudaykin & Ivanov 2019; Gualdi & Verde 2020; 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 freestreaming 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; VillaescusaNavarro et al. 2018; GarcíaFarieta 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 Lband Legacy Allsky 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 Cosmicflows3 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ándezMonteagudo et al. (2015), Planck Collaboration Int. XXXVII (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, 2008; Kosowsky & 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ándezMonteagudo 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 (ChavesMontero et al. 2021). The forthcoming CMB surveys including SO, CMBS4, and CMBHD (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 nonlinear scales was studied in Kuruvilla et al. (2020) 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 quasilinear 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 higherorder relative velocity statistics provide on the cosmological parameters, in particular with regard to summed neutrino mass. The threepoint 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 Fishermatrix 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 (VillaescusaNavarro et al. 2020). Thus, this work forms a parallel line of study to Hahn et al. (2020) and Hahn & VillaescusaNavarro (2021), which quantified the information content in the redshiftspace bispectrum of halos and galaxies, respectively, using the Fishermatrix formalism; moreover, they computed the necessary quantities directly from the Quijote simulations. It was shown that the redshiftspace 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 Fishermatrix 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.
2. Data and analysis
2.1. Quijote simulation suite
The Quijote^{6} (VillaescusaNavarro et al. 2020) simulation suite contains 44 100 Nbody simulations spanning more than a few thousand cosmological models, run using the treePM code GADGET3 (Springel 2005). The initial conditions (ICs) of the simulation, at redshift z = 127, were generated using the secondorder 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 VI (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 (presentday 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 friendsoffriends 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 Mpc^{−3} for the reference simulations. For further details regarding the Quijote suite of simulations, we refer the reader to VillaescusaNavarro et al. (2020).
2.2. Fishermatrix formalism
In this work, we used the Fishermatrix 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 ℒ 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 by , and to obtain it, we computed the covariance matrix of the desired statistics directly from the simulations as follows:
where . 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.
3. 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 v_{i} ≡ v(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}v_{2}⟩−⟨δ_{2}v_{1}⟩. Thus, the mean pairwise velocity can be written as (e.g. Fisher 1995; Juszkiewicz et al. 1998; Reid & White 2011)
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., . 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 semianalytical approaches to modelling mean pairwise velocity has been shown to achieve 10−20% precision at low redshift for fairly nonlinear 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
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:
where and . Thus, we obtain the radial (matter) component as
Similarly, the mean radial relative velocity can be written for the pair 23 in △_{123}:
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 (scaleindependent and massdependent) bias , and it can be written as
and similarly
where the superscript ‘h’ represents the mean statistics for the halos, and in our case implies all halos with mass greater than 5 × 10^{13} h^{−1} M_{⊙}. To quantify the information content using the Fishermatrix formalism, we measured the threepoint 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 from the simulations by measuring the average of the radial components of the relative velocities, i.e. , 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 and 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 quasilinear 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 largestscale 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 reference 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 and , 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 overpredicts by 2−4% at around 90−100 h^{−1} Mpc, and similarly underpredicts by 3−5% at around 50 h^{−1} Mpc. As the leading order prediction for and 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 and overpredict 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 nonlinear scales, one would need to improve the analytical predictions for 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 nonlinear scales where the leading order theoretical predictions would otherwise not be accurate enough.
Fig. 1. Top panel: measurement (solid line) from the halo catalogues and the analytical prediction (dashed line) using linear perturbation theory as given in Eq. (13). Bottom panel: same for . The plot shows all triangular configurations, in which the minimum and maximum separations for a triangular side are 50 and 120 h^{−1} Mpc, respectively. Light blue vertical band corresponds to triangular configurations where one leg of the triangle is fixed, r_{12} ∈ (100, 105) h^{−1} Mpc. 
4. Results
4.1. 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 (top panel) and (bottom panel) on M_{ν} in Fig. 2. The figure focuses on the triangular 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 secondorder 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 increases. The maximal difference as before is for the equilateral triangle, and minimal when r_{31} ≪ r_{23} < r_{12}. In the case of , 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 decreases.
Fig. 2. Top: effect of M_{ν} on . The length of the fixed leg of the triangle corresponds to r_{12} ∈ (115, 120) h^{−1} Mpc. Each panel corresponds to a different M_{ν}, and it is indicated in the bottom right of each panel (in units of eV). The initial conditions of the fiducial cosmology simulations used in this figure were generated using the Zel’dovich approximation. Bottom: same as top row, but for . 
For the particular triangular configurations as shown in Fig. 2, we checked how other cosmological parameters affect the relative velocity statistics. We see that increasing Ω_{m} decreases and viceversa. 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 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 as , 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}. and 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 threepoint 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 twopoint clustering information. To further understand this, we take a look at the derivatives for (left panel) and (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 largescale triangles, the shape for and looks similar. However, for the smallscale triangles, this is not the case. This points to the possibility that can be used to break the M_{ν} − σ_{8} degeneracy.
Fig. 3. Derivatives of (left), and (right) with respect to the different cosmological parameters for all triangular configurations, in which the minimum and maximum separation for a triangular side is 20 and 120 h^{−1} Mpc, respectively. The cosmological parameter under consideration is denoted on the right side of each row. 
4.2. Cosmological parameters
As can be seen from Eq. (2), one of the quantities entering the Fishermatrix 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 and , this allowed us to compute the covariance matrix reliably. We show the crosscorrelation coefficient () 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 and are highly correlated; this is not surprising, as in the leading order 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)).
Fig. 4. Left panel shows correlation matrix in the case of halo pairwise velocities, while the right panel shows it for the relative velocities between halo pairs in a triplet computed using the 15 000 realisations of the fiducial run of the Quijote suite of simulations. For the right panel, we used all triangular configurations with 20 ≤ r_{31} ≤ r_{23} ≤ r_{12} ≤ 120 h^{−1} Mpc. 
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 nonlinear 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 Fishermatrix 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 threepoint relative velocities. We used the nuisance parameter b as a scaling factor for , , and . For r_{min} = 20 h^{−1} Mpc, we find that 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 by a factor of {8.91, 11.81, 15.46, 20.97, 10.93, 13.40}. The threepoint 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 nonlinear 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 over those obtained from mean pairwise velocity increases as the minimum separation scale considered decreases.
Fig. 5. Joint 68.3% (dark shaded contour) and 95.4% (light shaded contour) credible region for all the pairs of model parameters (six cosmological and one nuisance parameter) at z = 0. Colourcoding is indicated in the legend. 
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 (Hahn et al. 2020) and galaxy (Hahn & VillaescusaNavarro 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 and 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 to compute 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 and , 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 and separately, and in the second case we consider the joint 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 alone, while the orange contour shows that from . In the figure, we have not shown the joint credible regions from because it is the same as . We can see that on all parameters is able to obtain tighter constraints. Quantitatively, they are 1.34−1.46 times better than those considering only or .
Fig. 6. Joint 68.3% (dark shaded contour) and 95.4% (light shaded contour) credible region for all the pairs of model parameters (six cosmological and one nuisance parameter) at z = 0 for (blue) and (orange), respectively. 
Marginalised 1σ errors for the cosmological parameters and the nuisance parameter (b).
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 kSZlike 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 threepoint mean relative velocity statistics perform. The constraints from Hahn et al. (2020) provide us with 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 redshiftspace power spectrum for k_{max} = 0.5 h^{−1} Mpc, 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 (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 redshiftspace bispectrum monopole (the threepoint clustering statistics in Fourier space). Our constraint on M_{ν} from the threepoint 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.
5. Conclusion
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 threepoint relative velocity statistics. Furthermore, we showed the fidelity of analytical the prediction for halos. We accounted for the biasing by introducing a simple (scaleindependent and massdependent) 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.
We checked the effect of cosmological parameters on the mean relative velocities between halo pairs in a triplet. We find distinctive shape dependence on and 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 Fishermatrix formalism, where the necessary derivatives and the covariance matrices were directly measured from the Quijote suite of simulations. For the covariance matrices, we utilised 15 000 realisations of the reference cosmology. The derivatives for the threepoint mean relative velocity were also computed directly from the simulations. Combining these, we obtained the Fisher constraints by utilising the mean halo pairwise velocities and . In the case of pairwise velocities, we used pair separations from 20 h^{−1} Mpc to 120 h^{−1} Mpc, with a bin width of 5 h^{−1} Mpc. In the case of and , we considered all triangular configurations with a minimum triangle side length of 20 h^{−1} Mpc and a maximum of 120 h^{−1} Mpc. This amounted to using 1168 total triangular configurations. was able to achieve a 1σ constraint of 0.0655 eV for M_{ν}, which is factor of 13.40 improvement over constraints from the mean pairwise velocities. The current lower limit of the sum of neutrino masses, M_{ν} = ∑m_{ν} ≳ 0.06 eV, comes from the neutrino oscillation experiments (e.g. Forero et al. 2014; GonzalezGarcia et al. 2016; Capozzi et al. 2017; de Salas et al. 2018). Thus, the threepoint mean relative velocity statistics could help us to pinpoint the summed mass of neutrinos in the future. Similarly, for the other cosmological parameters, the threepoint mean relative velocities made it possible to achieve similar information gains over pairwise velocity and can be seen quantitatively in Table 1. This information gain is possible as and , being a threepoint relative velocity statistics, make more shape dependence variation possible, with different cosmological parameters compared to the mean pairwise velocity. We did not study the effect of binning on the constraints from the mean relative velocity statistics, which uses a rather broad bin width and hence could lead to information loss, especially in the case of mean pairwise velocities. Though in this work our goal was to have a fair comparison of the information content of both the twopoint and threepoint mean relative velocities, and hence we kept the same bin width.
In summary, we quantitatively show that the mean radial relative velocity between pairs in a triplet as a cosmological observable at redshift zero could lead to a significant information gain in comparison to the mean radial pairwise velocity. The constraining power from is comparable to those obtained from the halo bispectrum measurement (Hahn et al. 2020) as detailed in Sect. 4.2. This represents the first step in understanding the viability of utilising the mean relative velocity between pairs in a triplet as a cosmological observable from either future peculiar velocity surveys or kinetic Sunyaev–Zel’dovich experiments.
Acknowledgments
We thank the referee for their constructive comments, which helped in improving the manuscript. JK likes to thank Francisco VillaescusaNavarro for all the valuable and helpful discussions regarding the Quijote simulation suite. We thank Julien Grain, Tony Bonnaire and Hideki Tanimura for carefully reading the manuscript, and providing comments. We acknowledge 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, JCAP, 2019, 056 [CrossRef] [Google Scholar]
 Agarwal, N., Desjacques, V., Jeong, D., & Schmidt, F. 2021, JCAP, 2021, 021 [Google Scholar]
 eBOSS Collaboration (Alam, S.), et al. 2021, Phys. Rev. D, 103, 083533 [Google Scholar]
 Anderson, T. 2003, An Introduction to Multivariate Statistical Analysis, Wiley Series in Probability and Statistics (Wiley) [Google Scholar]
 Aviles, A., & Banerjee, A. 2020, JCAP, 2020, 034 [Google Scholar]
 Battaglia, N. 2016, JCAP, 2016, 058 [Google Scholar]
 Behnel, S., Bradshaw, R., Citro, C., et al. 2011, Comput. Sci. Eng., 13, 13 [CrossRef] [Google Scholar]
 Bhattacharya, S., & Kosowsky, A. 2007, ApJ, 659, L83 [NASA ADS] [CrossRef] [Google Scholar]
 Bhattacharya, S., & Kosowsky, A. 2008, Phys. Rev. D, 77, 083004 [NASA ADS] [CrossRef] [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 [CrossRef] [Google Scholar]
 Carlson, J., Reid, B., & White, M. 2013, MNRAS, 429, 1674 [NASA ADS] [CrossRef] [Google Scholar]
 Castorina, E., Carbone, C., Bel, J., Sefusatti, E., & Dolag, K. 2015, JCAP, 2015, 043 [NASA ADS] [CrossRef] [Google Scholar]
 ChavesMontero, J., HernándezMonteagudo, C., Angulo, R. E., & Emberson, J. D. 2021, MNRAS, 503, 1798 [NASA ADS] [CrossRef] [Google Scholar]
 Chudaykin, A., & Ivanov, M. M. 2019, JCAP, 2019, 034 [CrossRef] [Google Scholar]
 da Cunha, E., Hopkins, A. M., Colless, M., et al. 2017, PASA, 34, e047 [NASA ADS] [CrossRef] [Google Scholar]
 De Bernardis, F., Aiola, S., Vavagiakis, E. M., et al. 2017, JCAP, 3, 008 [Google Scholar]
 de Salas, P. F., Forero, D. V., Ternes, C. A., Tortola, M., & Valle, J. W. F. 2018, Phys. Lett. B, 782, 633 [CrossRef] [Google Scholar]
 Djorgovski, S., & Davis, M. 1987, ApJ, 313, 59 [Google Scholar]
 Dressler, A., LyndenBell, D., Burstein, D., et al. 1987, ApJ, 313, 42 [Google Scholar]
 Dupuy, A., Courtois, H. M., & Kubik, B. 2019, MNRAS, 486, 440 [NASA ADS] [CrossRef] [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 [CrossRef] [Google Scholar]
 GonzalezGarcia, M. C., Maltoni, M., & Schwetz, T. 2016, Nucl. Phys. B, 908, 199 [CrossRef] [Google Scholar]
 Gualdi, D., & Verde, L. 2020, JCAP, 2020, 041 [Google Scholar]
 Hahn, C., & VillaescusaNavarro, F. 2021, JCAP, 2021, 029 [Google Scholar]
 Hahn, C., VillaescusaNavarro, F., Castorina, E., & Scoccimarro, R. 2020, JCAP, 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, 357 [CrossRef] [PubMed] [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]
 Heymans, C., Tröster, T., Asgari, M., et al. 2021, A&A, 646, A140 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Howlett, C., StaveleySmith, L., & Blake, C. 2017, MNRAS, 464, 2517 [CrossRef] [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]
 Kamalinejad, F., & Slepian, Z. 2020, ArXiv eprints [arXiv:2011.00899] [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]
 Koda, J., Blake, C., Davis, T., et al. 2014, MNRAS, 445, 4267 [NASA ADS] [CrossRef] [Google Scholar]
 Kodwani, D., Alonso, D., & Ferreira, P. 2019, Open J. Astrophys., 2, 3 [CrossRef] [Google Scholar]
 Koribalski, B. S., StaveleySmith, L., Westmeier, T., et al. 2020, Ap&SS, 365, 118 [CrossRef] [Google Scholar]
 Kosowsky, A., & Bhattacharya, S. 2009, Phys. Rev. D, 80, 062003 [CrossRef] [Google Scholar]
 Kuruvilla, J., & Porciani, C. 2018, MNRAS, 479, 2256 [CrossRef] [Google Scholar]
 Kuruvilla, J., & Porciani, C. 2020, JCAP, 2020, 043 [CrossRef] [Google Scholar]
 Kuruvilla, J., Aghanim, N., & McCarthy, I. G. 2020, A&A, 644, A170 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Li, Y.C., Ma, Y.Z., Remazeilles, M., & Moodley, K. 2018, Phys. Rev. D, 97, 023514 [NASA ADS] [CrossRef] [Google Scholar]
 McCarthy, I. G., Bird, S., Schaye, J., et al. 2018, MNRAS, 476, 2999 [NASA ADS] [CrossRef] [Google Scholar]
 Mueller, E.M., de Bernardis, F., Bean, R., & Niemack, M. D. 2015a, ApJ, 808, 47 [CrossRef] [Google Scholar]
 Mueller, E.M., de Bernardis, F., Bean, R., & Niemack, M. D. 2015b, Phys. Rev. D, 92, 063501 [NASA ADS] [CrossRef] [Google Scholar]
 Nguyen, N.M., Jasche, J., Lavaux, G., & Schmidt, F. 2020, JCAP, 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]
 Sehgal, N., Aiola, S., Akrami, Y., et al. 2019a, BAAS, 51, 6 [NASA ADS] [Google Scholar]
 Sehgal, N., Nguyen, H. N., Meyers, J., et al. 2019b, BAAS, 51, 43 [Google Scholar]
 Shirasaki, M., Huff, E. M., Markovic, K., & Rhodes, J. D. 2021, ApJ, 907, 38 [NASA ADS] [CrossRef] [Google Scholar]
 Soergel, B., Flender, S., Story, K. T., et al. 2016, MNRAS, 461, 3172 [Google Scholar]
 Springel, V. 2005, MNRAS, 364, 1105 [Google Scholar]
 Sunyaev, R. A., & Zeldovich, Y. B. 1972, Comments Astrophys. Space Phys., 4, 173 [NASA ADS] [EDP Sciences] [Google Scholar]
 Sunyaev, R. A., & Zeldovich, I. B. 1980, MNRAS, 190, 413 [NASA ADS] [CrossRef] [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]
 Tully, R. B., & Fisher, J. R. 1977, A&A, 500, 105 [Google Scholar]
 Tully, R. B., Courtois, H. M., & Sorce, J. G. 2016, AJ, 152, 50 [NASA ADS] [CrossRef] [Google Scholar]
 Uhlemann, C., Kopp, M., & Haugg, T. 2015, Phys. Rev. D, 92, 063004 [NASA ADS] [CrossRef] [Google Scholar]
 Valogiannis, G., Bean, R., & Aviles, A. 2020, JCAP, 2020, 055 [Google Scholar]
 Verde, L. 2010, Statistical Methods in Cosmology (Berlin/Heidelberg: SpringerVerlag), 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]
 Vlah, Z., Castorina, E., & White, M. 2016, JCAP, 2016, 007 [Google Scholar]
 Wang, L., Reid, B., & White, M. 2014, MNRAS, 437, 588 [NASA ADS] [CrossRef] [Google Scholar]
 Wong, Y. Y. Y. 2008, JCAP, 2008, 035 [NASA ADS] [CrossRef] [Google Scholar]
 Yankelevich, V., & Porciani, C. 2019, MNRAS, 483, 2078 [CrossRef] [Google Scholar]
Appendix A: Numerical stability
The Fisher matrix (Eq. 2) calculation depends on two elements: (i) the covariance matrix, and (ii) the derivatives of the summary statistics with respect to the cosmological parameters (and nuisance parameters). In this work, we calculated both directly using the Quijote suite of simulations. Hence it is of vital importance to verify that we are not biased by any numerical effects. We explain how we tested the convergence of both these ingredients in this section.
Fig. A.1. Convergence test of the marginalised onesigma constraints by varying the number of simulations to compute the covariance matrix. The parameters in consideration are shown in the legend. The grey horizontal band shows the region where the variation in the constraints are within 1% compared to the constraint obtained using N_{cov} = 15000. 
The covariance matrix was calculated using 15,000 realisations from the Quijote suite of simulations. However, the data vector of has 2 × 1168 triangular configurations. To verify that our analysis was not affected by numerical effects while calculating the covariance matrix, we varied the number of simulations and checked whether this impacted the final constraints from the Fisher forecast. In Fig. A.1, we showcase the ratio between the marginalised 1σ constraints (for all the cosmological parameters and the nuisance parameter considered in this work) obtained using covariance matrix calculated from N_{cov} simulations, and that obtained using covariance matrix calculated from 15,000 simulations. Using 5000 simulations to compute the covariance matrix leads to 23% on n_{s} and M_{ν} and roughly 1% for Ω_{m}, Ω_{b} and h. However, when using 10,000 simulations or more, it can be seen that there is a clear convergence and a variation well below one percent on all parameter constraints (including the nuisance parameter).
On the top panel of Fig. A.2, we show the ratio of the parameter constraint obtained using N_{deriv} realisations for derivatives to that obtained using N_{deriv} = 500 realisations. We see that only using 300 realisations results in about a 20% difference when compared to the constraints obtained using full 500 simulations. The variations are around 10% and 5% when using 400 and 450 realisations respectively. However, the nuisance parameter b, which is a scaling factor, has the least variation when changing the number of simulations used for computing the derivatives. The variation is only around 2.5% when using 450 simulations. In the bottom panel we show the relative improvement for the cosmological parameters form using over ⟨w_{12}r_{12}⟩_{p}. Using over 450 realisations of the simulations, the improvements remain mostly within 10%.
Fig. A.2. Marginalised onesigma constraints by varying the number of simulations to compute the derivatives (top panel). Bottom panel shows the relative improvement of over ⟨w_{12}r_{12}⟩_{p}. The parameters in consideration are the six cosmological parameters and one nuisance parameter (b) as shown in the legend. The grey horizontal band denotes the region where the variation in the constraints are within 5% compared to the constraint obtained using N_{deriv} = 500. 
Appendix B: Shape dependence of and
In this section, we show how the cosmological parameters affect the mean threepoint relative velocity statistics for the triangular configuration where one side of the leg, i.e. r_{12} is fixed to have a pair separation between 115 and 120 h^{−1}Mpc. Figs. B.1 and B.2 show the shape dependence of and , respectively. Each row corresponds to variation of a single cosmological parameter.
Fig. B.1. Ratio between for varying cosmology and the fiducial cosmology. Each row corresponds to a different cosmological parameter, as given in the bottom left of the first column in each row. The length of the fixed leg of the triangle corresponds to r_{12} ∈ (115, 120) h^{−1}Mpc. 
Fig. B.2. Same as Fig. B.1, but for . 
All Tables
Marginalised 1σ errors for the cosmological parameters and the nuisance parameter (b).
All Figures
Fig. 1. Top panel: measurement (solid line) from the halo catalogues and the analytical prediction (dashed line) using linear perturbation theory as given in Eq. (13). Bottom panel: same for . The plot shows all triangular configurations, in which the minimum and maximum separations for a triangular side are 50 and 120 h^{−1} Mpc, respectively. Light blue vertical band corresponds to triangular configurations where one leg of the triangle is fixed, r_{12} ∈ (100, 105) h^{−1} Mpc. 

In the text 
Fig. 2. Top: effect of M_{ν} on . The length of the fixed leg of the triangle corresponds to r_{12} ∈ (115, 120) h^{−1} Mpc. Each panel corresponds to a different M_{ν}, and it is indicated in the bottom right of each panel (in units of eV). The initial conditions of the fiducial cosmology simulations used in this figure were generated using the Zel’dovich approximation. Bottom: same as top row, but for . 

In the text 
Fig. 3. Derivatives of (left), and (right) with respect to the different cosmological parameters for all triangular configurations, in which the minimum and maximum separation for a triangular side is 20 and 120 h^{−1} Mpc, respectively. The cosmological parameter under consideration is denoted on the right side of each row. 

In the text 
Fig. 4. Left panel shows correlation matrix in the case of halo pairwise velocities, while the right panel shows it for the relative velocities between halo pairs in a triplet computed using the 15 000 realisations of the fiducial run of the Quijote suite of simulations. For the right panel, we used all triangular configurations with 20 ≤ r_{31} ≤ r_{23} ≤ r_{12} ≤ 120 h^{−1} Mpc. 

In the text 
Fig. 5. Joint 68.3% (dark shaded contour) and 95.4% (light shaded contour) credible region for all the pairs of model parameters (six cosmological and one nuisance parameter) at z = 0. Colourcoding is indicated in the legend. 

In the text 
Fig. 6. Joint 68.3% (dark shaded contour) and 95.4% (light shaded contour) credible region for all the pairs of model parameters (six cosmological and one nuisance parameter) at z = 0 for (blue) and (orange), respectively. 

In the text 
Fig. A.1. Convergence test of the marginalised onesigma constraints by varying the number of simulations to compute the covariance matrix. The parameters in consideration are shown in the legend. The grey horizontal band shows the region where the variation in the constraints are within 1% compared to the constraint obtained using N_{cov} = 15000. 

In the text 
Fig. A.2. Marginalised onesigma constraints by varying the number of simulations to compute the derivatives (top panel). Bottom panel shows the relative improvement of over ⟨w_{12}r_{12}⟩_{p}. The parameters in consideration are the six cosmological parameters and one nuisance parameter (b) as shown in the legend. The grey horizontal band denotes the region where the variation in the constraints are within 5% compared to the constraint obtained using N_{deriv} = 500. 

In the text 
Fig. B.1. Ratio between for varying cosmology and the fiducial cosmology. Each row corresponds to a different cosmological parameter, as given in the bottom left of the first column in each row. The length of the fixed leg of the triangle corresponds to r_{12} ∈ (115, 120) h^{−1}Mpc. 

In the text 
Fig. B.2. Same as Fig. B.1, but for . 

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.