11institutetext: ML Research at Allegro.pl
11email<EMAIL_ADDRESS>
# NeuralNDCG: Direct Optimisation of a Ranking Metric via Differentiable
Relaxation of Sorting

Przemysław Pobrotyn Equal contribution. Radosław Białobrzeski11footnotemark: 1

###### Abstract

Learning to Rank (LTR) algorithms are usually evaluated using Information
Retrieval metrics like Normalised Discounted Cumulative Gain (NDCG) or Mean
Average Precision. As these metrics rely on sorting predicted items’ scores
(and thus, on items’ ranks), their derivatives are either undefined or zero
everywhere. This makes them unsuitable for gradient-based optimisation, which
is the usual method of learning appropriate scoring functions. Commonly used
LTR loss functions are only loosely related to the evaluation metrics, causing
a mismatch between the optimisation objective and the evaluation criterion. In
this paper, we address this mismatch by proposing NeuralNDCG, a novel
differentiable approximation to NDCG. Since NDCG relies on the non-
differentiable sorting operator, we obtain NeuralNDCG by relaxing that
operator using NeuralSort, a differentiable approximation of sorting. As a
result, we obtain a new ranking loss function which is an arbitrarily accurate
approximation to the evaluation metric, thus closing the gap between the
training and the evaluation of LTR models. We introduce two variants of the
proposed loss function. Finally, the empirical evaluation shows that our
proposed method outperforms previous work aimed at direct optimisation of NDCG
and is competitive with the state-of-the-art methods.

###### Keywords:

Learning to Rank ranking metric optimisation NDCG approximation

##  1 Introduction

Ranking is the problem of optimising, conditioned on some context, the
ordering of a set of items in order to maximise a given metric. The metric is
usually an Information Retrieval (IR) criterion chosen to correlate with user
satisfaction. Learning to Rank (LTR) is a machine learning approach to
ranking, concerned with learning the function which optimises the items’ order
from supervised data. In this work, without loss of generality, we assume our
set of items are search results and the context in which we want to optimise
their order is the user query.

Essentially, one would like to learn a function from search results into
permutations. Since the space of all permutations grows factorially in the
size of the search results set, the task of learning such a function directly
becomes intractable. Thus, most common LTR algorithms resort to the approach
known as score & sort. That is, instead of directly learning the correct
permutation of the search results, one learns a scoring function which
predicts relevancies of individual items, in the form of real-valued scores.
Items are then sorted in the descending order of the scores and thus produced
ordering is evaluated using an IR metric of choice. Typically, scoring
functions are implemented as either gradient boosted trees [14] or Multilayer
Perceptrons (MLP) [25]. Recently, there has been work in using the Transformer
[28] architecture as a scoring function [21]. In order to learn a good scoring
function, one needs a tagged dataset of query-search results pairs together
with ground truth relevancy of each search result (in the context of a given
query), as well as a loss function. There has been extensive research into
constructing appropriate loss functions for LTR (see [19] for an overview of
the field). Such loss functions fall into one of three categories: pointwise,
pairwise or listwise. Pointwise approaches treat the problem as a simple
regression or classification of the ground truth relevancy for each individual
search result, foregoing possible interactions between items. In pairwise
approaches, pairs of items are considered as independent variables and the
function is learned to correctly indicate the preference among the pair.
Examples include RankNet [5], LambdaRank [6] or LambdaMART [7]. However, IR
metrics consider entire search results lists at once, unlike pointwise and
pairwise algorithms. This mismatch has motivated listwise approaches, which
compute the loss based on the scores of the entire list of search results. Two
popular listwise approaches are ListNet [8] and ListMLE [29].

What these loss functions have in common is that they are either not connected
or only loosely connected to the IR metrics used in the evaluation. The
performance of LTR models is usually assessed using Normalised Discounted
Cumulative Gain (NDCG) [16] or Mean Average Precision (MAP) [1]. Since such
metrics rely on sorting the ground truth labels according to the scores
predicted by the scoring function, they are either not differentiable or flat
everywhere and thus cannot be used for gradient-based optimisation of the
scoring function. As a result, there is a mismatch between objectives
optimised by the aforementioned pairwise or listwise losses and metrics used
for the evaluation, even though it can be shown that some of such losses
provide upper bounds of IR measures [31], [30]. On the other hand, as
demonstrated in [23], under certain assumptions on the class of the scoring
functions, direct optimisation of IR measures on a large training set is
guaranteed to achieve high test performance on the same IR measure. Thus,
attempts to bridge the gap between LTR optimisation objectives and
discontinuous evaluation metrics are an important research direction.

In this work, we propose a novel approach to directly optimise NDCG by
approximating the sorting operator with NeuralSort [15]. Since the sorting
operator is the source of discontinuity in NDCG (and other IR metrics), by
substituting it with a differentiable approximation we obtain a smooth variant
of the metric.

The main contributions of the paper are:

  * •

We introduce NeuralNDCG, a novel smooth approximation of NDCG based on
differentiable relaxation of the sorting operator. The variants of the
proposed loss are discussed in Section 4.

  * •

We evaluate a Context-Aware Ranker [21] trained with NeuralNDCG loss on Web30K
[22] and Istella [12] datasets. We demonstrate favourable performance of
NeuralNDCG as compared to baselines. In particular, NeuralNDCG outperforms
ApproxNDCG [23], a competing method for direct optimisation of NDCG.

  * •

We provide an open-source Pytorch [20] implementation allowing for the
reproduction of our results available as part of the open-source allRank
framework111https://github.com/allegro/allRank.

The rest of the paper is organised as follows. In Section 2, we review the
related literature. In Section 3 we formalise the problem of LTR. In Section
4, we recap NeuralSort and demonstrate how it can be used to construct a novel
loss function, NeuralNDCG. In Section 5 we discuss our experimental setup and
results. In the final Section 6 we summarise our findings and discuss the
possible future work.

##  2 Related work

As already mentioned in the introduction, most LTR approaches can be
classified into one of three categories: pointwise, pairwise or listwise. For
a comprehensive overview of the field and most common approaches, we refer the
reader to [19].

In this work, we are concerned with the direct optimisation of non-smooth IR
measures. Methods for optimisation of such metrics can be broadly grouped into
two categories. The methods in the first category try to optimise the upper
bounds of IR metrics as surrogate loss functions. Examples include
$\mathrm{SVM}^{\mathrm{map}}$ [31] and $\mathrm{SVM}^{\mathrm{NDCG}}$ [9]
which optimise upper bounds on $1-\mathrm{MAP}$ and $1-\mathrm{NDCG}$,
respectively. On the other hand, ListNet was originally designed to minimise
cross-entropy between predicted and ground truth top-one probability
distributions, and as such its relation to NDCG was ill-understood. Only
recently was it shown to bound NDCG and Mean Reciprocal Rank (MRR) for binary
labels [3]. Further, a modification to ListNet was proposed in [2] for which
it can be shown that it bounds NDCG also for the graded relevance labels.
Popular methods like LambdaRank and LambdaMART forgo explicit formulation of
the loss function and instead heuristically formulate the gradients based on
NDCG considerations. Since the exact loss function is unknown, its theoretical
relation to NDCG is difficult to analyse.

The second category of methods aims to approximate an IR measure with a smooth
function and directly optimise resulting surrogate function. Our method falls
into this category. We propose to smooth-out NDCG by approximating non-
continuous sorting operator used in the computation of that measure. Recent
works proposing continuous approximation to sorting are already mentioned
NeuralSort, SoDeep [13] and smooth sorting as an Optimal Transport problem
[11]. We use NeuralSort for its firm mathematical foundation, the possibility
to control the degree of approximation and ability to generalise beyond the
maximum list length seen in training. SoDeep uses a deep neural network (DNN)
and synthetic data to learn to approximate the sorting operator and as such
lacks the aforementioned properties. Smooth sorting as Optimal Transport
reports similar performance to NeuralSort at benchmark tasks and we aim to
explore the use of it in NeuralNDCG in the future. By replacing the sorting
operator with its continuous approximation, we obtain NeuralNDCG, a
differentiable approximation of the IR measure. Other notable methods for
direct optimisation of NDCG include:

  * •

ApproxNDCG in which authors reformulated NDCG formula to involve summation
over documents, not their ranks. As a result, they introduce a non-
differentiable position function, which they approximate using a sigmoid. This
loss has been recently revisited in a DNN setting in [4].

  * •

SoftRank [27], where authors propose to smooth scores returned by the scoring
function with equal variance Gaussian distributions: thus deterministic scores
become means of Gaussian score distributions. Subsequently, they derive an
$O(n^{3})$ algorithm to compute Rank-Binomial rank distributions using the
smooth scores. Finally, NDCG is approximated by taking its expectation w.r.t.
the rank distribution.

##  3 Preliminaries

In this section, we formalise the problem and introduce the notation used
throughout the paper. Let
$({\bm{x}},{\bm{y}})\in\mathcal{X}^{n}\times\mathbb{Z}_{\geq 0}^{n}$ be a
training example consisting of a vector ${\bm{x}}$ of $n$ items $x_{i}$,
$1\leq i\leq n$, and a vector ${\bm{y}}$ of corresponding non-negative integer
relevance labels. Note that each item $x_{i}$ is itself a $d$-dimensional
vector of numerical features, and should be thought of as representing a
query-document pair. The set $\mathcal{X}$ is the space of all such vectors
${\bm{x}}_{i}$. Thus, a pair $({\bm{x}},{\bm{y}})$ represents a list of
vectorised search results for a given query together with the corresponding
ground truth relevancies. The dataset of all such pairs is denoted $\Psi$. The
goal of LTR is to find a scoring function
$f:\mathcal{X}^{n}\rightarrow\mathbb{R}^{n}$ that maximises the chosen IR
metric on $\Psi$. The scoring function is learned by minimising the empirical
risk
$\mathcal{L}(f)=\frac{1}{|\Psi|}\sum_{({\bm{x}},{\bm{y}})\in\Psi}{\ell({\bm{y}},{\bm{s}})}$
where $\ell(\cdot)$ is a loss function and ${\bm{s}}=f({\bm{x}})$ is the
vector of predicted scores. As discussed earlier, in most LTR approaches there
is a mismatch between the loss function $\ell$ and the evaluation metric,
causing a discrepancy between the learning procedure and its assessment. In
this work, we focus on NDCG as our metric of choice and propose a new loss
called NeuralNDCG, which bridges the gap between the training and the
evaluation. Before we introduce NeuralNDCG, recall the definition of NDCG.

######  Definition 1.

Let $({\bm{x}},{\bm{y}})\in\mathcal{X}^{n}\times\mathbb{Z}_{\geq 0}^{n}$ be a
training example and assume the documents in ${\bm{x}}$ have been ranked in
the descending order of the scores computed using some scoring function $f$.
Let $r_{j}$ denote the relevance of the document ranked at $j$-th position,
$g(\cdot)$ denote a gain function and $d(\cdot)$ denote a discount function.
Then, the Discounted Cumulative Gain at $k$-th position ($k\leq n$) is defined
as

$\mathrm{DCG}@k=\sum_{j=1}^{k}g(r_{j})d(j)$ (1)

and Normalised Discounted Cumulative Gain at $k$ is defined as

$\mathrm{NDCG}@k=\frac{1}{\mathrm{maxDCG}@k}\mathrm{DCG}@k$ (2)

where $\mathrm{maxDCG}@k$ is the maximum possible value of $\mathrm{DCG}@k$,
computed by ordering the documents in ${\bm{x}}$ by their decreasing ground
truth relevancy.

Note that, typically, the discount function $d(\cdot)$ and the gain function
$g(\cdot)$ are given by $d(j)=\frac{1}{\log_{2}(j+1)}$ and
$g(r_{j})=2^{r_{j}}-1$, respectively.

##  4 Loss formulation

In this section we define NeuralNDCG, a novel differentiable approximation to
NDCG. It relies on NeuralSort, a smooth relaxation of the sorting operator. We
begin by recalling NeuralSort, proceed to define NeuralNDCG and discuss its
possible variants.

###  4.1 Sorting relaxation

Recall that sorting a list of scores ${\bm{s}}$ is equivalent to left-
multiplying a column vector of scores by the permutation matrix
$P_{\texttt{sort}({\bm{s}})}$ induced by permutation sort(${\bm{s}}$) sorting
the scores. Thus, in order to approximate the sorting operator, it is enough
to approximate the induced permutation matrix. In [15], the permutation matrix
is approximated via a unimodal row stochastic matrix
$\widehat{P}_{\texttt{sort}({\bm{s}})}(\tau)$ given by:

$\widehat{P}_{\texttt{sort}({\bm{s}})}[i,:](\tau)=\mathrm{softmax}[((n+1-2i){\bm{s}}-A_{{\bm{s}}}\mathds{1})/\tau]$
(3)

where $A_{{\bm{s}}}$ is the matrix of absolute pairwise differences of
elements of ${\bm{s}}$ such that $A_{{\bm{s}}}[i,j]=|s_{i}-s_{j}|$,
$\mathds{1}$ denotes the column vector of all ones and $\tau>0$ is a
temperature parameter controlling the accuracy of approximation. For brevity,
for the remainder of the paper we refer to
$\widehat{P}_{\texttt{sort}({\bm{s}})}(\tau)$ simply as $\widehat{P}$.

Note that the temperature parameter $\tau$ allows to control the trade-off
between the accuracy of the approximation and the variance of the gradients.
Generally speaking, the lower the temperature, the better the approximation at
the cost of a larger variance in the gradients. In fact, it is not difficult
to demonstrate that:

$\lim_{\tau\to
0}\widehat{P}_{\texttt{sort}({\bm{s}})}(\tau)=P_{\texttt{sort}({\bm{s}})}$ (4)

(see [15] for proof). This fact will come in handy once we define NeuralNDCG.

An approximation of a permutation matrix by Equation 3 is a deterministic
function of the predicted scores. Authors of NeuralSort proposed also a
stochastic version, by deriving a reparametrised sampler for a Plackett-Luce
family of distributions. Essentially, they propose to perturb scores $s$ with
a vector $g$ of i.i.d. Gumbel perturbations with zero mean and a fixed scale
$\beta$ to obtain perturbed scores $\tilde{s}=\beta\log s+g$. Perturbed scores
are then used in place of deterministic scores in the formula for
$\widehat{P}$.

We experimented with both deterministic and stochastic approximations to
sorting and found them to yield similar results. Thus, for brevity, in this
work we focus on the deterministic variant.

Table 1: Approximate sorting with NeuralSort. Given ground truth ${\bm{y}}=[4,2,1,0,4,3]$ and predicted scores ${\bm{s}}~{}=~{}[0.5,0.2,0.1,0.01,0.65,0.3]$, ${\bm{y}}$ is sorted by $\widehat{P}$ for different values of $\tau$. Exact sorting is shown in the first row. | Quasi-sorted ground truth | Sum after sorting
---|---|---
$\lim_{\tau\to 0}$ | 4 | 4 | 3 | 2 | 1 | 0 | 14
$\tau=0.01$ | 4 | 4 | 3 | 2 | 0.99992 | 0.00012339 | 14.00004339
$\tau=0.1$ | 3.9995 | 3.8909 | 2.8239 | 1.9730 | 0.9989 | 0.3136 | 13.9998
$\tau=1$ | 3.3893 | 2.9820 | 2.4965 | 2.0191 | 1.6097 | 1.2815 | 10.388

###  4.2 NeuralNDCG

If the ground truth permutation is known, one could minimise the cross-entropy
loss between the ground truth permutation matrix and its approximation given
by $\widehat{P}$, as done in the experiments section in [15]. However, for
many applications, including ranking, the exact ground truth permutation is
not known. Relevance labels of individual items produce many possible valid
ground truth permutation matrices. Thus, instead of optimising the cross-
entropy loss, we use NeuralSort to introduce NeuralNDCG, a novel loss function
appropriate for LTR.

Given a list of documents ${\bm{x}}$, its corresponding vector of scores
${\bm{s}}=f({\bm{x}})$ and the ground truth labels ${\bm{y}}$ we first find
the approximate permutation matrix $\widehat{P}$ induced by the scores
${\bm{s}}$ using Equation 3. We then apply the gain function $g$ to the vector
of ground truths ${\bm{y}}$ and obtain a vector $g({\bm{y}})$ of gains per
document. We then left-multiply the column vector $g({\bm{y}})$ of gains by
$\widehat{P}$ and obtain an “approximately” sorted version of the gains,
$\widehat{g({\bm{y}})}$. Another way to think of that approximate sorting is
that the $k$-th row of $\widehat{P}$ gives weights of documents $x_{i}$ in the
computation of gain at rank $k$ after sorting. Gain at rank $k$ is then the
weighted sum of ground truth gains, weighted by the entries in the $k$-th row
of $\widehat{P}$. Note that after the approximate sorting the actual integer
values of ground truth gains become ”distorted” and are not necessarily
integers anymore (See Table 1 for example). In particular, the sum of quasi-
sorted gains $\widehat{g({\bm{y}})}$ may differ from that of the original
vector $g({\bm{y}})$. This leads to a peculiar behaviour of NeuralNDCG near
the discontinuities of true NDCG (Figure 1), which may be potentially harmful
for optimisation using Stochastic Gradient Descent [24]. Since $\widehat{P}$
is row-stochastic but not-necessarily column-stochastic (i.e. each column does
not necessarily sum to one), an individual ground truth gain $g({\bm{y}})_{j}$
may have corresponding weights in the rows of $\widehat{P}$ that do not sum to
one (and, in particular, may sum to more than one), so it will overcontribute
to the total sum of $\widehat{g({\bm{y}})}$. To alleviate that problem, we
additionally perform Sinkhorn scaling [26] on $\widehat{P}$ (i.e. we
iteratively normalize all rows and columns until convergence222We stop the
procedure after 30 iterations or when the maximum difference between row or
column sum and one is less than $10^{-6}$, whatever happens first.) before
using it for quasi-sorting. This way, the columns also sum to one and the
approximate sorting is smoothed-out (again, see Figure 1). The remaining steps
are identical to the computation of $\mathrm{NDCG}@k$, with the exception that
the gain of the relevance function $r_{j}$ is replaced with the $j$-th
coordinate of quasi-sorted gains $\widehat{g({\bm{y}})}$. For the discount
function $d$, we use the usual inverse logarithmic discount and for the gain
function $g$ we used the usual power function. For the computation of the
maxDCG, we use original ground truth labels ${\bm{y}}$. To find NeuralNDCG at
rank $k$, we simply truncate quasi-sorted gains to the $k$-th position and
compute the maxDCG at $k$-th rank.

We thus obtain the following formula for NeuralNDCG:

$\mathrm{NeuralNDCG}_{k}(\tau)({\bm{s}},{\bm{y}})=N_{k}^{-1}\sum_{j=1}^{k}(\mathrm{scale}(\widehat{P})\cdot
g({\bm{y}}))_{j}\cdot d(j)$ (5)

where $N_{k}^{-1}$ is the maxDCG at $k$-th rank, $\mathrm{scale}(\cdot)$ is
Sinkhorn scaling and $g(\cdot)$ and $d(\cdot)$ are their gain and discount
functions. Note that the summation is over the first $k$ ranks.

Finally, since the popular autograd libraries provide means to minimise a
given loss functions, we use $(-1)\times\mathrm{NeuralNDCG}$ for optimisation.

Figure 1: Given ground truth $y=[2,1,0,0,0]$ and a list of scores
${\bm{s}}=[4,1,0,0,x]$, we vary the value of the score $x$ and plot resulting
NDCG induced by the scores along with NeuralNDCG ($\tau=1.0$) with and without
Sinkhorn scaling of $\widehat{P}$.

###  4.3 NeuralNDCG Transposed

In the above formulation of NeuralNDCG, the summation is done over the ranks
and gain at each rank is computed as a weighted sum of all gains, with weights
given by the rows of $\widehat{P}$. We now provide an alternative formulation
of NeuralNDCG, called NeuralNDCG Transposed ($\mathrm{NeuralNDCG}^{T}$ for
short), where the summation is done over the documents, not their ranks.

As previously, let ${\bm{x}}$ be a list of documents with corresponding scores
${\bm{s}}$ and ground truth relevancies ${\bm{y}}$. We begin by finding the
approximate permutation matrix $\widehat{P}$. Since we want to sum over the
documents and not their ranks, we need to find the weighted average of
discounts per document, not the weighted average of gains per rank as before.
To this end, we transpose $\widehat{P}$ to obtain an approximation
$\widehat{P}^{T}$ of the inverse of the permutation matrix corresponding to
sorting the documents ${\bm{x}}$ by their corresponding scores ${\bm{y}}$.
Thus, $\widehat{P}^{T}$ can be thought of as an approximate unsorting matrix -
when applied to sorted documents (ranks), it will (approximately) recover
their original ordering. Since $\widehat{P}$ is row-stochastic,
$\widehat{P}^{T}$ is column-stochastic. As we want to apply it by left-
multiplication, we want it to be row-stochastic. Thus, similarly to before, we
perform Sinkhorn scaling of $\widehat{P}^{T}$. After Sinkhorn scaling, the
$k$-th row of $\widehat{P}^{T}$ can be thought of as giving the weights of
different ranks when computing the discount of the $k$-th document. We can now
find the vector of the weighted averages of discounts per document by
computing $\widehat{P}^{T}{\bm{d}}$, where ${\bm{d}}$ is the vector of
logarithmic discounts per rank (${\bm{d}}_{j}=d(j)$). Note that since we want
to perform summation over the documents, not ranks, it is not enough to sum
the first $k$ elements to truncate NDCG to the $k$-th position. Instead, the
entries of the discounts vector ${\bm{d}}$ corresponding to ranks $j>k$ are
set to 0. This way, the documents which would end up at ranks $j>k$ after
sorting end up having weighted discounts being close to 0, and equal to 0 in
the limit of the temperature $\tau$. Thus, even though the summation is done
over all documents, we still recover $\mathrm{NDCG}@k$.

Hence, $\mathrm{NeuralNDCG}^{T}$ is given by the following formula:

${\mathrm{NeuralNDCG}^{T}}_{k}(\tau)({\bm{s}},{\bm{y}})=N_{k}^{-1}\sum_{i=1}^{n}g({\bm{y}}_{i})\cdot(\mathrm{scale}(\widehat{P}^{T})\cdot{\bm{d}})_{i}$
(6)

where $N_{k}^{-1}$ is the maxDCG at $k$-th rank, $\mathrm{scale}(\cdot)$ is
Sinkhorn scaling, $g(\cdot)$ is the gain function, ${\bm{d}}$ is the vector of
logarithmic discounts per rank set to 0 for ranks $j>k$, and the summation is
done over all $n$ documents.

###  4.4 Properties of NeuralNDCG

By Equation 4, in the limit of the temperature, the approximate permutation
matrix $\widehat{P}$ becomes the true permutation matrix $P$. Thus, as the
temperature approaches zero, NeuralNDCG approaches true NDCG in both its
variants. See Figure 2 for examples of the effect of the temperature on the
accuracy of the approximation.

Comparing to ApproxNDCG, our proposed approximation to NDCG showcases more
favourable properties. We can easily compute NDCG at any rank position $k$,
whereas in ApproxNDCG, one needs to further approximate the truncation
function using an approximation of the position function. This approximation
of an approximation leads to a compounding of errors. We deal away with that
problem by using a single approximation of the permutation matrix.
Furthermore, the approximation of the position function in ApproxNDCG is done
using a sigmoid function, which may lead to the vanishing gradient problem.

SoftRank suffers from a high computational complexity of $O(n^{3})$: in order
to compute all the derivatives required by the algorithm, a recursive
computation is necessary. Authors relieve that cost by approximating all but a
few of the Rank-Binomial distributions used, but at a cost of the accuracy of
their solution. On the other hand, computation of $\widehat{P}$ is of
$O(n^{2})$ complexity.

Figure 2: Given ground truth ${\bm{y}}=[1,2,3,4,5]$ and a list of scores
${\bm{s}}=[1,2,3,4,x]$, we vary the value of the score $x$ and plot resulting
NDCG induced by the scores along with Deterministic NeuralNDCG for different
temperatures $\tau$.

##  5 Experiments

This section describes the experimental setup used to empirically verify the
proposed loss functions.

###  5.1 Datasets

Experiments were conducted on two datasets: Web30K and Istella 333There are a
few variants of this dataset, we used the Istella full dataset.. Both datasets
consists of queries and associated search results. Each query-document pair is
represented by a real-valued feature vector and has an associated graded
relevance on the scale from 0 (irrelevant) to 4 (highly relevant). For both
datasets, we standardise the features, log-transforming selected ones, before
feeding them to the learning algorithm. Since the lengths of search results
lists in the datasets are unequal, we padded or subsampled to equal length for
training, but used the full list length for evaluation. Web30K comes split
into five folds. However, following the common practice in the field, we
report results obtained on Fold 1 of the data. We used 60% of the data for
training, 20% for validation and hyperparameter tuning and the remaining 20%
for testing. Istella datasets comes partition into a training and a test fold
according to a 80%-20% schema. We additionally split the training data into
training and validation data to obtain a 60%/20%/20% split, similarly to
Web30K. We tune the hyperparameters of our models on the validation data and
report performance on the test set, having trained the best models on the full
training fold. In both datasets there are a number of queries for which the
associated search results list contains no relevant documents (i.e. all
documents have label 0). We refer to these queries as empty queries. For such
queries, the NDCG of their list of results can be arbitrarily set to either 0
or 1. To allow for a fair comparison with the current state of the art, we
followed LightGBM [17] implementation of setting NDCG of such lists to 1
during the evaluation. Table 2 summaries the characteristics of the datasets
used.

Table 2: Dataset statistics Dataset | Features | Queries in training | Queries in test | Empty queries
---|---|---|---|---
Web30K | 136 | 18919 | 6306 | 982
Istella | 220 | 23219 | 9799 | 50

###  5.2 Scoring function

For the scoring function $f$, we used the Context-Aware Ranker, a ranking
model based on the Transformer architecture. The model can be thought of as
the encoder part of the Transformer, taking raw features of items present in
the same list as input and outputting a real-valued score for each item. Given
the ubiquity of Transformer-based models in the literature, we refer to reader
to [21] for the details of the architecture used. Compared to the original
network described in [21], we used smaller architectures. For both datasets,
we used an architecture consisting of 2 encoder blocks of a single attention
head each, with a hidden dimension of 384. The dimensionality of initial
fully-connected layer was set to 96 for models trained on Web30K and 128 for
models trained on Istella. We did not apply any activation on the output
except for NeuralNDCG and NeuralNDCGT. It exhibited suboptimal performance
without any nonlinear output activation function and, in this case, we applied
Tanh to the output. For both datasets, the same architectures were used across
all loss functions.

###  5.3 Training hyperparameters

In all cases, we used Adam [18] optimiser and set the learning rate to
$0.001$. The batch size was set to 64 (Web30K) or 110 (Istella) and search
results lists were truncated or padded to the length of 240 when training. We
trained the networks for 100 epochs, decaying the learning rate by the factor
of $0.1$ after 50 epochs.

###  5.4 Loss functions

We compared the performance of variants of NeuralNDCG against the following
loss functions. For a pointwise baseline, we used a simple RMSE of predicted
relevancy. Specifically, the output of the network $f$ is passed through a
sigmoid function and then multiplied by the number of relevance levels. The
root mean squared difference of this score and the ground truth relevance is
the loss. Pairwise losses we compared with consist of RankNet and LambdaRank.
Similarly to NeuralNDCG, these losses support training with a specific rank
cutoff. We thus train models with these losses at ranks $5$, $10$ and at the
maximum rank. The two most popular listwise losses are ListNet and ListMLE,
and we, too, included them in our evaluation. Finally, the other method of
direct optimisation of NDCG which we compared with was ApproxNDCG. We did not
compare with SoftRank, as its $O(n^{3})$ complexity proved prohibitive. We
tuned ApproxNDCG and NeuralNDCG smoothness hyperparameters for optimal
performance on the test set. Both ApproxNDCG’s $\alpha$ and NeuralNDCG’s
$\tau$ parameter were set to 1 as other values in the $[0.01;100]$ interval
did not show any improvement.

Table 3: Test NDCG on Web30K and Istella. Boldface is the best performing loss column-wise. Loss | WEB30K | Istella
---|---|---
| NDCG@5 | NDCG@10 | NDCG@5 | NDCG@10
NeuralNDCG@5 | 50.32 | 52.01 | 65.32 | 69.97
NeuralNDCG@10 | 50.89 | 52.77 | 65.65 | 70.68
NeuralNDCG@max | 51.56 | 53.46 | 65.69 | 70.55
$\mathrm{NeuralNDCG}^{T}@5$ | 50.50 | 52.14 | 65.46 | 69.95
$\mathrm{NeuralNDCG}^{T}@10$ | 50.85 | 52.70 | 66.02 | 71.02
$\mathrm{NeuralNDCG}^{T}$@max | 51.45 | 53.49 | 65.60 | 70.53
ApproxNDCG | 49.07 | 50.90 | 63.14 | 67.94
ListNet | 50.75 | 52.80 | 65.62 | 70.70
ListMLE | 49.81 | 51.82 | 59.85 | 66.24
RankNet@5 | 49.14 | 50.75 | 64.45 | 68.74
RankNet@10 | 50.95 | 52.69 | 65.75 | 70.68
RankNet@max | 49.84 | 51.82 | 64.57 | 70.37
LambdaRank@5 | 48.70 | 50.10 | 63.50 | 67.75
LambdaRank@10 | 49.66 | 51.34 | 65.21 | 69.82
LambdaRank@max | 51.55 | 53.47 | 65.90 | 71.09
RMSE | 50.51 | 52.46 | 65.62 | 70.76
XGBoost | 46.80 | 49.17 | 61.04 | 65.74

###  5.5 Results

For both datasets, we report models performance in terms of NDCG@5 and
NDCG@10. Results are collected in Table 3. Both NeuralNDCG variants in every
rank cutoff setting outperform ApproxNDCG on both datasets in all metrics
reported. Moreover, NeuralNDCG variants with specific rank cutoffs provide the
best performance among all losses in both metrics on the WEB30K dataset and
NDCG@5 on the Istella dataset. For reference, we also report the results of a
GBDT model trained with XGBoost [10] and objective rank:pairwise (as rank:ndcg
is known to yield suboptimal results444For details, please visit
https://github.com/dmlc/xgboost/issues/6352.).

##  6 Conclusions

In this work we introduced NeuralNDCG, a novel differentiable approximation of
NDCG. By substituting the discontinuous sorting operator with NeuralSort, we
obtain a robust, efficient and arbitrarily accurate approximation to NDCG. Not
only does it enjoy favourable theoretical properties, but also proves to be
effective in empirical evaluation, yielding competitive performance, on par
with LambdaRank. This work can easily be extended to other rank-based metrics
like MAP; a possibility we aim to explore in the future. Another interesting
extension of this work would be the substitution of NeuralSort with another
method of approximation of the sorting operator, most notably the method
treating sorting as an Optimal Transport problem [11].

## References

  * [1] Baeza-Yates, R.A., Ribeiro-Neto, B.: Modern Information Retrieval. Addison-Wesley Longman Publishing Co., Inc., Boston, MA, USA (1999)
  * [2] Bruch, S.: An alternative cross entropy loss for learning-to-rank (2019)
  * [3] Bruch, S., Wang, X., Bendersky, M., Najork, M.: An analysis of the softmax cross entropy loss for learning-to-rank with binary relevance. In: Proceedings of the 2019 ACM SIGIR International Conference on Theory of Information Retrieval. p. 75–78. ICTIR ’19, Association for Computing Machinery, New York, NY, USA (2019). https://doi.org/10.1145/3341981.3344221, https://doi.org/10.1145/3341981.3344221
  * [4] Bruch, S., Zoghi, M., Bendersky, M., Najork, M.: Revisiting approximate metric optimization in the age of deep neural networks. In: Proceedings of the 42nd International ACM SIGIR Conference on Research and Development in Information Retrieval. p. 1241–1244. SIGIR’19, Association for Computing Machinery, New York, NY, USA (2019). https://doi.org/10.1145/3331184.3331347, https://doi.org/10.1145/3331184.3331347
  * [5] Burges, C., Shaked, T., Renshaw, E., Lazier, A., Deeds, M., Hamilton, N., Hullender, G.: Learning to rank using gradient descent. In: Proceedings of the 22Nd International Conference on Machine Learning. pp. 89–96. ICML ’05, ACM, New York, NY, USA (2005). https://doi.org/10.1145/1102351.1102363, http://doi.acm.org/10.1145/1102351.1102363
  * [6] Burges, C.J., Ragno, R., Le, Q.V.: Learning to rank with nonsmooth cost functions. In: Schölkopf, B., Platt, J.C., Hoffman, T. (eds.) Advances in Neural Information Processing Systems 19, pp. 193–200. MIT Press (2007), http://papers.nips.cc/paper/2971-learning-to-rank-with-nonsmooth-cost-functions.pdf
  * [7] Burges, C.J.C.: From RankNet to LambdaRank to LambdaMART: An overview. Tech. rep., Microsoft Research (2010), http://research.microsoft.com/en-us/um/people/cburges/tech_reports/MSR-TR-2010-82.pdf
  * [8] Cao, Z., Qin, T., Liu, T.Y., Tsai, M.F., Li, H.: Learning to rank: From pairwise approach to listwise approach. In: Proceedings of the 24th International Conference on Machine Learning. pp. 129–136. ICML ’07, ACM, New York, NY, USA (2007). https://doi.org/10.1145/1273496.1273513, http://doi.acm.org/10.1145/1273496.1273513
  * [9] Chakrabarti, S., Khanna, R., Sawant, U., Bhattacharyya, C.: Structured learning for non-smooth ranking losses. In: Proceedings of the 14th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. p. 88–96. KDD ’08, Association for Computing Machinery, New York, NY, USA (2008). https://doi.org/10.1145/1401890.1401906, https://doi.org/10.1145/1401890.1401906
  * [10] Chen, T., Guestrin, C.: Xgboost: A scalable tree boosting system. In: Proceedings of the 22Nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. pp. 785–794. KDD ’16, ACM, New York, NY, USA (2016). https://doi.org/10.1145/2939672.2939785, http://doi.acm.org/10.1145/2939672.2939785
  * [11] Cuturi, M., Teboul, O., Vert, J.P.: Differentiable ranking and sorting using optimal transport. In: Wallach, H., Larochelle, H., Beygelzimer, A., d Alche-Buc, F., Fox, E., Garnett, R. (eds.) Advances in Neural Information Processing Systems 32, pp. 6858–6868. Curran Associates, Inc. (2019), http://papers.nips.cc/paper/8910-differentiable-ranking-and-sorting-using-optimal-transport.pdf
  * [12] Dato, D., Lucchese, C., Nardini, F.M., Orlando, S., Perego, R., Tonellotto, N., Venturini, R.: Fast ranking with additive ensembles of oblivious and non-oblivious regression trees. ACM Trans. Inf. Syst. 35(2) (Dec 2016). https://doi.org/10.1145/2987380, https://doi.org/10.1145/2987380
  * [13] Engilberge, M., Chevallier, L., Perez, P., Cord, M.: Sodeep: A sorting deep net to learn ranking loss surrogates. In: The IEEE Conference on Computer Vision and Pattern Recognition (CVPR) (June 2019)
  * [14] Friedman, J.H.: Greedy function approximation: A gradient boosting machine. Annals of Statistics 29, 1189–1232 (2000)
  * [15] Grover, A., Wang, E., Zweig, A., Ermon, S.: Stochastic optimization of sorting networks via continuous relaxations. In: International Conference on Learning Representations (2019), https://openreview.net/forum?id=H1eSS3CcKX
  * [16] Järvelin, K., Kekäläinen, J.: Cumulated gain-based evaluation of ir techniques. ACM Trans. Inf. Syst. 20(4), 422–446 (Oct 2002). https://doi.org/10.1145/582415.582418, http://doi.acm.org/10.1145/582415.582418
  * [17] Ke, G., Meng, Q., Finley, T., Wang, T., Chen, W., Ma, W., Ye, Q., Liu, T.Y.: Lightgbm: A highly efficient gradient boosting decision tree. In: Guyon, I., Luxburg, U.V., Bengio, S., Wallach, H., Fergus, R., Vishwanathan, S., Garnett, R. (eds.) Advances in Neural Information Processing Systems 30, pp. 3146–3154. Curran Associates, Inc. (2017), http://papers.nips.cc/paper/6907-lightgbm-a-highly-efficient-gradient-boosting-decision-tree.pdf
  * [18] Kingma, D.P., Ba, J.: Adam: A method for stochastic optimization (2014), http://arxiv.org/abs/1412.6980, cite arxiv:1412.6980Comment: Published as a conference paper at the 3rd International Conference for Learning Representations, San Diego, 2015
  * [19] Liu, T.Y.: Learning to rank for information retrieval. Found. Trends Inf. Retr. 3(3), 225–331 (Mar 2009). https://doi.org/10.1561/1500000016, http://dx.doi.org/10.1561/1500000016
  * [20] Paszke, A., Gross, S., Chintala, S., Chanan, G., Yang, E., DeVito, Z., Lin, Z., Desmaison, A., Antiga, L., Lerer, A.: Automatic differentiation in PyTorch. In: NIPS Autodiff Workshop (2017)
  * [21] Pobrotyn, P., Bartczak, T., Synowiec, Mikołaj an Białobrzeski, R., Bojar, J.: Context-aware learning to rank with self-attention. In: SIGIR eCom ’20. Virtual Event, China. (2020)
  * [22] Qin, T., Liu, T.M.: Introducing letor 4.0 datasets. ArXiv abs/1306.2597 (2013)
  * [23] Qin, T., Liu, T.Y., Li, H.: A general approximation framework for direct optimization of information retrieval measures. Inf. Retr. 13, 375–397 (08 2010). https://doi.org/10.1007/s10791-009-9124-x
  * [24] Robbins, H., Monro, S.: A stochastic approximation method. Annals of Mathematical Statistics 22, 400–407 (1951)
  * [25] Rumelhart, D.E., Hinton, G.E., Williams, R.J.: Learning internal representations by error propagation. Tech. rep., California Univ San Diego La Jolla Inst for Cognitive Science (1985)
  * [26] Sinkhorn, R.: A relationship between arbitrary positive matrices and doubly stochastic matrices. Ann. Math. Statist. 35(2), 876–879 (06 1964). https://doi.org/10.1214/aoms/1177703591, https://doi.org/10.1214/aoms/1177703591
  * [27] Taylor, M., Guiver, J., Robertson, S., Minka, T.: Softrank: optimizing non-smooth rank metrics. pp. 77–86 (01 2008). https://doi.org/10.1145/1341531.1341544
  * [28] Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A.N., Kaiser, L.u., Polosukhin, I.: Attention is all you need. In: Guyon, I., Luxburg, U.V., Bengio, S., Wallach, H., Fergus, R., Vishwanathan, S., Garnett, R. (eds.) Advances in Neural Information Processing Systems 30, pp. 5998–6008. Curran Associates, Inc. (2017), http://papers.nips.cc/paper/7181-attention-is-all-you-need.pdf
  * [29] Xia, F., Liu, T.Y., Wang, J., Zhang, W., Li, H.: Listwise approach to learning to rank: Theory and algorithm. In: Proceedings of the 25th International Conference on Machine Learning. pp. 1192–1199. ICML ’08, ACM, New York, NY, USA (2008). https://doi.org/10.1145/1390156.1390306, http://doi.acm.org/10.1145/1390156.1390306
  * [30] Xu, J., Li, H.: Adarank: a boosting algorithm for information retrieval. In: Proceedings of the 30th annual ACM SIGIR conference. pp. 391–398. ACM, Amsterdam, The Netherlands (2007). https://doi.org/10.1145/1277741.1277809
  * [31] Yue, Y., Finley, T., Radlinski, F., Joachims, T.: A support vector method for optimizing average precision. In: Proceedings of the 30th Annual International ACM SIGIR Conference on Research and Development in Information Retrieval. p. 271–278. SIGIR ’07, Association for Computing Machinery, New York, NY, USA (2007). https://doi.org/10.1145/1277741.1277790, https://doi.org/10.1145/1277741.1277790


11institutetext: Université Paris-Saclay, CNRS, Institut d’Astrophysique
Spatiale, 91405, Orsay, France

# Cosmology with the kinetic Sunyaev–Zeldovich effect: Independent of the
optical depth and $\sigma_{8}$

Joseph Kuruvilla

(Received 2021; accepted)

The cosmological constraints from the kinetic Sunyaev–Zeldovich experiments
are degenerate with the optical depth measurement, which 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 independent of the optical
depth, and $\sigma_{8}$. Using 22,000 $N$-body simulations from the Quijote
suite, we quantify 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 have a factor of 6.2–12.9, and 2.3–5.7 improvement on all
cosmological model parameters when compared to those obtained from the mean
pairwise velocity, and the redshift-space halo power spectrum, respectively.
Thus the new statistic paves a way forward to constrain cosmological
parameters independent of the optical depth and $\sigma_{8}$ using data from
future kinetic Sunyaev–Zeldovich experiments alone.

###### Key Words.:

large-scale structure of the Universe – cosmology: theory – cosmic background
radiation

††offprints<EMAIL_ADDRESS>
##  1 Introduction

In modern physics, one of the outstanding question is regarding the
determination of the mass of neutrinos which has fundamental implications to
both particle physics and cosmology. The neutrino oscillation experiments have
established that the neutrinos should 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 to measure the
absolute scale of the neutrino mass other experiments are required. Recently
the Karlsruhe Tritium Neutrino (KATRIN) experiment has 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 $\beta$-decay, and is model independent.

However stronger constraints on the summed neutrino mass ($M_{\nu}$) can
obtained through combining various cosmological probes, as the massive
neutrinos leave an imprint on various cosmological observables (e.g. Wong,
2011; Lesgourgues & Pastor, 2012). But these constraints have an additional
model dependence. Assuming the standard ‘lambda cold dark matter’
($\Lambda$CDM) model, one of the strongest constraint on the summed neutrino
mass has been obtained by combining 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_{\nu}<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 Observatory111https://simonsobservatory.org/ (SO, Ade et al.,
2019), CMB-S4222https://cmb-s4.org/ (Abazajian et al., 2016)], and the stage
IV galaxy redshift surveys [e.g. the ‘Dark Energy Spectroscopic
Instrument’333https://www.desi.lbl.gov/ (DESI, DESI Collaboration et al.,
2016), Euclid444https://www.euclid-ec.org/ (Laureijs et al., 2011) and Nancy
Grace Roman space telescope555https://roman.gsfc.nasa.gov/ (Spergel et al.,
2015)].

Currently the $M_{\nu}$ constraints from galaxy clustering are mainly obtained
using two-point statistics, i.e. the power spectrum in 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 real- (e.g. Saito et al.,
2008; Wong, 2008; Castorina et al., 2015) and redshift-space (e.g.
Villaescusa-Navarro et al., 2018; García-Farieta et al., 2019). However they
are affected by the $M_{\nu}$-$\sigma_{8}$ degeneracy, and thus acts as a
limitation in measuring the summed neutrino mass. The three-point clustering
statistics in 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 contains additional
cosmological information compared to its two-point counterpart, and thus able
to obtain substantial improvements on constraining other cosmological
parameters also (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, like 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 avenue is to use velocity statistics like 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) has shown that the mean pairwise velocity can
utilised to constrain the summed neutrino mass, and Kuruvilla et al. (2020)
has 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 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 will develop a new
statistic which is independent of the optical depth using the first moment of
the three-point relative velocities, i.e. the mean relative velocities between
pairs in a triplet, and thus circumvent the problem of optical depth as a
limitation factor in kSZ experiments. This forms 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, which we use 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, and we finally
conclude 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 which is
known as the kinetic Sunyaev–Zeldovich (kSZ; Sunyaev & Zeldovich, 1972, 1980)
which arises if the scattering medium is moving relative to the Hubble flow.
The fractional temperature fluctuation caused due to kSZ is

$\displaystyle\left.\frac{\Delta
T(\hat{n})}{T_{\mathrm{cmb}}}\right|_{\mathrm{kSZ}}$
$\displaystyle=-\int\mathrm{d}l\
\sigma_{\mathrm{T}}\left(\frac{\bm{v}_{\mathrm{e}}\cdot\hat{n}}{c}\right)n_{\mathrm{e}}\,,$
$\displaystyle=-\tau\left(\frac{\bm{v}_{\mathrm{e}}\cdot\hat{n}}{c}\right)\,,$
(1)

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

The kSZ signal detection is challenging because of its small amplitude and its
spectrum being identical to that of primary CMB temperature fluctuations. One
of the approaches 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 et al., 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.,
2021). There have been evidences of kSZ signal using other techniques also
(e.g. Hernández-Monteagudo et al., 2015; Schaan et al., 2016; Nguyen et al.,
2020; Tanimura et al., 2021; Chaves-Montero et al., 2021; Schaan et al.,
2021).

In the case of 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^{\mathrm{kSZ}}(r_{12})\rangle}{T_{\mathrm{cmb}}}\simeq-\tau\frac{\bar{w}(r_{12})}{c}\,,$
(2)

where $\langle\Delta T^{\mathrm{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 which can be defined in the single
streaming regime as follows

$\displaystyle\langle\bm{w}_{12}|\bm{r}_{12}\rangle_{\mathrm{p}}=$
$\displaystyle\
\displaystyle\frac{\langle(1+\delta_{1})(1+\delta_{2})(\bm{v}_{2}-\bm{v}_{1})\rangle}{\langle(1+\delta_{1})(1+\delta_{2})\rangle}\,,$
(3)

where $\delta_{i}\equiv\delta(\bm{x}_{i})$ is the density contrast,
$\bm{v}_{i}\equiv\bm{v}(\bm{x}_{i})\equiv\bm{u}(\bm{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\bm{w}_{12}|\bm{r}_{12}\rangle_{\mathrm{p}}=\bar{w}(r_{12})\,\hat{\bm{r}}_{12}\simeq\displaystyle-\frac{f}{\pi^{2}}\,\hat{\bm{r}}_{12}\int_{0}^{\infty}k\,j_{1}(kr_{12})\,P(k)\,\mathrm{d}k\,,$
(4)

where $\hat{\bm{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^{\mathrm{kSZ}}\propto\tau f\sigma^{2}_{8}\,$ (5)

and thus implying that the growth rate measurement from the pairwise kSZ is
perfectly degenerate with optical depth (Keisler & Schmidt, 2013). Here we
have presented the argument for the matter, but in observations there will be
an additional bias dependence entering in the above equation.

###  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
$\triangle_{123}=(r_{12},r_{23},r_{31})$: (i)
$\langle\bm{w}_{12}|\triangle_{123}\rangle_{\mathrm{t}}$, and (ii)
$\langle\bm{w}_{23}|\triangle_{123}\rangle_{\mathrm{t}}$. The subscript t here
implies that the averages are computed over all triplets with separations
$(r_{12},r_{23},r_{31})$. Similar 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)

$\displaystyle\langle\bm{w}_{12}|\triangle_{123}\rangle_{\mathrm{t}}=$
$\displaystyle\
\displaystyle\frac{\langle(1+\delta_{1})(1+\delta_{2})(1+\delta_{3})(\bm{v}_{2}-\bm{v}_{1})\rangle}{\langle(1+\delta_{1})(1+\delta_{2})(1+\delta_{3})\rangle}$
$\displaystyle\simeq$ $\displaystyle\
\langle\delta_{1}\bm{v}_{2}\rangle-\langle\delta_{2}\bm{v}_{1}\rangle+\langle\delta_{3}\bm{v}_{2}\rangle-\langle\delta_{3}\bm{v}_{1}\rangle$
$\displaystyle=$ $\displaystyle\
\bar{w}(r_{12})\,\hat{\bm{r}}_{12}-\frac{1}{2}\left[\bar{w}(r_{23})\,\hat{\bm{r}}_{23}+\bar{w}(r_{31})\,\hat{\bm{r}}_{31}\right]\,.$
(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\bm{w}_{12}|\triangle_{123}\rangle_{\mathrm{t}}$, it is as follows

$\displaystyle\langle\bm{w}_{12}|\triangle_{123}\rangle_{\mathrm{t}}$
$\displaystyle=\langle\bm{w}_{12}\cdot\hat{\bm{r}}_{12}|\triangle_{123}\rangle_{\mathrm{t}}\,\hat{\bm{r}}_{12}+\langle\bm{w}_{12}\cdot\hat{\bm{t}}|\triangle_{123}\rangle_{\mathrm{t}}\,\hat{\bm{t}}$
$\displaystyle=R_{12}(\triangle_{123})\,\hat{\bm{r}}_{12}+T_{12}(\triangle_{123})\,\hat{\bm{t}}\;,$
(7)

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

$\displaystyle R_{12}(\triangle_{123})=\bar{w}(r_{12})$
$\displaystyle-\frac{1}{2}\Bigg{[}\bar{w}(r_{23})\,\cos\chi$
$\displaystyle-\bar{w}(r_{31})\,\frac{r_{12}+r_{23}\cos\chi}{\sqrt{r_{12}^{2}+r_{23}^{2}+2r_{12}r_{23}\cos\chi}}\Bigg{]}\;.$
(8)

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

$\displaystyle R_{23}(\triangle_{123})=\bar{w}(r_{23})$
$\displaystyle-\frac{1}{2}\Bigg{[}\bar{w}(r_{12})\,\cos\chi$
$\displaystyle-\bar{w}(r_{31})\,\frac{r_{23}+r_{12}\cos\chi}{\sqrt{r_{12}^{2}+r_{23}^{2}+2r_{12}r_{23}\cos\chi}}\Bigg{]}\;.$
(9)

Similar to Eq. (2), the three-point mean relative temperature difference from
kSZ can be written down as

$\frac{\Delta
T^{\mathrm{kSZ}}_{ij}(\triangle_{123})}{T_{\mathrm{cmb}}}\simeq-\tau\frac{R^{\mathrm{h}}_{ij}(\triangle_{123})}{c}\,,$
(10)

where $R^{\mathrm{h}}_{ij}(\triangle_{123})$ respresents 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}(\triangle_{123})$ (Kuruvilla & Aghanim, 2021). Based on the radial
mean relative velocities between pairs in a triplet, we can introduce a new
ratio statistic

$\displaystyle\mathcal{R}(\triangle_{123})$
$\displaystyle=\frac{\langle\bm{w}_{12}\cdot\hat{\bm{r}}_{12}|\triangle_{123}\rangle_{\mathrm{t}}}{\langle\bm{w}_{23}\cdot\hat{\bm{r}}_{23}|\triangle_{123}\rangle_{\mathrm{t}}}\,,$
(11)

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

$\displaystyle\mathcal{R}(\triangle_{123})$
$\displaystyle=\frac{R_{12}(\triangle_{123})}{R_{23}(\triangle_{123})}=\frac{R^{\mathrm{h}}_{12}(\triangle_{123})}{R^{\mathrm{h}}_{23}(\triangle_{123})}\equiv\frac{\Delta
T^{\mathrm{kSZ}}_{12}(\triangle_{123})}{\Delta
T^{\mathrm{kSZ}}_{23}(\triangle_{123})}\,.$ (12)

The above introduced 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}(\triangle_{123})$.

##  3 Data and analysis

###  3.1 Quijote simulation suite

In this work, we make use of the Quijote666https://quijote-
simulations.readthedocs.io/ (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}\,\mathrm{Gpc}$, and tracks the evolution of $512^{3}$ 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_{\mathrm{m}}=0.3175$, the baryonic matter density:
$\Omega_{\mathrm{b}}=0.049$, the primordial spectral index of the density
perturbations: $n_{\mathrm{s}}=0.9624$, the amplitude of the linear power
spectrum on the scale of $8\ h^{-1}\mathrm{Mpc}$: $\sigma_{8}=0.834$, and the
present-day value of the Hubble constant: $H_{0}\equiv
H(z=0)=100\,h\,\mathrm{km}\,\mathrm{s}^{-1}\mathrm{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 the fiducial cosmology. The variations are as follows:
{$\Omega^{+}_{\mathrm{m}},\Omega^{-}_{\mathrm{m}},\Omega^{+}_{\mathrm{b}},\Omega^{-}_{\mathrm{b}},n^{+}_{\mathrm{s}},n^{-}_{\mathrm{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), and has $512^{3}$ neutrino particles in
addition to the CDM particles. To compute the numerical derivatives with
respect to massive neutrinos, the Quijote suite provides an addition 500
random realisations for the fiducial cosmology, in which the ICs were also
generated using ZA.

In this work, we use 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 $M_{\mathrm{h}}>5\times
10^{13}\ h^{-1}\mathrm{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}\,\mathrm{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 use 30
realisations of the particle data (randomly down-sampled to $100^{3}$
particles) to compute $\mathcal{R}(\triangle_{123})$.

###  3.2 Fisher-matrix formalism

To quantify the error estimates on the cosmological parameters, we use 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\,,$
(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 write the Fisher information matrix as

$F_{\alpha\beta}=\frac{\partial\mkern
1.0mu\bm{\mathcal{R}}}{\partial\theta_{\alpha}}\cdot\hat{\mathbf{C}}^{-1}\cdot\frac{\partial\mkern
1.0mu\bm{\mathcal{R}}^{\mathsf{T}}}{\partial\theta_{\beta}}\,,$ (14)

where $\bm{\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 have neglected a term which 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

$\widetilde{\mathbf{C}}=\frac{1}{N_{\mathrm{sims}}-1}\sum_{i=1}^{N_{\mathrm{sims}}}\left(\bm{\mathcal{R}}_{i}-\overline{\bm{\mathcal{R}}}\right)\left(\bm{\mathcal{R}}_{i}-\overline{\bm{\mathcal{R}}}\right)^{\mathsf{T}}\,,$
(15)

where
$\overline{\bm{\mathcal{R}}}=N_{\mathrm{sims}}^{-1}\sum_{i=1}^{N_{\mathrm{sims}}}\bm{\mathcal{R}}_{i}$,
and $N_{\mathrm{sims}}$ denotes the total number of simulations used to
compute the covariance matrix (in this work $N_{\mathrm{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_{\mathrm{sims}}-N_{\mathrm{bins}}-2}{N_{\mathrm{sims}}-1}\
\widetilde{\mathbf{C}}^{-1}\,,$ (16)

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

We numerically compute 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 follows:
$\theta\equiv\\{\Omega_{\mathrm{m}},\Omega_{\mathrm{b}},h,n_{\mathrm{s}},\sigma_{8}\\}$,
we make use of the central difference approximation to compute the derivative
numerically

$\frac{\partial\mkern
1.0mu\bm{\mathcal{R}}}{\partial\theta}\simeq\frac{\bm{\mathcal{R}}(\theta+\mathrm{d}\theta)-\bm{\mathcal{R}}(\theta-\mathrm{d}\theta)}{2\
\mathrm{d}\theta}\,.$ (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\mkern 1.0mu\bm{\mathcal{R}}}{\partial
M_{\nu}}\simeq\frac{-\bm{\mathcal{R}}(M_{\nu}=0.4)+4\bm{\mathcal{R}}(M_{\nu}=0.2)-\bm{\mathcal{R}}(M_{\nu}=0)}{0.4}\,.$
(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 be consistent to compute the partial derivative,
we make use of another 500 realisations of fiducial cosmology (with
$M_{\nu}=0$ eV) simulation in which the initial conditions were also generated
using ZA.

##  4 Results

Figure 1: Top: comparing theoretical prediction (orange dashed line) for
$\mathcal{R}(\triangle_{123})$ using perturbation theory at LO against the
direct measurement (blue solid line) from the halo catalogs of the Quijote
suite of simulation. Bottom: residual showing the deviation of the theoretical
prediction from the direct measurement from the simulations. The blue shaded
region denotes the 5% region. The triangle configuration go from the smallest
being $\\{(40,45),(40,45),(40,45)\\}\,h^{-1}\mathrm{Mpc}$ to the largest which
corresponds to $\\{(115,120),(115,120),(115,120)\\}\,h^{-1}\mathrm{Mpc}$.

In Fig. 1 we show the direct measurements of $\mathcal{R}(\triangle_{123})$
from the 15,000 reference halo catalogs (solid blue line), and compare it
against the LO prediction (dashed orange line). We consider all triangular
configurations with $r_{\mathrm{min}}\in(40,45)$ and
$r_{\mathrm{max}}\in(115,120)$, and such that $r_{12}\geq r_{23}\geq r_{31}$.
All the separation scales have a bin width of 5 $h^{-1}\mathrm{Mpc}$. It thus
corresponds to a total of 766 triangular configurations, spanning from
configuration ‘0’ being the smallest (i.e.
$\triangle_{123}\in\\{(40,45),(40,45),(40,45)\\}\,h^{-1}\mathrm{Mpc}$) to
configuration ‘765’ being the largest
($\triangle_{123}\in\\{(115,120),(115,120),(115,120)\\}\,h^{-1}\mathrm{Mpc}$).
One can see from Eqs. (8) and (9) that 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}=1$ when this condition is met.
When comparing the theoretical predictions, we see that it is overall accurate
within 4–5% for configurations with all separation lengths greater than 55
$h^{-1}\mathrm{Mpc}$. As 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}\mathrm{Mpc}$. This thus motivates us to directly measure $\mathcal{R}$
from the simulations to compute the derivatives for the Fisher information
matrix.

Figure 2: Ratio of $\mathcal{R}$ at $\sigma_{8}^{+}=0.849$ to $\mathcal{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.

In addition $\mathcal{R}(\triangle_{123})$ is unaffected by variation in
$\sigma_{8}$ as mentioned earlier in Sect. 2.2. To demonstrate this, we
compute the ratio of $\mathcal{R}(\triangle_{123})$ for $\sigma^{+}_{8}=0.849$
and $\sigma^{-}_{8}=0.819$ using dark matter particles from 30 realisations,
and showcase it 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}(\triangle_{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^{\mathrm{h}}_{12}$,
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 halo we utilised the 15,000 catalogues. The shaded
regions represents the $1\sigma$ errors from the propagation of uncertainties
of the mean relative velocity statistics for matter, and halo. One 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 $\mathcal{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%, and hence for the purpose of Fisher
matrix formalism we consider $\mathcal{R}$ being independent of a (constant)
linear bias term.

Figure 3: The dashed (blue) line shows the bias for $R_{12}$, i.e. the mean
radial relative velocity between pairs 1 and 2 in a triplet. While the solid
(black) line shows the (weak) bias dependence of the ratio statistics
$\mathcal{R}$, which is equal to one within 1–2$\%$ for all triangular
configurations considered here.

Since $\mathcal{R}$ is found to be independent of $\sigma_{8}$ at the scales
we are probing (i.e. $r_{\mathrm{min}}\geq 40\ h^{-1}\mathrm{Mpc}$), this
renders the statistic at an unique position of being unaffected by the
degeneracy in the $M_{\nu}$-$\sigma_{8}$ parameter plane. We check the impact
of summed neutrino mass on $\mathcal{R}$ utilising three non-zero neutrino
mass in Fig. 4. The solid (blue) line shows the impact of $M_{\nu}=0.1$ eV on
$\mathcal{R}$ when compared to zero neutrino mass cosmology. Similarly the
dashed (orange) and dash-dotted (green) lines showcases the impact of
$M_{\nu}=0.2$ and $M_{\nu}=0.4$ eV, respectively. As can be seen, when the
neutrino mass increases there is a decrease in the infall velocity between a
pair in most of the triangular configurations. This is related to the free-
streaming of neutrinos as a result of them having large thermal velocities. As
a result, below the free-streaming scale neutrinos does not cluster, which
further slows down the collapse of the matter in general. 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 $\mathcal{R}$ for all the
configurations we measured, the maximal effect of suppression is seen in the
case when $M_{\nu}=0.4$ eV, and for the triangular configuration
$\\{(100,105),(50,55),(50,55)\\}\ h^{-1}\mathrm{Mpc}$ when compared to the
zero neutrino mass cosmology.

Figure 4: Effect of summed neutrino mass on $\mathcal{R}(\triangle_{123})$, as
measured directly from the simulations, when compared to zero neutrino mass
fiducial cosmology. The summed neutrino mass considered here are denoted in
the legend, and its units are given in eV.

###  4.1 Cosmological parameters

Figure 5: The correlation matrix (i.e. the covariance matrix of $\mathcal{R}$
normalised by its diagonal elements) computed using 15,000 realisations of the
Quijote simulations. The triangle configurations are same as in Fig. 1. Figure
6: Joint 68.3% (dark shaded contour) and 95.4% (light shaded contour) credible
region for all the pairs of cosmological model parameters at $z=0$.

We now turn our attention to the information content in the new ratio
statistic $\mathcal{R}$, and 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 for it are the partial derivatives of
$\mathcal{R}$ with respect to the cosmological model parameters and the
covariance matrix ($\mathbf{C}$). We show case the correlation matrix in Fig.
5, which is given as $\mathrm{C}_{ij}/\sqrt{\mathrm{C}_{ii}\mathrm{C}_{jj}}$,
wherein the covariance matrix is directly measured from the simulations using
15,000 realisations. We notice the presence of non-diagonal terms in it, being
positively correlated at similar triangular configurations while tending to be
negatively correlated as the configurations differ substantially.

We now take a look at the information content in $\mathcal{R}$ using the
Fisher information matrix formalism, as defined in Sect. 3.2. As mentioned, we
have computed both the elements of it directly from the simulations. Since the
bias dependence of $\mathcal{R}$ was shown to being very weak even at small
scales ($\sim$ 40–50 $h^{-1}\mathrm{Mpc}$) and being independent at large
scales ($\geq 80$ $h^{-1}\mathrm{Mpc}$), we do not consider the bias parameter
in the Fisher-matrix formalism. Thus the model parameters are
$\Omega_{\mathrm{m}}$, $\Omega_{\mathrm{b}}$, $h$, $n_{\mathrm{s}}$, and
$M_{\nu}$. We show the results of our Fisher forecast in Fig. 6, where the
dark and light shaded contours denote the 68.3% and 95.4% joint credible
region for all possible model parameters, respectively. The $1\sigma$
marginalised error for any model parameter $\theta_{\alpha}$ is given by
$\sqrt{F_{\alpha\alpha}^{-1}}$, and are as follows for the parameters we
considered:
$\\{\Omega_{\mathrm{m}},\Omega_{\mathrm{b}},h,n_{\mathrm{s}},M_{\nu}\\}\equiv\\{0.0158,0.0041,0.0391,0.0394,0.1175\\}$.

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). As a fair comparison, we use the
constraints obtained from them using $r_{\mathrm{min}}=40\
h^{-1}\mathrm{Mpc}$. It presents a factor of improvement of {6.2, 7.6, 9.8,
12.9, 8.87} for {$\Omega_{\mathrm{m}}$, $\Omega_{\mathrm{b}}$, $h$,
$n_{\mathrm{s}}$, $M_{\nu}$}, respectively over the mean pairwise velocity.
However when comparing against the constraints from
$R^{\mathrm{h}}_{12}(\triangle_{123})+R^{\mathrm{h}}_{23}(\triangle_{123})$
(i.e. the mean three-point relative velocities), $\mathcal{R}$ has its
constraining power shrunk by a factor of 1.34–1.44 for all the model
parameters. This could also due to the fact that the length of the data vector
is twice in $R^{\mathrm{h}}_{12}+R^{\mathrm{h}}_{23}$ when compared to
$\mathcal{R}$. This shrinkage in constraining power was also seen when
considering $R^{\mathrm{h}}_{12}$ and $R^{\mathrm{h}}_{23}$ separately as
compared to its combination in Kuruvilla & Aghanim (2021).

It is informative to ask 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). Comparing against the constraints from the redshift-space power
multipoles for $k_{\mathrm{max}}=0.2\,h\,\mathrm{Mpc}^{-1}$ (which is closest
to the $r_{\mathrm{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_{\mathrm{m}}$,
$\Omega_{\mathrm{b}}$, $h$, $n_{\mathrm{s}}$, $M_{\nu}$}, respectively.
However it slightly reduces to {1.4, 2.9, 3.1, 3.3, 2.5} when comparing
against the constraints from power spectrum when
$k_{\mathrm{max}}=0.5\,h\,\mathrm{Mpc}^{-1}$. This improvement over 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_{\mathrm{max}}=0.2\,h\,\mathrm{Mpc}^{-1}$, $\mathcal{R}$
still obtains a factor of improvement of {1.8, 2.9, 3.2, 3.1, 1.8} for
{$\Omega_{\mathrm{m}}$, $\Omega_{\mathrm{b}}$, $h$, $n_{\mathrm{s}}$,
$M_{\nu}$}, respectively. But when considering a larger set of triangular
configurations for the bispectrum with
$k_{\mathrm{max}}=0.5\,h\,\mathrm{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_{\mathrm{m}}$, $\Omega_{\mathrm{b}}$, $h$, $n_{\mathrm{s}}$,
$M_{\nu}$}, respectively. This is not surprising, as the bispectrum monopole
in this case, probes further into the nonlinear scales while $\mathcal{R}$ was
analysed for triangular configurations with separation scales of
$40\,h^{-1}\mathrm{Mpc}$ and above.

As it currently stands one of the limitations to applying the new statistic,
$\mathcal{R}(\triangle_{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 which can be measured
from either the peculiar velocity surveys or kSZ experiments). In the case of
the mean pairwise velocity, Ferreira et al. (1999) has shown how to construct
such an estimator. Similarly for the case of $R_{ij}(\triangle_{123})$, and
$\mathcal{R}(\triangle_{123})$, we will be constructing such an estimator in a
future work. In the case of the mean pairwise velocities, an alternative
estimator exists using 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. And
furthermore the analysis we presented here took 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 accurately.

Another caveat which we has not discussed 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}(\triangle_{123})$, as long as the ratio is
done using the same mass bin. On the other hand assuming a fixed cosmology, we
could consider a scenario where measuring $R^{\mathrm{h}}_{23}$ is fixed to a
high halo mass bin while measuring $R^{\mathrm{h}}_{12}$ for various mass bins
in Eq. (11). Thus, it could lead to potentially measuring the (scaled) mass
dependence of optical depth (and degenerate with the bias factor) from direct
kSZ experiment directly.

##  5 Conclusions

Determination of neutrino mass using cosmological observables have become one
of the main goals for the forthcoming cosmological surveys. However the two-
point statistics in general is affected by the $M_{\nu}$-$\sigma_{8}$
degeneracy, whether using clustering or relative velocity statistics which
limits the potential of constraining neutrino mass from cosmology. With
regards 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 in them. 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 extended the applications with the mean three-point relative
velocity statistics, and introduced a new ratio statistic $\mathcal{R}$ [Eq.
(11)] which is unaffected by $\sigma_{8}$. This enables to constrain neutrino
mass, in addition to other cosmological parameters, independent of
$\sigma_{8}$. Moreover, in the context of kSZ experiments this statistic is
independent of optical depth, hence circumventing 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 $\mathcal{R}$ will be bias
independent on linear scales. We verified it by measuring $\mathcal{R}$ for
both halos and matter, and found that the bias is consistent with one at 1–2%
for all triangular configurations we probed in this work
($r_{\mathrm{min}}=40\ h^{-1}\mathrm{Mpc}$ and $r_{\mathrm{max}}=120\
h^{-1}\mathrm{Mpc}$).

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

We used the Fisher-matrix formalism to quantify the information content in
$\mathcal{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 $\mathcal{R}$ has a factor
of 6.2–12.9 improvement when compared against the constraints obtained from
the mean pairwise velocity. When compared against the power spectrum and
bispectrum, it still achieves an improvement in the constraints with 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 shown that it can act as
robust cosmological observable which could lead to sizeable information gain
in comparison to the mean radial pairwise velocity. One of the limitation of
the kSZ experiments is the optical depth degeneracy, and breaking this
degeneracy requires some form of 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 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 et al. (2016) Abazajian, K. N., Adshead, P., Ahmed, Z., et al. 2016, arXiv e-prints, arXiv:1610.02743
  * Ade et al. (2019) Ade, P., Aguirre, J., Ahmed, Z., et al. 2019, J. Cosmology Astropart. Phys., 2019, 056
  * Agarwal et al. (2021) Agarwal, N., Desjacques, V., Jeong, D., & Schmidt, F. 2021, J. Cosmology Astropart. Phys., 2021, 021
  * Aker et al. (2021) Aker, M., Beglarian, A., Behrens, J., et al. 2021, arXiv e-prints, arXiv:2105.08533
  * Alam et al. (2021) Alam, S., Aubert, M., Avila, S., et al. 2021, Phys. Rev. D, 103, 083533
  * Anderson (2003) Anderson, T. 2003, An Introduction to Multivariate Statistical Analysis, Wiley Series in Probability and Statistics (Wiley)
  * Battaglia (2016) Battaglia, N. 2016, J. Cosmology Astropart. Phys., 2016, 058
  * Bayer et al. (2021) Bayer, A. E., Villaescusa-Navarro, F., Massara, E., et al. 2021, ApJ, 919, 24
  * Behnel et al. (2011) Behnel, S., Bradshaw, R., Citro, C., et al. 2011, Computing in Science & Engineering, 13, 13
  * Calafut et al. (2021) Calafut, V., Gallardo, P. A., Vavagiakis, E. M., et al. 2021, Phys. Rev. D, 104, 043502
  * Capozzi et al. (2017) Capozzi, F., Di Valentino, E., Lisi, E., et al. 2017, Phys. Rev. D, 95, 096014
  * Castorina et al. (2015) Castorina, E., Carbone, C., Bel, J., Sefusatti, E., & Dolag, K. 2015, J. Cosmology Astropart. Phys., 2015, 043
  * Chaves-Montero et al. (2021) Chaves-Montero, J., Hernández-Monteagudo, C., Angulo, R. E., & Emberson, J. D. 2021, MNRAS, 503, 1798
  * Chen et al. (2021) Chen, Z., Zhang, P., Yang, X., & Zheng, Y. 2021, arXiv e-prints, arXiv:2109.04092
  * Choudhury & Hannestad (2020) Choudhury, S. R. & Hannestad, S. 2020, J. Cosmology Astropart. Phys., 2020, 037
  * Chudaykin & Ivanov (2019) Chudaykin, A. & Ivanov, M. M. 2019, J. Cosmology Astropart. Phys., 2019, 034
  * De Bernardis et al. (2017) De Bernardis, F., Aiola, S., Vavagiakis, E. M., et al. 2017, J. Cosmology Astropart. Phys., 3, 008
  * de Salas et al. (2017) 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 et al. (2016) DESI Collaboration, Aghamousa, A., Aguilar, J., et al. 2016, arXiv e-prints, arXiv:1611.00036
  * Ferreira et al. (1999) Ferreira, P. G., Juszkiewicz, R., Feldman, H. A., Davis, M., & Jaffe, A. H. 1999, ApJ, 515, L1
  * Fisher (1995) Fisher, K. B. 1995, ApJ, 448, 494
  * Flender et al. (2017) Flender, S., Nagai, D., & McDonald, M. 2017, ApJ, 837, 124
  * Forero et al. (2014) Forero, D. V., Tórtola, M., & Valle, J. W. F. 2014, Phys. Rev. D, 90, 093006
  * García-Farieta et al. (2019) García-Farieta, J. E., Marulli, F., Veropalumbo, A., et al. 2019, MNRAS, 488, 1987
  * Gonzalez-Garcia et al. (2016) Gonzalez-Garcia, M. C., Maltoni, M., & Schwetz, T. 2016, Nuclear Physics B, 908, 199
  * Gualdi & Verde (2020) Gualdi, D. & Verde, L. 2020, J. Cosmology Astropart. Phys., 2020, 041
  * Hahn & Villaescusa-Navarro (2021) Hahn, C. & Villaescusa-Navarro, F. 2021, J. Cosmology Astropart. Phys., 2021, 029
  * Hahn et al. (2020) Hahn, C., Villaescusa-Navarro, F., Castorina, E., & Scoccimarro, R. 2020, J. Cosmology Astropart. Phys., 2020, 040
  * Hand et al. (2012) Hand, N., Addison, G. E., Aubourg, E., et al. 2012, Physical Review Letters, 109, 041101
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 585
  * Hartlap et al. (2007) Hartlap, J., Simon, P., & Schneider, P. 2007, A&A, 464, 399
  * Heavens (2009) Heavens, A. 2009, arXiv e-prints, arXiv:0906.0664
  * Hernández-Monteagudo et al. (2015) Hernández-Monteagudo, C., Ma, Y.-Z., Kitaura, F. S., et al. 2015, Phys. Rev. Lett., 115, 191301
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 9
  * Juszkiewicz et al. (1998) Juszkiewicz, R., Fisher, K. B., & Szapudi, I. 1998, ApJ, 504, L1
  * Kaufmann (1967) Kaufmann, G. M. 1967, Some Bayesian Moment Formulae, Report No. 6710. Centre for Operations Research and Econometrics (Catholic University of Louvain, Heverlee)
  * Keisler & Schmidt (2013) Keisler, R. & Schmidt, F. 2013, ApJ, 765, L32
  * Kodwani et al. (2019) Kodwani, D., Alonso, D., & Ferreira, P. 2019, The Open Journal of Astrophysics, 2, 3
  * Kuruvilla & Aghanim (2021) Kuruvilla, J. & Aghanim, N. 2021, A&A, 653, A130
  * Kuruvilla et al. (2020) Kuruvilla, J., Aghanim, N., & McCarthy, I. G. 2020, A&A, 644, A170
  * Kuruvilla & Porciani (2020) Kuruvilla, J. & Porciani, C. 2020, J. Cosmology Astropart. Phys., 2020, 043
  * Laureijs et al. (2011) Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, arXiv e-prints, arXiv:1110.3193
  * Lesgourgues & Pastor (2012) Lesgourgues, J. & Pastor, S. 2012, Adv. High Energy Phys., 2012, 608515
  * Li et al. (2018) Li, Y.-C., Ma, Y.-Z., Remazeilles, M., & Moodley, K. 2018, Phys. Rev. D, 97, 023514
  * Madhavacheril et al. (2019) Madhavacheril, M. S., Battaglia, N., Smith, K. M., & Sievers, J. L. 2019, Phys. Rev. D, 100, 103532
  * Mueller et al. (2015) Mueller, E.-M., de Bernardis, F., Bean, R., & Niemack, M. D. 2015, Phys. Rev. D, 92, 063501
  * Nguyen et al. (2020) Nguyen, N.-M., Jasche, J., Lavaux, G., & Schmidt, F. 2020, J. Cosmology Astropart. Phys., 2020, 011
  * Planck Collaboration et al. (2016) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 586, A140
  * Planck Collaboration et al. (2020) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, A&A, 641, A6
  * Reid & White (2011) Reid, B. A. & White, M. 2011, MNRAS, 417, 1913
  * Saito et al. (2008) Saito, S., Takada, M., & Taruya, A. 2008, Phys. Rev. Lett., 100, 191301
  * Samushia et al. (2021) Samushia, L., Slepian, Z., & Villaescusa-Navarro, F. 2021, MNRAS, 505, 628
  * Schaan et al. (2021) Schaan, E., Ferraro, S., Amodeo, S., et al. 2021, Phys. Rev. D, 103, 063513
  * Schaan et al. (2016) Schaan, E., Ferraro, S., Vargas-Magaña, M., et al. 2016, Phys. Rev. D, 93, 082002
  * Smith et al. (2018) Smith, K. M., Madhavacheril, M. S., Münchmeyer, M., et al. 2018, arXiv e-prints, arXiv:1810.13423
  * Soergel et al. (2016) Soergel, B., Flender, S., Story, K. T., et al. 2016, MNRAS, 461, 3172
  * Soergel et al. (2018) Soergel, B., Saro, A., Giannantonio, T., Efstathiou, G., & Dolag, K. 2018, MNRAS, 478, 5320
  * Spergel et al. (2015) Spergel, D., Gehrels, N., Baltay, C., et al. 2015, arXiv e-prints, arXiv:1503.03757
  * Springel (2005) Springel, V. 2005, MNRAS, 364, 1105
  * Sunyaev & Zeldovich (1980) Sunyaev, R. A. & Zeldovich, I. B. 1980, MNRAS, 190, 413
  * Sunyaev & Zeldovich (1972) Sunyaev, R. A. & Zeldovich, Y. B. 1972, Comments on Astrophysics and Space Physics, 4, 173
  * Tanimura et al. (2021) Tanimura, H., Zaroubi, S., & Aghanim, N. 2021, A&A, 645, A112
  * Tegmark et al. (1997) Tegmark, M., Taylor, A. N., & Heavens, A. F. 1997, ApJ, 480, 22
  * Uhlemann et al. (2020) Uhlemann, C., Friedrich, O., Villaescusa-Navarro, F., Banerjee, A., & Codis, S. 2020, MNRAS, 495, 4006
  * Vagnozzi et al. (2018) Vagnozzi, S., Dhawan, S., Gerbino, M., et al. 2018, Phys. Rev. D, 98, 083501
  * Verde (2010) Verde, L. 2010, Statistical Methods in Cosmology, Vol. 800 (Springer Verlag), 147–177
  * Villaescusa-Navarro et al. (2018) Villaescusa-Navarro, F., Banerjee, A., Dalal, N., et al. 2018, ApJ, 861, 53
  * Villaescusa-Navarro et al. (2020) Villaescusa-Navarro, F., Hahn, C., Massara, E., et al. 2020, ApJS, 250, 2
  * Wong (2008) Wong, Y. Y. Y. 2008, J. Cosmology Astropart. Phys., 2008, 035
  * Wong (2011) Wong, Y. Y. Y. 2011, Annual Review of Nuclear and Particle Science, 61, 69
  * Yankelevich & Porciani (2019) Yankelevich, V. & Porciani, C. 2019, MNRAS, 483, 2078
  * Yasini et al. (2019) Yasini, S., Mirzatuny, N., & Pierpaoli, E. 2019, ApJ, 873, L23