© J. Kuruvilla and N. Aghanim 2021

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¹ and future CMB experiments such as the Simons Observatory² (SO, Ade et al. 2019) and CMB-S4³ (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 & 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 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⁴ (da Cunha et al. 2017), and the Widefield ASKAP L-band Legacy All-sky Blind Survey⁵ (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 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á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 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 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 (Villaescusa-Navarro et al. 2020). 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.

2. Data and analysis

2.1. Quijote simulation suite

The Quijote⁶ (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⁻³ Gpc³ and follows the cosmological evolution of 512³ cold dark matter (CDM) particles (and, additionally, in the case of massive neutrino simulations, 512³ 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⁻¹ Mpc) σ₈ = 0.834, and (present-day value of the Hubble constant) H₀ ≡ H(z = 0) = 100 h km s⁻¹ Mpc⁻¹ 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σ₈ = 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 N-body simulations of the Quijote suite. These halos were identified using a friends-of-friends algorithm. We selected halos that have M_h > 5 × 10¹³ h⁻¹ M_⊙ (corresponding to groups and clusters of galaxies) at z = 0, which gives a mean number density of $\overline{n}\sim 0.92\times {10}^{-4}{h}^3$ Mpc⁻³ for the reference simulations. For further details regarding the Quijote suite of simulations, we refer the reader to Villaescusa-Navarro et al. (2020).

2.2. 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)

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

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

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

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 ${\widehat{\mathbf{C}}}^{-1}$, and to obtain it, we computed the covariance matrix of the desired statistics directly from the simulations as follows:

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

where $\overline{\mathit{S}}=\frac{1}{N_{\mathrm{sims}}}\sum _{i=1}^{N_{\mathrm{sims}}}{\mathit{S}}_i$. 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):

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

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, σ₈}, we made use of the central difference approximation:

$$\begin{array}{c}\hfill \frac{\partial \mathit{S}}{\partial \theta }\simeq \frac{\mathit{S}(\theta +\mathrm{d}\theta )-\mathit{S}(\theta -\mathrm{d}\theta )}{2\mathrm{d}\theta }·\end{array}$$(5)

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

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

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

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

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₁₂|r₁₂⟩_p ≃ ⟨δ₁v₂⟩−⟨δ₂v₁⟩. Thus, the mean pairwise velocity can be written as (e.g. Fisher 1995; Juszkiewicz et al. 1998; Reid & White 2011)

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

where j₁(x) = sin(x)/x² − 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₁₂, 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., $\overline{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⁻¹ 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⁻¹ Mpc and above (Shirasaki et al. 2021).

The relative velocity between pairs can be generalised to triplets with separations △₁₂₃ = (r₁₂, r₂₃, r₃₁) 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₁₂|△₁₂₃⟩_t and ⟨w₂₃|△₁₂₃⟩_t (where the subscript t implies that the averages are computed over all particle triplets with separations △₁₂₃). The relative velocity between pair 12 in a triplet, in the single stream fluid approximation, can be written as

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

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:

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

where $\widehat{\mathit{t}}=({\widehat{\mathit{r}}}_{23}-cos\chi {\widehat{\mathit{r}}}_{12})/sin\chi$ and $\chi =arccos({\widehat{\mathit{r}}}_{12}·{\widehat{\mathit{r}}}_{23})$. Thus, we obtain the radial (matter) component as

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

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

$$\begin{array}{cc}\hfill R_{23}({▵}_{123})=\overline{w}(r_{23})& -\frac{1}{2}[\overline{w}(r_{12})cos\chi \hfill \\ \hfill & -\overline{w}(r_{31})\frac{r_{23}+r_{12}cos\chi }{\sqrt{r_{12}^2+r_{23}^2+2r_{12}{r}_{23}cos\chi }}]·\hfill \end{array}$$(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({\overline{M}}_{\mathrm{h}})$, and it can be written as

$$\begin{array}{cc}\hfill R_{12}^{\mathrm{h}}({▵}_{123},M_{\mathrm{h}})& \simeq ⟨(1+b{\delta }_1)(1+b{\delta }_2)(1+b{\delta }_3)({\mathit{v}}_2-{\mathit{v}}_1)⟩\hfill \\ \hfill & =b({\overline{M}}_{\mathrm{h}})R_{12}({▵}_{123}),\hfill \end{array}$$(13)

and similarly

$$\begin{array}{c}\hfill R_{23}^{\mathrm{h}}({▵}_{123},M_{\mathrm{h}})=b({\overline{M}}_{\mathrm{h}})R_{23}({▵}_{123}),\end{array}$$(14)

where the superscript ‘h’ represents the mean statistics for the halos, and in our case ${\overline{M}}_{\mathrm{h}}$ implies all halos with mass greater than 5 × 10¹³ h⁻¹ 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₁₂ ≥ r₂₃ ≥ r₃₁ ≥ 20 h⁻¹ Mpc and separation lengths less than 120 h⁻¹ Mpc. All separations use a bin width of 5 h⁻¹ Mpc and thus have 1168 total triangular configurations.

We directly computed $R_{\mathit{ij}}^{\mathrm{h}}$ from the simulations by measuring the average of the radial components of the relative velocities, i.e. $({\mathit{v}}_j-{\mathit{v}}_i)·{\widehat{\mathit{r}}}_{\mathit{ij}}$, 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_{12}^{\text{h}}$ and $R_{23}^{\text{h}}$ for triangle configurations with r₁₂ ≥ r₂₃ ≥ r₃₁ ≥ 50 h⁻¹ Mpc and separation lengths of less than 120 h⁻¹ 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⁻¹ Mpc -wide bins. Thus, the smallest triangle (triangle configuration ‘0’) considered in Fig. 1 corresponds to r₁₂, r₂₃, r₃₁ ∈ {(50, 55),(50, 55),(50, 55)} h⁻¹ Mpc, while the largest-scale triangle (triangle configuration ‘552’) has r₁₂, r₂₃, r₃₁ ∈ {(115, 120),(115, 120),(115, 120)} h⁻¹ Mpc. In the top and bottom rows, we plot the relative velocity between pairs 12 and 23 in △₁₂₃, respectively. The bias factor as given in Eq. (13) was determined by fitting the prediction using least-squares method, and we obtain b(M_h > 5 × 10¹³) = 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₁₂ and R₂₃ 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₁₂ < 105 h⁻¹ Mpc) and let r₂₃ and r₃₁ 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 left-most data point in the vertical shaded region has its second and third legs of the triangle fixed as r₂₃ = r₃₁ ∈ (50, 55) h⁻¹ Mpc. Thus, as the halo ‘3’ is closer to the pair 12, the mean velocity R₁₂ experiences a greater infall speed. On the other hand, in the case of R₂₃ 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₁₂ and R₂₃ will be equal to each other when r₁₂ = r₂₃, 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_{12}^{\mathrm{h}}$ and $R_{23}^{\mathrm{h}}$, 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⁻¹ 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⁻¹ Mpc, and similarly under-predicts by 3−5% at around 50 h⁻¹ Mpc. As the leading order prediction for $R_{12}^{\mathrm{h}}$ and $R_{23}^{\mathrm{h}}$ 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_{12}^{\mathrm{h}}$ and $R_{23}^{\mathrm{h}}$ over-predict the measurement by 4.6%, and this configuration corresponds to the equilateral triangle with a separation length between 90 and 95 h⁻¹ 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⁻¹ 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_{\mathit{ij}}^{\mathrm{h}}$ 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.

Fig. 1. Top panel: $R_{12}^{\mathrm{h}}$ 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 $R_{23}^{\mathrm{h}}$. 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, r12 ∈ (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 $R_{12}^{\mathrm{h}}$ (top panel) and $R_{23}^{\mathrm{h}}$ (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₁₂ ∈ (115, 120) h⁻¹ Mpc. In each panel, the top right part where r₃₁/r₁₂ = r₂₃/r₁₂ = 1 corresponds to equilateral configuration. Bottom bins (i.e. r₃₁/r₁₂ = r₂₃/r₁₂ ∼ 0.5) contain the folded configurations. Since our analysis has a minimum separation length of 20 h⁻¹ Mpc, we do not explore squeezed configurations (r₁₂ ∼ r₂₃ ≫ r₃₁). 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_{12}^{\mathrm{h}}$ increases. The maximal difference as before is for the equilateral triangle, and minimal when r₃₁ ≪ r₂₃ < r₁₂. In the case of $R_{23}^{\mathrm{h}}$, as seen from the bottom panel, the trend is similar when r₃₁/r₁₂ > 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_{23}^{\mathrm{h}}$ decreases.

Fig. 2. Top: effect of Mν on $R_{12}^{\mathrm{h}}$. The length of the fixed leg of the triangle corresponds to r12 ∈ (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 $R_{23}^{\mathrm{h}}$.

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 $R_{12}^{\mathrm{h}}$ and vice-versa. This overall behaviour is seen with variations in the Hubble parameter also, albeit with a different amplitude. However, for Ω_b and σ₈, the behaviour is reversed, and the amplitude of the change of $R_{12}^{\mathrm{h}}$ with respect to its fiducial counterpart is within 6% for these configurations. For the same triangular configurations, at r₃₁/r₁₂ ≳ 0.4, the variation in cosmological parameters yields a similar response for $R_{23}^{\mathrm{h}}$ as $R_{12}^{\mathrm{h}}$, albeit with a different amplitude. However, this changes in the case of r₃₁/r₁₂ < 0.4, where in the behaviour reverses. The reverse is seen least in the case of variation in σ₈. $R_{12}^{\mathrm{h}}$ and $R_{23}^{\mathrm{h}}$ show the maximal difference in response to cosmological parameters in the case of equilateral triangles and at r₃₁/r₁₂ < 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 σ₈. This should facilitate breaking the M_ν − σ₈ degeneracy in two-point clustering information. To further understand this, we take a look at the derivatives for $R_{12}^{\mathrm{h}}$ (left panel) and $R_{23}^{\mathrm{h}}$ (right panel) in Fig. 3 for all triangular configurations with r₁₂ ≥ r₂₃ ≥ r₃₁ ≥ 20 h⁻¹ Mpc and a maximum separation of 120 h⁻¹ 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 $\partial R_{23}^{\mathrm{h}}/\partial {\sigma }_8$ and $\partial R_{23}^{\mathrm{h}}/\partial M_{\nu }$ looks similar. However, for the small-scale triangles, this is not the case. This points to the possibility that $R_{23}^{\mathrm{h}}+R_{23}^{\mathrm{h}}$ can be used to break the M_ν − σ₈ degeneracy.

Fig. 3. Derivatives of $R_{12}^{\mathrm{h}}$ (left), and $R_{23}^{\mathrm{h}}$ (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 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_{12}^{\mathrm{h}}$ and $R_{23}^{\mathrm{h}}$, this allowed us to compute the covariance matrix reliably. We show the cross-correlation coefficient (${\rho }_{\mathit{ij}}=C_{\mathit{ij}}/\sqrt{C_{\mathit{ii}}{C}_{\mathit{jj}}}$) 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 non-diagonal contributions. We especially notice that $R_{12}^{\mathrm{h}}$ and $R_{23}^{\mathrm{h}}$ 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 ≤ r31 ≤ r23 ≤ r12 ≤ 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⁻¹ 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 ${⟨{\mathit{w}}_{12}^{\mathrm{h}}⟩}_{\mathrm{p}}$, $R_{12}^{\mathrm{h}}$, and $R_{23}^{\mathrm{h}}$. For r_min  =  20 h⁻¹ Mpc, we find that $R_{12}^{\mathrm{h}}+R_{23}^{\mathrm{h}}$ has a 1σ constraints for {Ω_m, Ω_b, h, n_s, σ₈, M_ν} = {0.0091, 0.0024, 0.0226, 0.0221, 0.0156, 0.0655 eV}, which presents an improvement over constraints from ${⟨{\mathit{w}}_{12}^{\mathrm{h}}⟩}_{\mathrm{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_{12}^{\mathrm{h}}+R_{23}^{\mathrm{h}}$ 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. Colour-coding 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 & 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_{12}^{\mathrm{h}}$ and $R_{23}^{\mathrm{h}}$ 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 $R_{12}^{\mathrm{h}}$ and $R_{23}^{\mathrm{h}}$, 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⁻¹ Mpc and 120 h⁻¹ Mpc. In one case, we considered the constraints from $R_{12}^{\mathrm{h}}$ and $R_{23}^{\mathrm{h}}$ separately, and in the second case we consider the joint $R_{12}^{\mathrm{h}}+R_{23}^{\mathrm{h}}$ 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_{12}^{\mathrm{h}}$ alone, while the orange contour shows that from $R_{12}^{\mathrm{h}}+R_{23}^{\mathrm{h}}$. In the figure, we have not shown the joint credible regions from $R_{23}^{\mathrm{h}}$ because it is the same as $R_{12}^{\mathrm{h}}$. We can see that on all parameters $R_{12}^{\mathrm{h}}+R_{23}^{\mathrm{h}}$ is able to obtain tighter constraints. Quantitatively, they are 1.34−1.46 times better than those considering only $R_{12}^{\mathrm{h}}$ or $R_{23}^{\mathrm{h}}$.

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 $R_{12}^{\mathrm{h}}$ (blue) and $R_{12}^{\mathrm{h}}+R_{23}^{\mathrm{h}}$ (orange), respectively.

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 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¹³ h⁻¹ M_⊙, while in Hahn et al. (2020) they considered M_h > 3.2 × 10¹³ h⁻¹ M_⊙). Using the monopole and quadrupole of the redshift-space power spectrum for k_max = 0.5 h⁻¹ 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 $R_{12}^{\mathrm{h}}+R_{23}^{\mathrm{h}}$ (for r_min = 20 h⁻¹ 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, σ₈), respectively.

Hahn et al. (2020) also provided the cosmological constraints from using the redshift-space bispectrum monopole (the three-point clustering statistics in Fourier space). Our constraint on M_ν from the three-point relative velocity statistics is slightly less than the 1-σ₈ 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 three-point mean relative velocity when compared to those from bispectrum monopole with a factor of (1.2, 1.7, 1.7, 1.5), respectively. Albeit σ₈ produces a slight decrease in the constraining power with a factor of 0.9. This shows the potential of three-point 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 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⁻¹ Mpc and r_max = 120 h⁻¹ 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 $R_{12}^{\mathrm{h}}$ and $R_{23}^{\mathrm{h}}$ 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 three-point mean relative velocity statistics, we used the Fisher-matrix 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 three-point 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 $R_{12}^{\mathrm{h}}+R_{23}^{\mathrm{h}}$. In the case of pairwise velocities, we used pair separations from 20 h⁻¹ Mpc to 120 h⁻¹ Mpc, with a bin width of 5 h⁻¹ Mpc. In the case of $R_{12}^{\mathrm{h}}$ and $R_{23}^{\mathrm{h}}$, we considered all triangular configurations with a minimum triangle side length of 20 h⁻¹ Mpc and a maximum of 120 h⁻¹ Mpc. This amounted to using 1168 total triangular configurations. $R_{12}^{\mathrm{h}}+R_{23}^{\mathrm{h}}$ 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; Gonzalez-Garcia et al. 2016; Capozzi et al. 2017; de Salas et al. 2018). Thus, the three-point mean relative velocity statistics could help us to pinpoint the summed mass of neutrinos in the future. Similarly, for the other cosmological parameters, the three-point 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 $R_{12}^{\mathrm{h}}$ and $R_{23}^{\mathrm{h}}$, being a three-point 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 two-point and three-point 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 $R_{12}^{\mathrm{h}}+R_{23}^{\mathrm{h}}$ 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 Villaescusa-Navarro 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 ERC-2015-AdG 695561. We are thankful to the community for developing and maintaining open-source software packages extensively used in our work, namely CYTHON (Behnel et al. 2011), MATPLOTLIB (Hunter 2007), and NUMPY (Harris et al. 2020).

References

Appendix A: 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 one-sigma 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 Ncov = 15000.

The covariance matrix was calculated using 15,000 realisations from the Quijote suite of simulations. However, the data vector of $R_{12}^{\mathrm{h}}+R_{23}^{\mathrm{h}}$ 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 2-3% 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 $R_{12}^{\mathrm{h}}+R_{23}^{\mathrm{h}}$ over ⟨w₁₂|r₁₂⟩_p. Using over 450 realisations of the simulations, the improvements remain mostly within 10%.

Fig. A.2. Marginalised one-sigma constraints by varying the number of simulations to compute the derivatives (top panel). Bottom panel shows the relative improvement of $R_{12}^{\mathrm{h}}+R_{23}^{\mathrm{h}}$ over ⟨w12|r12⟩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 Nderiv = 500.

Appendix B: Shape dependence of $R_{12}^{\mathrm{h}}$ and $R_{23}^{\mathrm{h}}$

In this section, we show how the cosmological parameters affect the mean three-point relative velocity statistics for the triangular configuration where one side of the leg, i.e. r₁₂ is fixed to have a pair separation between 115 and 120 h⁻¹Mpc. Figs. B.1 and B.2 show the shape dependence of $R_{12}^{\mathrm{h}}$ and $R_{23}^{\mathrm{h}}$, respectively. Each row corresponds to variation of a single cosmological parameter.

Fig. B.1. Ratio between $R_{12}^{\mathrm{h}}({▵}_{123})$ 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 r12 ∈ (115, 120) h−1Mpc.

Fig. B.2. Same as Fig. B.1, but for $R_{23}^{\mathrm{h}}({▵}_{123})$.

All Tables

All Figures

Fig. 1. Top panel: $R_{12}^{\mathrm{h}}$ 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 $R_{23}^{\mathrm{h}}$. 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, r12 ∈ (100, 105) h−1 Mpc.

In the text

Fig. 2. Top: effect of Mν on $R_{12}^{\mathrm{h}}$. The length of the fixed leg of the triangle corresponds to r12 ∈ (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 $R_{23}^{\mathrm{h}}$.

In the text

	Fig. 3. Derivatives of $R_{12}^{\mathrm{h}}$ (left), and $R_{23}^{\mathrm{h}}$ (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 ≤ r31 ≤ r23 ≤ r12 ≤ 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. Colour-coding 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 $R_{12}^{\mathrm{h}}$ (blue) and $R_{12}^{\mathrm{h}}+R_{23}^{\mathrm{h}}$ (orange), respectively.
In the text

	Fig. A.1. Convergence test of the marginalised one-sigma 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 Ncov = 15000.
In the text

Fig. A.2. Marginalised one-sigma constraints by varying the number of simulations to compute the derivatives (top panel). Bottom panel shows the relative improvement of $R_{12}^{\mathrm{h}}+R_{23}^{\mathrm{h}}$ over ⟨w12|r12⟩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 Nderiv = 500.

In the text

	Fig. B.1. Ratio between $R_{12}^{\mathrm{h}}({▵}_{123})$ 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 r12 ∈ (115, 120) h−1Mpc.
In the text

	Fig. B.2. Same as Fig. B.1, but for $R_{23}^{\mathrm{h}}({▵}_{123})$.
In the text