# Cosmology with the kinetic Sunyaev-Zeldovich effect: Independent of the optical depth and $\sigma_8$

Joseph Kuruvilla

► To cite this version:

Joseph Kuruvilla. Cosmology with the kinetic Sunyaev-Zeldovich effect: Independent of the optical depth and $\sigma_8$. Astronomy and Astrophysics - A&A, 2022, 660, pp.A113. 10.1051/0004-6361/202142325. [hal-03377469](https://hal.science/hal-03377469)

HAL Id: hal-03377469

https://hal.science/hal-03377469

Submitted on 31 Mar 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Cosmology with the kinetic Sunyaev–Zeldovich effect: Independent of the optical depth and $\sigma_8$

Joseph Kuruvilla

Université Paris-Saclay, CNRS, Institut d'Astrophysique Spatiale, 91405 Orsay, France
e-mail: joseph.kuruvilla@universite-paris-saclay.fr

Received 28 September 2021 / Accepted 3 February 2022

## ABSTRACT

Cosmological constraints obtained by the kinetic Sunyaev–Zeldovich experiments are degenerate with the optical depth measurement – an effect that is commonly known as the optical-depth degeneracy. In this work, we introduce a new statistic based on the first moment of relative velocity between pairs in a triplet, which is capable of constraining cosmological parameters independently of the optical depth and of $\sigma_8$. Using 22 000 N-body simulations from the Quijote suite, we quantified the information content in the new statistic using Fisher matrix forecast. We find that it is able to obtain strong constraints on the cosmological parameters, particularly on the summed neutrino mass. The constraints bring an improvement on all cosmological model parameters by a factor of 6.2–12.9 and 2.3–5.7 when compared to those obtained from the mean pairwise velocity and from the redshift-space halo power spectrum, respectively. Thus, this new statistic paves a way forward in constraining cosmological parameters independent of the optical depth and $\sigma_8$ when using data from future kinetic Sunyaev–Zeldovich experiments alone.

**Key words.** large-scale structure of Universe – cosmology: theory – cosmic background radiation

## 1. Introduction

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

All the while, stronger constraints on the summed neutrino mass ($M_ν$) can be obtained by combining various cosmological probes since the massive neutrinos leave an imprint on various cosmological observables (e.g., Wong 2011; Lesgourgues & Pastor 2012). Yet these constraints have an additional model dependence. Assuming the standard ‘lambda cold dark matter’ (ΛCDM) model, one of the strongest constraints on the summed neutrino mass has been obtained by combining the cosmic microwave background (CMB, Planck Collaboration VI 2020), baryonic acoustic oscillation, and redshift-space galaxy clustering (eBOSS Collaboration 2021) to obtain an upper limit of $M_ν < 0.102$ eV. However, by considering extensions of the standard cosmological model, the upper limit becomes less stringent (e.g., Vagnozzi et al. 2018; Choudhury & Hannestad 2020). The community can expect the summed neutrino mass to be measured with increased precision from cosmological probes in the foreseeable future with the advent of next generation of

CMB surveys (e.g., the Simons Observatory<fnref n="1" /> (SO, Ade et al. 2019), CMB-S4² (Abazajian et al. 2016)) as well as the stage IV galaxy redshift surveys (e.g., the ‘Dark Energy Spectroscopic Instrument’<fnref n="3" /> (DESI, DESI Collaboration 2016), Euclid<fnref n="4" /> (Laureijs et al. 2011) and Nancy Grace Roman space telescope<fnref n="5" /> (Spergel et al. 2015)).

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

<footnote n="1">¹ https://simonsobservatory.org/</footnote>

<footnote n="2">² https://cmb-s4.org/</footnote>

<footnote n="3">³ https://www.desi.lbl.gov/</footnote>

<footnote n="4">⁴ https://www.euclid-ec.org/</footnote>

<footnote n="5">⁵ https://roman.gsfc.nasa.gov/</footnote>
the one-point probability distribution function of the matter density (e.g., Uhlemann et al. 2020) and the void size function (e.g., Bayer et al. 2021).

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

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

## 2. Cosmology using the kSZ effect

### 2.1. Kinetic Sunyaev-Zeldovich effect

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

$$ \frac{\Delta T(\hat{n})}{T_{\text{cmb}}}|_{\text{kSZ}} = - \int dl \sigma_{\text{T}} \left( \frac{\mathbf{v}_{\text{e}} \cdot \hat{n}}{c} \right) n_{\text{e}}, \\ = -\tau \left( \frac{\mathbf{v}_{\text{e}} \cdot \hat{n}}{c} \right), \qquad (1) $$

where $\sigma_{\text{T}}$ is the Thomson scattering cross-section, $T_{\text{cmb}}$ is the CMB temperature, c is the speed of light, $\mathbf{v}_{\text{e}}$ is the peculiar velocity of free electrons, and $n_{\text{e}}$ is the physical free electron number density. The integral $\int dl$ is computed along the LOS which is given by $\hat{n}$. The optical depth is defined as $\tau = \int dl \sigma_{\text{T}} n_{\text{e}}$, that is, the integrated electron density.

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

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

$$ \frac{\langle \Delta T^{\text{kSZ}}(r_{12}) \rangle}{T_{\text{cmb}}} \approx -\tau \frac{\bar{w}(r_{12})}{c}, \qquad (2) $$

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

$$ \langle w_{12} | \mathbf{r}_{12} \rangle_{\text{p}} = \frac{\langle (1 + \delta_1)(1 + \delta_2)(\mathbf{v}_2 - \mathbf{v}_1) \rangle}{\langle (1 + \delta_1)(1 + \delta_2) \rangle}, \qquad (3) $$

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

$$ \langle w_{12} | \mathbf{r}_{12} \rangle_{\text{p}} = \bar{w}(r_{12}) \hat{\mathbf{r}}_{12} \approx -\frac{f}{\pi^2} \hat{\mathbf{r}}_{12} \int_0^{\infty} k \, j_1(kr_{12}) P(k) dk, \qquad (4) $$

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

$$ \Delta T^{\text{kSZ}} \propto \tau f \sigma_{8}^{2}, \qquad (5) $$

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

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

In the previous section, we mentioned about the mean relative velocity between two tracers (i.e., between a pair) or the mean pairwise velocity. However, this can be generalised to the case of three tracers, in which we can consider two mean relative velocity between pairs in a triplet with separations $\Delta_{123} = (r_{12}, r_{23}, r_{31})$: (i) $\langle w_{12} | \Delta_{123} \rangle_{\text{t}}$, and (ii) $\langle w_{23} | \Delta_{123} \rangle_{\text{t}}$. The subscript t here implies that the averages are computed over all triplets with separations $(r_{12}, r_{23}, r_{31})$. Similarly to Eq. (3), in the single
stream fluid approximation, the mean relative velocity between pair 12 in a triplet can be written as (Kuruvilla & Porciani 2020):

$$ 
\begin{aligned}
\langle \mathbf{w}_{12}|\Delta_{123} \rangle_t &= \frac{\langle (1 + \delta_1)(1 + \delta_2)(1 + \delta_3)(\mathbf{v}_2 - \mathbf{v}_1) \rangle}{\langle (1 + \delta_1)(1 + \delta_2)(1 + \delta_3) \rangle} \\
&\approx \langle \delta_1 \mathbf{v}_2 \rangle - \langle \delta_2 \mathbf{v}_1 \rangle + \langle \delta_3 \mathbf{v}_2 \rangle - \langle \delta_3 \mathbf{v}_1 \rangle \\
&= \bar{w}(r_{12}) \hat{\mathbf{r}}_{12} - \frac{1}{2} [\bar{w}(r_{23}) \hat{\mathbf{r}}_{23} + \bar{w}(r_{31}) \hat{\mathbf{r}}_{31}] .
\end{aligned}
\quad (6) $$

The three-point mean relative velocity statistics can be composed into both its radial ($R_{ij}$) and transverse ($T_{ij}$) component in the plane of the triangle defined by the particles. This is in contrast to the mean pairwise velocity for which the transverse component is zero. In the case of $\langle \mathbf{w}_{12}|\Delta_{123} \rangle_t$, it follows as:

$$ 
\begin{aligned}
\langle \mathbf{w}_{12}|\Delta_{123} \rangle_t &= \langle \mathbf{w}_{12} \cdot \hat{\mathbf{r}}_{12}|\Delta_{123} \rangle_t \hat{\mathbf{r}}_{12} + \langle \mathbf{w}_{12} \cdot \hat{\mathbf{t}}|\Delta_{123} \rangle_t \hat{\mathbf{t}} \\
&= R_{12}(\Delta_{123}) \hat{\mathbf{r}}_{12} + T_{12}(\Delta_{123}) \hat{\mathbf{t}},
\end{aligned}
\quad (7) $$

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

$$ R_{12}(\Delta_{123}) = \bar{w}(r_{12}) - \frac{1}{2} \left[ \bar{w}(r_{23}) \cos \chi - \bar{w}(r_{31}) \frac{r_{12} + r_{23} \cos \chi}{\sqrt{r_{12}^2 + r_{23}^2 + 2r_{12}r_{23} \cos \chi}} \right]. \quad (8) $$

Similarly, the mean radial relative velocity between the pair 23 in $\Delta_{123}$ can be written as

$$ R_{23}(\Delta_{123}) = \bar{w}(r_{23}) - \frac{1}{2} \left[ \bar{w}(r_{12}) \cos \chi - \bar{w}(r_{31}) \frac{r_{23} + r_{12} \cos \chi}{\sqrt{r_{12}^2 + r_{23}^2 + 2r_{12}r_{23} \cos \chi}} \right]. \quad (9) $$

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

$$ \frac{\Delta T_{ij}^{\text{kSZ}}(\Delta_{123})}{T_{\text{cmb}}} \approx -\tau \frac{R_{ij}^{\text{h}}(\Delta_{123})}{c}, \quad (10) $$

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

$$ \mathcal{R}(\Delta_{123}) = \frac{\langle \mathbf{w}_{12} \cdot \hat{\mathbf{r}}_{12}|\Delta_{123} \rangle_t}{\langle \mathbf{w}_{23} \cdot \hat{\mathbf{r}}_{23}|\Delta_{123} \rangle_t}, \quad (11) $$

which tells us how quickly the average infall velocity of pair 12 performs in comparison to the average infall velocity of pair 23 for a specific triangular configuration $\Delta_{123}$. On linear scales, using perturbation theory at LO, $\mathcal{R}(\Delta_{123})$ can be written as the ratio between Eqs. (8) and (9)

$$ \mathcal{R}(\Delta_{123}) = \frac{R_{12}(\Delta_{123})}{R_{23}(\Delta_{123})} = \frac{R_{12}^{\text{h}}(\Delta_{123})}{R_{23}^{\text{h}}(\Delta_{123})} \equiv \frac{\Delta T_{12}^{\text{kSZ}}(\Delta_{123})}{\Delta T_{23}^{\text{kSZ}}(\Delta_{123})}. \quad (12) $$

The above statistic is thus independent of optical depth, $\sigma_8$, and linear bias. In the following sections, we take a detailed look at whether the Ansatz of $\sigma_8$ and linear bias independence holds up. Additionally, we study the cosmological information content in $\mathcal{R}(\Delta_{123})$.

# 3. Data and analysis

## 3.1. Quijote simulation suite

In this work, we make use of the Quijote<fnref n="6" /> (Villaescusa-Navarro et al. 2020) suite of simulations, which was run using the tree-PM code GADGET-3 (Springel 2005). Spanning more than a few thousand cosmological models, it contains 44 100 N-body simulations. These simulations have a box length of $1 h^{-1}$ Gpc, tracking the evolution of 512³ cold dark matter (CDM) particles. The initial conditions (ICs) were generated at redshift $z = 127$ using the second-order Lagrangian perturbation theory. The fiducial cosmological parameters (assuming zero-summed neutrino mass) for the simulation is as follows: the total matter density: $\Omega_m = 0.3175$; baryonic matter density: $\Omega_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: $\sigma_8 = 0.834$; and the present-day value of the Hubble constant: $H_0 \equiv H(z=0) = 100 h \text{ km s}^{-1} \text{ Mpc}^{-1}$, with $h = 0.6711$. This is broadly consistent with the Planck 2018 result (Planck Collaboration VI 2020). The suite consists of 15 000 random realisations for the fiducial cosmology. For the purpose of calculating derivatives, Quijote provides a set of 500 random realisations wherein only one parameter is varied with respect to the fiducial cosmology. The variations are as follows: $\{\Omega_m^+, \Omega_m^-, \Omega_b^+, \Omega_b^-, n_s^+, n_s^-, \sigma_8^+, \sigma_8^-\} = \{0.3275, 0.3075, 0.051, 0.047, 0.9824, 0.9424, 0.849, 0.819\}$ and $\{h^+, h^-\} = \{0.6911, 0.6511\}$.

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

In this work, we use the halo catalog data from 22 000 N-body simulations of the Quijote suite. These halos were identified using a friends-of-friends algorithm. We selected halos that have a halo mass of $M_h > 5 \times 10^{13} h^{-1} M_\odot$ (corresponding to groups and clusters of galaxies) at $z = 0$, which gives a mean number density of $\bar{n} \sim 0.92 \times 10^{-4} h^3 \text{ Mpc}^{-3}$ for the reference simulations. Additionally, in the case of (i) fiducial cosmology, and (ii) for variations in $\sigma_8$ (both $\sigma_8^+$ and $\sigma_8^-$), we used 30 realisations of the particle data (randomly down-sampled to $100^3$ particles) to compute $\mathcal{R}(\Delta_{123})$.

## 3.2. Fisher-matrix formalism

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

$$ F_{\alpha\beta} = \left\langle -\frac{\partial^2 \ln \mathcal{L}}{\partial \theta_{\alpha} \partial \theta_{\beta}} \right\rangle, \quad (13) $$

where $\theta_\alpha$ and $\theta_\beta$ are two of the cosmological model parameters and $\mathcal{L}$ is the likelihood of the data given a model. Assuming a Gaussian likelihood, we can express the Fisher information matrix as

$$ F_{\alpha\beta} = \frac{\partial \mathcal{R}}{\partial \theta_{\alpha}} \cdot \hat{\mathbf{C}}^{-1} \cdot \frac{\partial \mathcal{R}^{\top}}{\partial \theta_{\beta}}, \quad (14) $$

<footnote n="6">⁶ https://quijote-simulations.readthedocs.io/</footnote>
where $\mathcal{R}$ represents the data vector for the ratio statistic we introduced in Eq. (11) and $\hat{\mathbf{C}}^{-1}$ is the precision matrix (i.e., the inverse covariance matrix). It should be noted that in the definition of $F_{\alpha\beta}$, we neglected a term that appears due to the cosmology dependence of the covariance matrix. However, the correction has been shown to have a negligible effect (Kodwani et al. 2019). We compute the covariance matrix of $\mathcal{R}$ directly from the simulations as follows:

$$ \tilde{\mathbf{C}} = \frac{1}{N_{\text{sims}} - 1} \sum_{i=1}^{N_{\text{sims}}} (\mathbf{R}_i - \bar{\mathbf{R}})(\mathbf{R}_i - \bar{\mathbf{R}})^{\text{T}}, \quad (15) $$

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

$$ \hat{\mathbf{C}}^{-1} = \frac{N_{\text{sims}} - N_{\text{bins}} - 2}{N_{\text{sims}} - 1} \tilde{\mathbf{C}}^{-1}, \quad (16) $$

where $N_{\text{bins}}$ is the number of bins in $\mathcal{R}$.

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

$$ \frac{\partial \mathcal{R}}{\partial \theta} \approx \frac{\mathcal{R}(\theta + d\theta) - \mathcal{R}(\theta - d\theta)}{2 d\theta}. \quad (17) $$

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

$$ \frac{\partial \mathcal{R}}{\partial M_{\nu}} \approx \frac{-\mathcal{R}(M_{\nu} = 0.4) + 4\mathcal{R}(M_{\nu} = 0.2) - \mathcal{R}(M_{\nu} = 0)}{0.4}. \quad (18) $$

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

## 4. Results

In Fig. 1, we show the direct measurements of $\mathcal{R}(\Delta_{123})$ from the 15 000 reference halo catalogs (solid blue line) and compare it with the LO predictions (dashed orange line). We consider all triangular configurations with $r_{\text{min}} \in (40, 45)$ and $r_{\text{max}} \in (115, 120)$, and such that $r_{12} \ge r_{23} \ge r_{31}$. All the separation scales have a bin width of $5 h^{-1}$ Mpc. It thus corresponds to a total of 766 triangular configurations, spanning from configuration ‘0’ being the smallest (i.e., $\Delta_{123} \in \{(40, 45), (40, 45), (40, 45)\} h^{-1}$ Mpc) to configuration ‘765’ being the largest ($\Delta_{123} \in \{(115, 120), (115, 120), (115, 120)\} h^{-1}$ Mpc). In Appendix A, we show the

<figure><img src="image_2.png" /><figcaption>Fig. 1. Comparison of theoretical prediction (orange dashed line) for $\mathcal{R}(\Delta_{123})$ using perturbation theory at LO against the direct measurement (blue solid line) from the halo catalogs of the Quijote suite of simulation is shown in the top panel. While the residual showing the deviation of the theoretical prediction from the direct measurement from the simulations is shown in the bottom panel. The blue shaded region denotes the 5% region. The triangle configuration go from the smallest being $\{(40, 45), (40, 45), (40, 45)\} h^{-1}$ Mpc to the largest which corresponds to $\{(115, 120), (115, 120), (115, 120)\} h^{-1}$ Mpc. The vertical dashed lines correspond to a change (increase) in the value of $r_{12}$ bin in $\Delta_{123}$.</figcaption></figure>

comparison between the theoretical predictions and the direct measurements when $r_{12}$ is fixed to $(115, 120) h^{-1}$ Mpc for the triangular configurations. We can see from Eqs. (8) and (9) that the LO prediction for the mean three-point relative velocities ($R_{12}$ and $R_{23}$) will be equal to each other when $r_{12} = r_{23}$, irrespective of the length of the third side. This is directly visible in Fig. 1, where $\mathcal{R}^{\text{theory}} = 1$ when this condition is met. When comparing the theoretical predictions with the direct measurements from the halo catalogs, we see that it is generally accurate within 4–5% for configurations with all separation lengths greater than $55 h^{-1}$ Mpc. As is expected when the separation length decreases, the fidelity of the LO prediction also decreases with the maximum deviation at about 27% for the triangular configuration $\{(100, 105), (50, 55), (50, 55)\} h^{-1}$ Mpc. This motivates us to directly measure $\mathcal{R}$ from the simulations to compute the derivatives for the Fisher information matrix.

In addition, $\mathcal{R}(\Delta_{123})$ is unaffected by variation in $\sigma_8$ as mentioned earlier in Sect. 2.2. To demonstrate this, we computed the ratio of $\mathcal{R}(\Delta_{123})$ for $\sigma_8^+ = 0.849$ and $\sigma_8^- = 0.819$ using dark matter particles from 30 realisations and showcase this result in Fig. 2. The (blue) dot represents the mean of the measurement, while the scatter is shown using the (orange) error bars. Thus we can conclude that $\mathcal{R}$ is independent of $\sigma_8$. Similarly in Fig. 3, we take a look at the bias dependence of $\mathcal{R}(\Delta_{123})$ (black solid line), and $R_{12}$ (blue dashed line), where we show the ratio of each between the halo and matter component for each of the summary statistics. The bias term for the mean relative velocity between pair ‘23’ in a triplet is similar to $R_{12}^h$, and is thus not shown in the figure. As reported in Kuruvilla & Aghanim (2021), for these triangular configurations (assuming a scale independent bias), it yields a bias factor around 1.85 for $R_{12}$ and $R_{23}$. For the purpose of computing these ratios in the figure, we used the mean relative velocity information for matter from 30 realisations of the dark-matter only simulations; and for the halo, we utilised the 15 000 catalogues. The shaded regions represent the $1\sigma$ errors from the propagation of uncertainties of the mean relative velocity
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Ratio of $R$ at $\sigma_8^+ = 0.849$ to $R$ at $\sigma_8^- = 0.819$. The (blue) solid dot, and the (orange) error bar represents the mean and the relative error on the mean, respectively.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Dashed (blue) line: bias for $R_{12}$, i.e. the mean radial relative velocity between pairs 1 and 2 in a triplet. The solid (black) line shows the (weak) bias dependence of the ratio statistics $R$, which is equal to one within 1–2% for all triangular configurations considered here.</figcaption></figure>

statistics for matter and the halo. We can see that on large separation scales, the newly introduced statistic (black solid line) is bias-independent, while for the smallest triangle configuration, there is a very weak dependence of bias when considering the newly introduced statistic. This thus supports the Ansatz presented in Eq. (12), where the LO in perturbation theory renders $R$ to be bias independent on linear scales. For all triangular configurations considered in this work, the bias is found to be equal to one within 1–2%; hence, for the purposes of Fisher matrix formalism, we consider $R$ being independent of a (constant) linear bias term.

Since $R$ is found to be independent of $\sigma_8$ at the scales we are probing (i.e., $r_{min} \ge 40 \, h^{-1} \, \text{Mpc}$), this sets the statistic in the unique position of being unaffected by the degeneracy in the $M_v-\sigma_8$ parameter plane. We checked the impact of summed neutrino mass on $R$ utilising three non-zero neutrino mass in Fig. 4. The solid (blue) line shows the impact of $M_v = 0.1 \, \text{eV}$ on $R$ when compared to zero neutrino-mass cosmology. Similarly, the dashed (orange) and dash-dotted (green) lines showcases the impact of $M_v = 0.2$ and $M_v = 0.4 \, \text{eV}$, respectively.

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Effect of summed neutrino mass on $R(\Delta_{123})$, as measured directly from the simulations, when compared to zero neutrino mass fiducial cosmology. The summed neutrino mass considered here are denoted in the legend, and its units are given in eV.</figcaption></figure>

As can be seen, when the neutrino mass increases there is a decrease in the infall velocity between a pair in most of the triangular configurations. This is related to the free-streaming of neutrinos as a result of them having large thermal velocities. As a result, below the free-streaming scale, neutrinos do not cluster, which further slows down the collapse of the matter overall. This leads to an overall reduction in the growth of overall density perturbations at scales below free-streaming scale, and thus causes a suppression of power on large Fourier modes when looking at the matter power spectrum (e.g., Wong 2011; Lesgourgues & Pastor 2012). When looking at $R$ for all the configurations we measured, the maximal effect of suppression is seen in the case when $M_v = 0.4 \, \text{eV}$, and for the triangular configuration {(100, 105), (50, 55), (50, 55)} $h^{-1} \, \text{Mpc}$ when compared to the zero neutrino mass cosmology.

## 4.1. Cosmological parameters

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

We now take a look at the information content in $R$ using the Fisher information matrix formalism, as defined in Sect. 3.2. As mentioned, we computed both the elements of it directly from the simulations. Since the bias dependence of $R$ was shown to being very weak even at small scales (~40–50 $h^{-1} \, \text{Mpc}$) and independent at large scales (≥80 $h^{-1} \, \text{Mpc}$), we did not consider the bias parameter in the Fisher-matrix formalism. Thus, the model parameters are $\Omega_m, \Omega_b, h, n_s$, and $M_v$. We show the results of our Fisher forecast obtained from $R$ (orange colour) in Fig. 6, where the contour denote the 68.3% joint credible region for all possible model parameters. The $1\sigma$ marginalised error for any
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Correlation matrix (i.e., the covariance matrix of $\mathcal{R}$ normalised by its diagonal elements) computed using 15 000 realisations of the Qujote simulations. The triangle configurations are same as in Fig. 1.</figcaption></figure>

model parameter $\theta_\alpha$ is given by $\sqrt{F_{\alpha\alpha}^{-1}}$, and they are as follows for the parameters we considered: $\{\sigma_{\Omega_m}, \sigma_{\Omega_b}, \sigma_{h}, \sigma_{n_s}, \sigma_{M_v}\} \equiv \{0.0158, 0.0041, 0.0391, 0.0394, 0.1175\}$ (where $\sigma_{M_v}$ is given in units of eV).

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

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

## 4.2. Discussions

It is useful to consider how the constraints from $\mathcal{R}$ fares against those obtained from clustering statistics. In order to answer

that question, we compare the constraints obtained in this work with those obtained from the redshift space halo power spectrum and the halo bispectrum in Hahn et al. (2020). Compared with the constraints from the redshift-space power multipoles for $k_{\max} = 0.2 \, h \, \text{Mpc}^{-1}$ (which is closest to the $r_{\min}$ considered in this work), $\mathcal{R}$ obtains a factor of improvement of $\{2.3, 3.6, 4.5, 5.4, 5.7\}$ for $\{\Omega_m, \Omega_b, h, n_s, M_v\}$, respectively. However it is slightly reduced to $\{1.4, 2.9, 3.1, 3.3, 2.5\}$ when comparing with the constraints from power spectrum when $k_{\max} = 0.5 \, h \, \text{Mpc}^{-1}$. This improvement over the power spectrum (a two-point summary statistics) is not surprising as $\mathcal{R}$ is based on the first moment of the three-point relative velocity statistics. Hence, it is interesting to compare the constraints against those obtained from the redshift-space bispectrum monopole, and when using all triangular configurations for $k_{\max} = 0.2 \, h \, \text{Mpc}^{-1}$, $\mathcal{R}$ still enjoys an improvement by a factor of $\{1.8, 2.9, 3.2, 3.1, 1.8\}$ for $\{\Omega_m, \Omega_b, h, n_s, M_v\}$, respectively. However, when considering a larger set of triangular configurations for the bispectrum with $k_{\max} = 0.5 \, h \, \text{Mpc}^{-1}$, there is less constraining power for $\mathcal{R}$ with a factor of $\{0.7, 1.0, 1.0, 0.9, 0.4\}$ for $\{\Omega_m, \Omega_b, h, n_s, M_v\}$, respectively. This is not surprising, as the bispectrum monopole in this case probes further into the nonlinear scales, whereas $\mathcal{R}$ was analysed for triangular configurations with separation scales of $40 \, h^{-1} \, \text{Mpc}$ and above.

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

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

Another caveat that we do not discuss in this work is the mass dependence of the optical depth parameter, which has been shown to increase as the halo mass increases (e.g., Battaglia 2016). We considered it as an averaged quantity, as shown in Eqs. (2) and (10). However, we do not envision the mass dependence to affect $\mathcal{R}(\Delta_{123})$, as long as the ratio is measured using the same mass bin. On the other hand, in assuming a fixed cosmology, we could consider a scenario where measuring $R_{23}^h$ is fixed to a high halo mass bin while measuring $R_{12}^h$ for various mass bins in Eq. (11). Thus, it could lead to making potential
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Joint 68.3% credible region for all the pairs of cosmological model parameters at $z = 0$. Different colours indicate the different summary statistics: $R_{12}^h + R_{23}^h$ (blue contour, taken from Kuruvilla & Aghanim 2021) and $R$ (orange contour). The fiducial values are indicated by the black dashed lines, and $M_v$ is given in units of eV.</figcaption></figure>

measurements of the (scaled) mass dependence of optical depth (and degenerate with the bias factor) from direct kSZ experiment directly.

**5. Conclusions**

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

In this paper, we extend the application of the mean three-point relative velocity statistics and we introduce a new ratio statistic, $R$, in Eq. (11), which is unaffected by $\sigma_8$. This enables us to constrain the neutrino mass, in addition to other cosmological parameters, independently of $\sigma_8$. Moreover, in the context of kSZ experiments, this statistic is independent of optical depth, hence it circumvents the optical depth degeneracy which currently acts a limiting factor in the determination of cosmological parameters from the kSZ experiments. Furthermore, the leading order perturbation theory prediction suggests that $R$ will be bias-independent on linear scales. We verified this assumption by measuring $R$ values for both halos and matter, finding

that the bias is consistent with one at 1–2% for all triangular configurations we probed in this work ($r_{\text{min}} = 40 h^{-1}$ Mpc and $r_{\text{max}} = 120 h^{-1}$ Mpc).

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

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

In summary, we have introduced a new statistic based on the mean radial relative velocity between pairs in a triplet and we have shown that it can act as robust cosmological observable that could lead to sizeable information gain in comparison to the mean radial pairwise velocity. One of the limitations of the kSZ experiments is the optical depth degeneracy and, thus, breaking this degeneracy requires some form of an external data set (for e.g., using fast radio bursts as suggested in Madhavacheril et al. 2019). This new statistic thus provides a way forward in which the cosmological parameters can be constrained using data from future kinetic Sunyaev–Zeldovich experiments alone, without being affected by the optical depth parameter.
**Acknowledgements.** We would like to thank the referee for their comments, Nabila Aghanim and Francisco Villaescusa-Navarro for useful discussions. JK acknowledges funding for the ByoPiC project from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program grant agreement ERC-2015-AdG 695561. We are thankful to the community for developing and maintaining open-source software packages extensively used in our work, namely CYTHON (Behnel et al. 2011), MATPLOTLIB (Hunter 2007) and NUMPY (Harris et al. 2020).

## References

Abazajian, K. N., Adshead, P., Ahmed, Z., et al. 2016, ArXiv e-prints [arXiv:1610.02743]

Ade, P., Aguirre, J., Ahmed, Z., et al. 2019, J. Cosmol. Astropart. Phys., 2019, 056

Agarwal, N., Desjacques, V., Jeong, D., & Schmidt, F. 2021, J. Cosmol. Astropart. Phys., 2021, 021

Aker, M., Beglarian, A., Behrens, J., et al. 2021, ArXiv e-prints [arXiv:2105.08533]

Alam, S., Aubert, M., Avila, S., et al. 2021, Phys. Rev. D, 103, 083533

Anderson, T. 2003, An Introduction to Multivariate Statistical Analysis, Wiley Series in Probability and Statistics (New York: Wiley)

Battaglia, N. 2016, J. Cosmol. Astropart. Phys., 2016, 058

Bayer, A. E., Villaescusa-Navarro, F., Massara, E., et al. 2021, ApJ, 919, 24

Behnel, S., Bradshaw, R., Citro, C., et al. 2011, Comput. Sci. Eng., 13, 13

Calafut, V., Gallardo, P. A., Vavagiakis, E. M., et al. 2021, Phys. Rev. D, 104, 043502

Capozzi, F., Di Valentino, E., Lisi, E., et al. 2017, Phys. Rev. D, 95, 096014

Castorina, E., Carbone, C., Bel, J., Sefusatti, E., & Dolag, K. 2015, J. Cosmol. Astropart. Phys., 2015, 043

Chaves-Montero, J., Hernández-Monteagudo, C., Angulo, R. E., & Emberson, J. D. 2021, MNRAS, 503, 1798

Chen, Z., Zhang, P., Yang, X., & Zheng, Y. 2022, MNRAS, 510, 5916

Choudhury, S. R., & Hannestad, S. 2020, J. Cosmol. Astropart. Phys., 2020, 037

Chudaykin, A., & Ivanov, M. M. 2019, J. Cosmol. Astropart. Phys., 2019, 034

De Bernardis, F., Aiola, S., Vavagiakis, E. M., et al. 2017, J. Cosmol. Astropart. Phys., 3, 008

de Salas, P. F., Forero, D. V., Ternes, C. A., Tortola, M., & Valle, J. W. F. 2017, ArXiv e-prints [arXiv:1708.01186]

DESI Collaboration (Aghamousa, A., et al.) 2016, ArXiv e-prints [arXiv:1611.00036]

Ferreira, P. G., Juszkiewicz, R., Feldman, H. A., Davis, M., & Jaffe, A. H. 1999, ApJ, 515, L1

Fisher, K. B. 1995, ApJ, 448, 494

Flender, S., Nagai, D., & McDonald, M. 2017, ApJ, 837, 124

Forero, D. V., Tórtola, M., & Valle, J. W. F. 2014, Phys. Rev. D, 90, 093006

García-Farieta, J. E., Marulli, F., Veropalumbo, A., et al. 2019, MNRAS, 488, 1987

Gonzalez-Garcia, M. C., Maltoni, M., & Schwetz, T. 2016, Nucl. Phys. B, 908, 199

Gualdi, D., & Verde, L. 2020, J. Cosmol. Astropart. Phys., 2020, 041

Hahn, C., & Villaescusa-Navarro, F. 2021, J. Cosmol. Astropart. Phys., 2021, 029

Hahn, C., Villaescusa-Navarro, F., Castorina, E., & Scoccimarro, R. 2020, J. Cosmology Astropart. Phys., 2020, 040

Hand, N., Addison, G. E., Aubourg, E., et al. 2012, Phys. Rev. Lett., 109, 041101

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 585

Hartlap, J., Simon, P., & Schneider, P. 2007, A&A, 464, 399

Heavens, A. 2009, ArXiv e-prints [arXiv:0906.0664]

Hernández-Monteagudo, C., Ma, Y.-Z., Kitaura, F. S., et al. 2015, Phys. Rev. Lett., 115, 191301

Hunter, J. D. 2007, Comput. Sci. Eng., 9, 9

Juszkiewicz, R., Fisher, K. B., & Szapudi, I. 1998, ApJ, 504, L1

Kaufmann, G. M. 1967, Some Bayesian Moment Formulae, Report No. 6710. Centre for Operations Research and Econometrics (Heverlee: Catholic University of Louvain)

Keisler, R., & Schmidt, F. 2013, ApJ, 765, L32

Kodwani, D., Alonso, D., & Ferreira, P. 2019, Open J. Astrophys., 2, 3

Kuruvilla, J., & Aghanim, N. 2021, A&A, 653, A130

Kuruvilla, J., & Porciani, C. 2020, J. Cosmol. Astropart. Phys., 2020, 043

Kuruvilla, J., Aghanim, N., & McCarthy, I. G. 2020, A&A, 644, A170

Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, ArXiv e-prints [arXiv:1110.3193]

Lesgourgues, J., & Pastor, S. 2012, Adv. High Energy Phys., 2012, 608515

Li, Y.-C., Ma, Y.-Z., Remazeilles, M., & Moodley, K. 2018, Phys. Rev. D, 97, 023514

Madhavacheril, M. S., Battaglia, N., Smith, K. M., & Sievers, J. L. 2019, Phys. Rev. D, 100, 103532

Mueller, E.-M., de Bernardis, F., Bean, R., & Niemack, M. D. 2015, Phys. Rev. D, 92, 063501

Nguyen, N.-M., Jasche, J., Lavaux, G., & Schmidt, F. 2020, J. Cosmol. Astropart. Phys., 2020, 011

Planck Collaboration VI. 2020, A&A, 641, A6

Planck Collaboration Int. XXXVII. 2016, A&A, 586, A140

Reid, B. A., & White, M. 2011, MNRAS, 417, 1913

Saito, S., Takada, M., & Taruya, A. 2008, Phys. Rev. Lett., 100, 191301

Samushia, L., Slepian, Z., & Villaescusa-Navarro, F. 2021, MNRAS, 505, 628

Schaan, E., Ferraro, S., Vargas-Magaña, M., et al. 2016, Phys. Rev. D, 93, 082002

Schaan, E., Ferraro, S., Amodeo, S., et al. 2021, Phys. Rev. D, 103, 063513

Smith, K. M., Madhavacheril, M. S., Münchmeyer, M., et al. 2018, ArXiv e-prints [arXiv:1810.13423]

Soergel, B., Flender, S., Story, K. T., et al. 2016, MNRAS, 461, 3172

Soergel, B., Saro, A., Giannantonio, T., Efstathiou, G., & Dolag, K. 2018, MNRAS, 478, 5320

Spergel, D., Gehrels, N., Baltay, C., et al. 2015, ArXiv e-prints [arXiv:1503.03757]

Springel, V. 2005, MNRAS, 364, 1105

Sunyaev, R. A., & Zeldovich, I. B. 1980, MNRAS, 190, 413

Sunyaev, R. A., & Zeldovich, Y. B. 1972, Comments Astrophys. Space Phys., 4, 173

Tanimura, H., Zaroubi, S., & Aghanim, N. 2021, A&A, 645, A112

Tegmark, M., Taylor, A. N., & Heavens, A. F. 1997, ApJ, 480, 22

Uhlemann, C., Friedrich, O., Villaescusa-Navarro, F., Banerjee, A., & Codis, S. 2020, MNRAS, 495, 4006

Vagnozzi, S., Dhawan, S., Gerbino, M., et al. 2018, Phys. Rev. D, 98, 083501

Verde, L. 2010, Statistical Methods in Cosmology (Berlin: Springer Verlag), 800, 147

Villaescusa-Navarro, F., Banerjee, A., Dalal, N., et al. 2018, ApJ, 861, 53

Villaescusa-Navarro, F., Hahn, C., Massara, E., et al. 2020, ApJS, 250, 2

Wong, Y. Y. Y. 2008, J. Cosmol. Astropart. Phys., 2008, 035

Wong, Y. Y. Y. 2011, Ann. Rev. Nucl. Part. Sci., 61, 69

Yankelevich, V., & Porciani, C. 2019, MNRAS, 483, 2078

Yasini, S., Mirzatuny, N., & Pierpaoli, E. 2019, ApJ, 873, L23
**Appendix A: Shape dependence of R**

<figure><img src="image_8.png" /><figcaption>**Fig. A.1.** Comparison of $\mathcal{R}$ measured directly from the halo catalogs (left panel), and the LO analytical prediction (right panel). The triangular configurations corresponds to the set where one of the leg of the triangle ($r_{12}$) is fixed to (115, 120) $h^{-1}$ Mpc.</figcaption></figure>

We compare the direct measurement of $\mathcal{R}$ from the halo catalogs (left panel) and the theoretical predictions (right panel) for a

<figure><img src="image_9.png" /><figcaption>**Fig. A.2.** Ratio between $\mathcal{R}$ measured directly from the halo catalogs, and the LO analytical prediction. The triangular configurations are the same as in Fig. A.1.</figcaption></figure>

small subset of triangular configurations in Fig. A.1. In this case, we fix one side of the triangle, namely, $r_{12}$, and vary the other two sides. The top right bin in both panels ($r_{23}/r_{12} = r_{31}/r_{12} = 1$) corresponds to the equilateral configuration. In Fig. A.2, we show the ratio between them and see that the largest discrepancy can be seen at $r_{31}/r_{12} \sim 0.36$ and $r_{23}/r_{12} \sim 0.68$. This is expected as the separation decreases the fidelity of the LO prediction also decreases.