# Eigen-decomposition of Covariance matrices : An application to the BAO
Linear Point

Jaemyoung (Jason) Lee  0000-0001-6633-9793<EMAIL_ADDRESS>Department
of Physics and Astronomy, University of Pennsylvania, Philadelphia, PA 19104,
U.S.A.  Farnik Nikakhtar  0000-0002-3641-4366 Department of Physics, Yale
University, New Haven, CT 06511, U.S.A  Aseem Paranjape  0000-0001-6832-9273
Inter-University Centre for Astronomy and Astrophysics, Ganeshkhind, Pune
411007, India  Ravi K. Sheth  0000-0002-2330-0917 Center for Particle
Cosmology, University of Pennsylvania, Philadelphia, PA 19104 – USA

###### Abstract

The Baryon Acoustic Oscillation (BAO) feature in the two-point correlation
function (TPCF) of discrete tracers such as galaxies is an accurate standard
ruler. The covariance matrix of the TPCF plays an important role in
determining how the precision of this ruler depends on the number density and
clustering strength of the tracers, as well as the survey volume. An eigen-
decomposition of this matrix provides an objective way to separate the
contributions of cosmic variance from those of shot-noise to the statistical
uncertainties. For the signal-to-noise levels that are expected in ongoing and
next-generation surveys, the cosmic variance eigen-modes dominate. These modes
are smooth functions of scale, meaning that: they are insensitive to the
modest changes in binning that are allowed if one wishes to resolve the BAO
feature in the TPCF; they provide a good description of the correlated
residuals which result from fitting smooth functional forms to the measured
TPCF; they motivate a simple but accurate approximation for the uncertainty on
the Linear Point (LP) estimate of the BAO distance scale. This approximation
allows one to quantify the precision of the BAO distance scale estimate
without having to generate a large ensemble of mock catalogs and explains why:
the uncertainty on the LP does not depend on the functional form fitted to the
TPCF or the binning used; the LP is more constraining than the peak or dip
scales in the TPCF; the evolved TPCF is less constraining than the initial
one, so that reconstruction schemes can yield significant gains in precision.

large-scale structure of Universe

###### pacs:

##  I Introduction

Some of the tightest constraints on the cosmological distance - redshift
relation come from the BAO feature in the pair correlation function Eisenstein
et al. (2005); Anselmi et al. (2018a); Alam et al. (2021); Abbott et al.
(2024); Adame et al. (2024a). This has led to interest in the precision with
which the pair correlation function (TPCF) can be measured, and how this
precision translates into uncertainties on the distance scale estimate. As a
result, there is significant interest in understanding the covariance between
pair counts on different scales.

On the $\sim 140$ Mpc scales of most interest to BAO cosmology, the Gauss-
Poisson approximation to the covariance is rather accurate Smith et al.
(2008a); Grieb et al. (2016). In this approximation, three different terms
contribute: one is a purely cosmological term, the other is a pure shot-noise
term, and the third is a combination of the two. The first part of our paper
is devoted to a study of the relative importance of these terms. We address
this by rotating the covariance matrix into diagonal form, checking how many
eigenvectors contribute significantly to the total covariance, and then
looking at those eigenvectors. It turns out that this provides a simple way of
seeing which term dominates and when, as well as for understanding the shapes
of the eigenvectors. We use this insight to explore how binning of the pair
counts (width of the rectangular bins, or different bin shapes) affects the
structure of the covariance matrix.

The second half of this paper applies these insights to a particular estimator
of the BAO scale: the Linear Point (LP). The LP feature in the correlation
function of dark matter or galaxies can be used as a standard cosmological
ruler Anselmi et al. (2016, 2019, 2018a); Parimbelli et al. (2021); O’Dwyer et
al. (2020); Nikakhtar et al. (2021a); Lee et al. (2024); Paranjape and Sheth
(2022); He et al. (2023). The LP lies midway between the peak and dip values
in $\xi(r)$, the two-point correlation function:

$r_{\rm LP}\equiv\frac{r_{\rm peak}+r_{\rm dip}}{2}.$ (1)

Evidently, the precision with which $r_{\rm LP}$ can be estimated from data
depends on the covariance between the $r_{\rm peak}$ and $r_{\rm dip}$
estimates. In turn, this depends on the covariance matrix of the measurements,
which depends on the widths of bins in which pairs were counted (or, more
generally, the bin shapes themselves). This has led to significant
computational efforts simply to determine the optimal bin width Vargas-Magaña
et al. (2014); Anselmi et al. (2018b).

Moreover, $r_{\rm LP}$ is typically estimated by fitting a pre-determined
functional form to the measured $\xi$ (e.g. polynomials, Chebyshev
polynomials, generalized half-integer Laguerre functions). The associated
error bars would then appear to be closely tied to this functional form (e.g.
Equations 2.6 and 9 in Parimbelli et al. (2021) and Nikakhtar et al. (2021a)).
However, in practice, provided that the fits are good, neither the $r_{\rm
LP}$ estimates nor their error bars depend strongly on which functional form
is fit. Our analysis of the covariance matrix allows us to provide a rather
general estimate of the expected precision that is not tied to a particular
functional form. It also allows us to address a closely related question. In
principle, the inflection point $r_{\rm infl}$, the scale on which
$d^{2}\xi/dr^{2}=0$, could also be used as a standard rod (Anselmi et al.,
2016). Previous work has suggested that it is less robust than the LP (Anselmi
et al., 2018b; Nikakhtar et al., 2021b); our analysis provides some insight
into why this is so.

This paper is organized as follows. Section II describes how the eigenvalues
and eigenvectors of the TPCF covariance matrix change as the shot-noise
increases, and then use this to provide a simple estimate of the error on the
LP. Section III shows how our results depend on the binning. Section IV
summarizes our conclusions.

##  II Methods and Results

Because the neighboring bins of the TPCF amplitudes are correlated, the
covariance matrix of the bin counts is not diagonal. Here below, we describe
in detail how we use the structure of the covariance matrix to estimate
realistic error bars on the BAO distance scale.

Where necessary to illustrate our results, we will use a comoving volume of
$5\times(1.024)^{3}$ $h^{-3}$Gpc3 in a flat $\Lambda$CDM model with
$(\Omega_{m,0},\Omega_{b,0})=(0.281,0.046)$, and
$(h,n_{s},\sigma_{8})=(0.697,0.971,0.842)$ as in Lee et al. (2024). The
associated linear theory values of $r_{\rm LP}$ and $r_{\rm infl}$ for the
dark matter are $97.154~{}h^{-1}$Mpc and $97.635~{}h^{-1}$Mpc, respectively.
For easy comparison with Refs.Nikakhtar et al. (2021a); Lee et al. (2024), we
focus on biased tracers at $z=0.5$ (also denoted as $z=0.5057$ following Lee
et al. (2024)). Although we explore other combinations of number density and
clustering strength, our fiducial choice has $\bar{n}=3.2\times
10^{-4}/(h^{-1}{\rm Mpc})^{3}$ and linear bias factor $b\approx 1.97$, which
is similar to the Baryon Oscillation Spectroscopic Survey (BOSS) and the Dark
Energy Spectroscopic Instrument (DESI) survey Alam et al. (2015); Adame et al.
(2024b). The combination $\bar{n}P(k_{\rm max})$, where $k_{\rm max}$ is the
scale on which $P(k)=b^{2}\,P_{\rm Lin}(k)$ is maximum, is sometimes used as a
crude measure of whether the BAO clustering signal is dominated by shot-noise.
Our fiducial choice has $\bar{n}P(k_{\rm max})\sim 5$; shot-noise dominates
for values smaller than unity. While we provide our formalism in terms of the
real-space correlation function, we use the redshift-space monopole in our
figures to show that our methodology is valid even under redshift-space
distortions. This makes $b_{\rm eff}=\sqrt{(b^{2}+2bf/3+f^{2}/5)}$ with
$f=d\ln D/d\ln a$ being the linear theory growth rate Kaiser (1987). With
$b\sim 1.97$ at $z\sim 0.5$, $b_{\rm eff}\sim 2.23$.

###  II.1 Gauss-Poisson approximation to covariance matrix

We begin with the two-point correlation function, which is related to the
power spectrum $P(k)$ by

$\xi(s)\equiv\int\frac{dk}{k}\frac{k^{3}P(k)}{2\pi^{2}}\,j_{0}(ks).$ (2)

A crude model for nonlinear evolution sets $P(k)=P_{\rm Lin}(k)\,{\rm
e}^{-k^{2}\Sigma^{2}}$ (e.g. Crocce and Scoccimarro, 2008), so the nonlinear
correlation function is a ‘smeared’ version of the linear one. At $z=0.5$ in
our fiducial cosmology, $\Sigma\approx 4.7~{}h^{-1}$Mpc for the real-space
TPCF; it is slightly larger, $\Sigma_{\rm eff}\approx 6.3~{}h^{-1}$Mpc, for
the biased redshift-space monopole Nikakhtar et al. (2021b); Paranjape and
Sheth (2023).

If the correlation function is estimated by counting pairs separated by
$s\pm\Delta s/2$ in a discrete set of $N_{\rm tot}$ particles distributed in a
volume $V_{s}$, then the TPCF covariance matrix described by the ‘Gauss-
Poisson’ approximation is given by Smith et al. (2008b); Grieb et al. (2016);
Parimbelli et al. (2021):

$C_{ij}=\frac{2}{V_{s}}\int_{0}^{\infty}\frac{dk~{}k^{2}}{2\pi^{2}}\bar{j_{0}}(ks_{i})\bar{j_{0}}(ks_{j})\big{[}P(k)+\frac{1}{\bar{n}}\big{]}^{2}$
(3)

where $\bar{n}\equiv N_{\rm tot}/V_{s}$ is the survey number density and
$\bar{j_{0}}$ is the bin-averaged spherical Bessel function:

$\bar{j_{0}}(ks_{i})\equiv\frac{4\pi}{V_{s_{i}}}\int_{s_{i}-\Delta
s/2}^{s_{i}+\Delta s/2}s^{2}j_{0}(ks)ds,$ (4)

with the volume $V_{s_{i}}=4\pi/3(s_{i,\rm{max}}^{3}-s_{i,\rm{min}}^{3})$,
$s_{i}$ being the midpoint of bin $i$, and $\Delta s$ being the bin-size. The
shot-noise only term proportional to $1/\bar{n}^{2}$ only contributes when
$i=j$, i.e., to the error bar in a single bin. The other two terms describe
the covariance between bins, and come from ‘cosmic variance.’ This covariance
must be accounted for when estimating the uncertainty on the BAO scale.

Figure 1: Fractional contribution of the sorted eigenvalues $\lambda_{j}$ to
the total variance, $\sum_{1}^{i}{\lambda_{j}}/\sum_{1}^{30}{\lambda_{j}}$,
for a range of choices of the ratio of ‘signal’ ($P(k)$ amplitude) to ‘noise’
($1/\bar{n}$). Fiducial and idealized no shot-noise cases are very similar;
more than 95% of the total variance comes from the first 3 modes. More modes
matter when the noise dominates. Smooth dashed curve shows our analytic
expression for the pure noise case: $(2/\bar{n}V)/(\bar{n}V_{i})$. Dotted
curve is for $b\to b-1$ but fiducial shot-noise. Figure 2: The first 9
eigenvectors $\sqrt{\lambda_{i}}\Lambda_{i}(s)$ at $z=0.5057$, for default
shot-noise $\bar{n}=3.2\times 10^{-4}/(h^{-1}{\rm Mpc})^{3}$ (solid) and no
shot-noise (dashed). The first 3 modes (left hand panel), which contribute
most to the total variance, are much less sensitive to shot-noise. Higher
order modes (right hand panel, note the factor of 4 difference in the $y$-axis
scale) are essentially zero when there is no shot-noise, and resemble delta-
functions which peak at different scales, followed by oscillations, when shot-
noise is present.

###  II.2 Eigen-decomposition of the covariance matrix

To illustrate our results, we now use Eq. 3 to generate $C_{ij}$, for 30 non-
overlapping bins of $\Delta r=2~{}h^{-1}$ Mpc, running from
$[60-120]~{}h^{-1}$Mpc, with the fiducial values of background cosmology,
redshift, survey volume, biased tracer number density, and clustering strength
mentioned at the start of this section.

Next, we diagonalize $C_{ij}$. The eigenvectors, which we denote
$\Lambda_{i}(s)$, provide a set of orthogonal shape functions, whose relative
importance is set by the eigenvalues $\lambda_{i}$. Before we consider the
interplay of $P$ and $\bar{n}$ in determining these eigen-modes, note that the
survey volume $V_{s}$ only appears as an overall scaling. It scales the
eigenvalues up and down but keeps their ratios fixed, and does not affect the
eigen-shapes. That said, $V_{s}$ is important because it does not enter in the
definition of $\xi$ itself. Therefore, larger $V_{s}$ means that the eigen-
modes will have smaller amplitudes compared to $\xi$. This will be important
below.

Before we look at the shapes, Fig. 1 shows the fractional contribution of the
eigenvalues $\lambda_{i}$ (ordered from largest to smallest) to the total
variance. The various curves show different choices for the relative
contributions of ‘signal’ $P(k)$ and ‘shot-noise’ $1/\bar{n}$.

The lowermost curve is for the pure shot-noise limit (we have set $P(k)=0$).
This case is analytic: $C_{ii}$ is diagonal, with entries
$(2/\bar{n}V)/(\bar{n}V_{i})$, where $V_{i}$ is the volume of the $i$th bin.
For bins of width $\Delta r=r_{i+1}-r_{i}$ that are equally spaced with
spacing $\epsilon\Delta r$ (typically $\epsilon=1$ but we will see later why
the more general case is interesting) where $r_{i}$ is the lower bound of the
$i$th bin, $V_{i}=(4\pi/3)(\Delta
r)^{3}[(j\epsilon-(\epsilon-1)/2)^{3}-(j\epsilon-1-(\epsilon-1)/2)^{3}]$ where
$j=i+60/(\epsilon\Delta r)$ and the bin centers are given by
$\frac{(r_{i}+r_{i+1})}{2}=60+(i-1)\epsilon\Delta r+\Delta r/2$ since our bins
start at $60~{}h^{-1}$Mpc. Note that all eigenvectors matter. These
eigenvectors are delta functions, one for each bin, centered on the middle of
the bin.

In contrast, the uppermost curve shows the case in which $1/\bar{n}=0$: here
$C_{ij}$ is completely determined by the $P^{2}$ term in Eq. 3. Notice that
now the variance is dominated by just a few eigenvalues/eigenvectors. The
intermediate curves show different choices for the $1/\bar{n}$. The ‘fiducial’
choice ($\bar{n}P(k_{\rm max})\approx 5$) is very similar to the no-noise
case. However, as the noise increases, more modes begin to matter.111The cross
term in $C_{ij}$ is also analytic: it is a smoothed version of $\xi$, with
smoothing depending on scales $i$ and $j$, but the expression is lengthy so we
have not written it explicitly here.

Fig. 2 shows the corresponding eigenvectors. Notice that the $n$th eigenvector
has $n-1$ zero-crossings, at least for the first few $n$. It is striking that
the mode with 4 zero-crossings divides the 60-120$~{}h^{-1}$Mpc range up into
patches that are approximately the size of the BAO feature itself. Presumably
this is because the same $P(k)$ appears in both $\xi$ and $C_{ij}$. Dashed
curves are for the case of no shot-noise, and solid curves have the fiducial
shot-noise. Recall that these are the cases that are dominated by the first
few eigenvalues, and the associated eigenvectors are very similar and very
smooth. The smoothness is consistent with the expectation that terms
contributing to cosmic variance should be smooth functions of scale. However,
the similarity is particularly interesting here: it suggests that the
eigenvectors for the no-noise case remain interesting even in the presence of
fiducial noise.

The other higher-order eigenvectors, which contribute little to the total
variance, are more strongly modified by the presence of shot-noise. To gain
some insight, recall that the pure shot-noise eigenvectors would be a set of
delta functions, each centered on a bin. But, when $P$ is significant, these
delta functions are now approximately rotated into the basis in which the
$P^{2}$ term is diagonal. This mixes the delta functions, and is why these
higher order modes display oscillations. We will exploit this relatively clean
separation into cosmic variance vs. shot-noise dominated eigen-modes in the
next section.

###  II.3 Eigen-decomposition of correlation function realizations

We can write one realization $\Xi$ of the real-space TPCF as:

$\Xi(s)=\xi(s)+\sum_{i}g_{i}\Lambda_{i}(s),$ (5)

where $\xi(s)$ is given by Eq. 2, and the $g_{i}$’s are independent Gaussian
random variates with variance $\lambda_{i}$ and mean zero, so the other terms
represent the (correlated) scatter around the mean.

The symbols in Fig. 3 show how the shape of $\Xi$ changes as more modes are
added to $\xi$, for one realization where we have assumed the fiducial noise
and bias. The changes are relatively mild because $V_{s}$ is sufficiently
large that the $\lambda_{i}$ are small. To highlight the differences as more
modes are added, the symbols in Fig. 4 show the total residual from the mean,
$\Xi(s)-\xi(s)=\sum_{i}g_{i}\Lambda_{i}(s)$, for this same realization. The
other curves show the contribution from modes 1 to 4
($\sum_{i=1}^{i=4}g_{i}\Lambda_{i}(s)$), and from 5 onwards
($\sum_{i=5}^{i=30}g_{i}\Lambda_{i}(s)$). Clearly, the first 4 modes capture
most of the residual, including the small change in shape, while the sum of
modes 5 and onwards is mostly uncorrelated noise (small amplitude oscillations
around zero). This just illustrates what Fig. 1 showed: the higher order modes
are not particularly important.

Figure 3: $\xi(s)$, $\Xi(s)$, and Eq. 5 with the first 2 and 4 modes. For a
survey volume smaller than our fiducial $V_{s}$ of $\sim
5~{}(h^{-1}\rm{Gpc})^{3}$, the deviations of the other 3 curves from the blue
$\xi(s)$ curve would be larger.  Figure 4: Total residual $\Xi(s)-\xi(s)$ for
one TPCF realization (symbols), and the contributions to it from modes 1 to 4
(dot-dashed), and from 5 onwards (dashed).

###  II.4 Truncation of modes and basis functions for fitting the TPCF

Because we see that the sum of the first 4 modes captures most of the
residual, while the remaining modes are mostly ‘noise,’ it is interesting to
consider restricting the sum in Eq. 5 to include only the first $\sim 4$
modes. Evidently, this removes the shot-noise dominated fluctuations from the
realization of $\Xi$, leaving a smoother curve. In essence, this is the smooth
curve one is after when ‘fitting’ the correlation. This is a non-trivial
statement, since $\xi$ $+$ the first few eigenvectors, while smoother than the
full $\Xi$, will not generically have the same shape as $\xi$. Nevertheless,
since typical datasets were designed (i.e. $V_{s}$ is large enough) so that
cosmic variance does not dominate, the amplitude of these ‘cosmic variance’
modes is small compared to the amplitude of $\xi$ itself, so the correction to
the shape is small (c.f. Fig. 3).

In the same vein, suppose one is interested in derivatives of the correlation
function. Although

$\Xi^{\prime}=\xi^{\prime}+\sum_{i}g_{i}\Lambda_{i}^{\prime},$ (6)

if we include all 30 terms in the sum, then this will be like differentiating
a single measurement of the correlation function. However, it is well known
that one should not differentiate a noisy measurement; rather, one should
first fit a smooth functional form to the measurement and then differentiate
the fit. In the present context, our model for this procedure is to assert
that one is not interested in $d\Xi/ds$ when the sum includes all 30 terms;
rather, one should only include the first 4 (really, the ones which account
for, say, 90% of the variance).

The virtue of this point of view is that this shape is clearly determined
solely by the shape of $P(k)$; in particular, it makes no reference to the set
of basis functions which one wishes to fit to $\Xi$ (simple polynomials?
Laguerre functions? etc.). This is attractive, since a reasonable concern is
whether the set of basis functions which worked for one underlying $P(k)$ will
also work for another (polynomials for one, Laguerres for another?). Here, the
point is that one should think of the eigenvectors $\Lambda_{i}$ as being the
most appropriate set of basis functions, since these are clearly determined by
the shape of $P(k)$.

###  II.5 Model-independent error estimates on BAO distance scale

We will now use this insight to discuss how one might quantify uncertainties
on estimates of the BAO scales $r_{\rm{LP}}$, $r_{\rm{Dip}}$ and
$r_{\rm{Peak}}$. We also consider, $r_{\rm{infl}}$, the scale on which the
second derivative vanishes, as an alternative to $r_{\rm LP}$.

We start with Eq. 6, but restrict the sum to the first few terms (the ones
which account for, say, 90% of the variance). Next, note that the scale
$s_{\rm{max}}$ where $\Xi^{\prime}=0$ is not necessarily the same as
$r_{\rm{max}}$, where $\xi^{\prime}=0$. Assuming $\Delta_{\rm{max}}\equiv
s_{\rm{max}}-r_{\rm{max}}\ll 1$, we have when $\Xi^{\prime}=0$,

$\begin{split}0&=\xi^{\prime}(s_{\rm{max}})+\sum_{i}g_{i}\Lambda_{i}^{\prime}(s_{\rm{max}})\\\
&\approx\xi^{\prime}(r_{\rm{max}})+\Delta_{\rm{max}}\xi^{\prime\prime}(r_{\rm{max}})+\sum_{i}g_{i}\Lambda_{i}^{\prime}(s_{\rm{max}})\\\
&\approx\Delta_{\rm{max}}\xi^{\prime\prime}(r_{\rm{max}})+\sum_{i}g_{i}\big{[}\Lambda_{i}^{\prime}(r_{\rm{max}})+\Delta_{\rm{max}}\Lambda_{i}^{\prime\prime}(r_{\rm{max}})\big{]},\end{split}$
(7)

where we have used that $\xi^{\prime}=0$ at $r_{\rm{max}}$. Isolating
$\Delta_{\rm{max}}$ yields

$\frac{\Delta_{\rm{max}}}{r_{\rm
max}}=-\frac{\sum_{i}g_{i}\,(\Lambda_{i}^{\prime}/r_{\rm{max}})/\xi^{\prime\prime}}{1+\sum_{i}g_{i}(\Lambda_{i}^{\prime\prime}/\xi^{\prime\prime})}.$
(8)

In practice, on the peak and dip scales,
$(\Lambda_{i}^{\prime\prime}/\xi^{\prime\prime})\leq 0.1$ or so. Since
$\langle g_{i}^{2}\rangle=\lambda_{i}$ is also small, we can neglect the term
in the denominator and approximate

$\bigg{\langle}\frac{\Delta_{\rm{max}}^{2}}{r_{\rm{max}}^{2}}\bigg{\rangle}\approx\sum_{i}\lambda_{i}\,\bigg{(}\frac{r_{\rm
max}\Lambda_{i}^{\prime}(r_{\rm{max}})}{r_{\rm{max}}^{2}\xi^{\prime\prime}(r_{\rm{max}})}\bigg{)}^{2}.$
(9)

Eq. 9 shows that the root-mean-square (RMS) of $\Delta_{\rm{max}}$ increases
as $|\xi^{\prime\prime}|$ decreases. Since the nonlinear TPCF is more smeared
(i.e. less curved at the peak and dip scales) than the linear theory TPCF
(recall discussion of Eq. 2), we expect the uncertainty in the peak scale to
be larger in the evolved field (at lower redshifts). We return to this point
in the next subsection.

The same logic can be applied to the dip scale, so

$\bigg{\langle}\frac{\Delta_{\rm{min}}^{2}}{r_{\rm{min}}^{2}}\bigg{\rangle}\approx\sum_{i}\lambda_{i}\,\bigg{(}\frac{\Lambda_{i}^{\prime}(r_{\rm{min}})}{r_{\rm{min}}\xi^{\prime\prime}(r_{\rm{min}})}\bigg{)}^{2}.$
(10)

Hence, just as for $\Delta_{\rm{min}}$, the RMS of $\Delta_{\rm{min}}$ in the
evolved TPCF is larger than the linear theory value.

Finally, for the error on the LP scale,

$\Delta_{\rm{LP}}\equiv\frac{\Delta_{\rm{max}}+\Delta_{\rm{min}}}{2},$ (11)

we have:

$\langle\Delta_{\rm{LP}}^{2}\rangle\approx\frac{\langle\Delta_{\rm{max}}^{2}\rangle}{4}+\frac{\langle\Delta_{\rm{min}}^{2}\rangle}{4}+\sum_{i}\frac{\lambda_{i}}{2}\frac{\Lambda_{i}^{\prime}(r_{\rm{max}})}{\xi^{\prime\prime}(r_{\rm{max}})}\frac{\Lambda_{i}^{\prime}(r_{\rm{min}})}{\xi^{\prime\prime}(r_{\rm{min}})}.$
(12)

Because $\xi^{\prime\prime}(r_{\rm{max}})$ and
$\xi^{\prime\prime}(r_{\rm{min}})$ have opposite signs (by definition), the
variance of $\Delta_{\rm{LP}}$ is smaller than either $\Delta_{\rm{Dip}}$ or
$\Delta_{\rm{Peak}}$. This demonstrates why the LP is a more precise probe
than either $r_{\rm{Dip}}$ and $r_{\rm{Peak}}$.

Similarly for the inflection point: if $r_{\rm infl}$ is the scale where
$\xi^{\prime\prime}=0$, and $s_{\rm infl}$ is the scale where
$\Xi^{\prime\prime}=0$ (when the sum which defines $\Xi$ is truncated to only
include the eigenvalues whose eigenvectors contribute $\sim 90\%$ of the
variance), we would set

$\begin{split}0&=\xi^{\prime\prime}(s_{\rm{infl}})+\sum_{i}g_{i}\Lambda_{i}^{\prime\prime}(s_{\rm{infl}})\\\
&\approx\Delta_{\rm{infl}}\xi^{\prime\prime\prime}(r_{\rm{infl}})+\sum_{i}g_{i}\big{[}\Lambda_{i}^{\prime\prime}(r_{\rm{infl}})+\Delta_{\rm{infl}}\Lambda_{i}^{\prime\prime\prime}(r_{\rm{infl}})\big{]},\end{split}$
(13)

where $\Delta_{\rm infl}=s_{\rm infl}-r_{\rm infl}$. Thus,

$\frac{\Delta_{\rm infl}}{r_{\rm
infl}}=-\frac{\sum_{i}g_{i}(\Lambda_{i}^{\prime\prime}/r_{\rm
infl})/\xi^{\prime\prime\prime}}{1+\sum_{i}g_{i}(\Lambda_{i}^{\prime\prime\prime}/\xi^{\prime\prime\prime})}\approx-\sum_{i}g_{i}\,\frac{\Lambda_{i}^{\prime\prime}}{r_{\rm
infl}\,\xi^{\prime\prime\prime}}.$ (14)

so

$\bigg{\langle}\frac{\Delta_{\rm{infl}}^{2}}{r_{\rm{infl}}^{2}}\bigg{\rangle}\approx\sum_{i}\lambda_{i}\,\bigg{(}\frac{\Lambda_{i}^{\prime\prime}(r_{\rm{infl}})}{r_{\rm{infl}}\,\xi^{\prime\prime\prime}(r_{\rm{infl}})}\bigg{)}^{2}.$
(15)

###  II.6 Comparison with standard method

Figure 5: Fractional uncertainties on the dip, peak, inflection point and LP
scales (square root of Eqs. 10, 9, 14, 12), as a function of the number of
eigen-modes that are included (symbols). Horizontal dashed lines show the
corresponding uncertainties estimated in the standard way (see main text).

To see how well this works, we first estimate the four scales (peak, dip, LP
and inflection point) in the standard way (e.g. Anselmi et al., 2016, 2019;
Parimbelli et al., 2021; Nikakhtar et al., 2021a; Paranjape and Sheth, 2023;
Lee et al., 2024): we made 100 mock realizations of the measurement (Eq. 5),
fitted a 7th order polynomial to each, and estimated the various scales by
differentiating the fit. The RMS scatter of each scale satisfies $\sigma_{\rm
dip}>\sigma_{\rm peak}>\sigma_{\rm infl}>\sigma_{\rm LP}$: the LP is the most
precise, followed by the inflection point, and etc. This is consistent with
previous work (e.g. Tables 4 to 7 of Lee et al. (2024)). The actual values are
shown as the horizontal dashed lines in Fig. 5. Previous work has shown that
it does not matter if one fits 7th order generalized Laguerre functions
instead (Nikakhtar et al., 2021a).

Now we turn to our eigen-mode-based estimates. Eqs. 10, 9, 14, 12 show that
these depend on the number of modes that are included in the relevant sums.
The symbols in Fig. 5 show how these estimates increase as more modes are
added. We previously argued that one should not include the higher order modes
(because one should not differentiate noisy data); these are the ones for
which the error estimate starts to diverge. The plateau at intermediate values
indicates that the error estimate is not very sensitive to exactly how many
modes are included, provided we have enough modes, and are not including the
ones which are dominated by shot-noise. (This plateau is not an artifact
produced by approximation Eq. 9; it is present even if we use Eq. 8.) In
effect, this plateau provides an objective measure of how many modes should be
included to accurately model the scatter between realizations, analogous to
how Fig. 1 provides an objective way of deciding which modes are most
important for a single realization. Indeed, previous work has noted that the
inflection point is slightly less robust than the LP: here, this is indicated
by the fact that it has a shorter plateau.

Figure 6: Same as Fig. 5, but with $\Sigma=0$ (no smearing), for both
$\xi_{0}(s)$ and the covariance matrix. Comparison with Fig. 5 (note the
difference in $y$-axis range) shows the potential gains in precision which
come from working with the reconstructed BAO signal.

It is reassuring that, not only do the plateau values reproduce the
qualitative trends shown by the standard method, they are within 80% of the
fitted RMS values in all cases. Some of this discrepancy arises from the
assumption that the scale which the standard method identifies as being the
peak, say $t_{\rm max}$ ($t$ for template), may differ slightly from $s_{\rm
max}$ (our eigen-mode based estimate of the peak scale). As a result, the RMS
of $t_{\rm max}-r_{\rm max}$ differs from the RMS of $s_{\rm max}-r_{\rm
max}$; evidently, this difference is small.

Recall that the standard method results do not depend on the functional form
that was fit to the binned TPCF. In effect, our analysis shows why: the low
order eigenvectors $\Lambda_{i}$ represent the covariance around any good fit
to the measurements which is not ‘fitting the noise.’

We noted that, because the BAO feature in $\xi$ becomes more smeared at late
times, $\xi^{\prime\prime}$ decreases, so we expect the uncertainty on the LP
scale to increase (c.f. Eqs. 9–12). Fig. 6 tests this: it shows the same
comparison as in Fig. 5, but now with $P(k)\propto P_{\rm Lin}(k)$ (i.e.
$\Sigma=0$, no smearing) for determining both $\xi_{0}(s)$ and $C_{ij}$ when
using Eq. 5 to produce 100 realizations of $\Xi(s)$. Setting $\Sigma\to 0$
changes the BAO feature in $\xi_{0}$ dramatically, and $C_{ij}$ less so: e.g.
the total variance is about 15% larger, and eigenvalues 3 to 7 contribute
considerably more to the total variance compared to when $\Sigma\neq 0$. The
change in $\xi_{0}$ means that a 7th order polynomial is no longer a good fit,
so, for the ‘standard’ analysis, we used a 9th order polynomial. Comparison of
the dashed lines here with those in the previous figure shows that the RMS is
decreased by a factor of about 2 to 3. This is especially true for the peak
which is most affected by the smearing.222It may seem surprising that, in
contrast to the evolved field, the peak and LP scales in the linear field are
measured with similar precision. This is mainly a signal-to-noise issue.
Recall that the original reason for working with the LP was not its precision,
but its robustness to evolution/smearing Anselmi et al. (2016). Notice that
this decrease is reproduced by our eigen-mode estimates. The onset of the
plateau is delayed from about mode 4 to about mode 7 or so, since now it is
the first 7 modes which contribute to most of the variance. This decrease
demonstrates the potential gains which come from using the reconstructed TPCF
rather than the smeared one to measure the BAO scales: for the LP, this
results in a precision of better than 0.4% as opposed to 0.8%. Our analysis
has allowed us to estimate this improvement without having to run simulations.

Some reconstruction methods move the observed biased tracers back to (an
estimate of) their initial positions Nikakhtar et al. (2022, 2023, 2024). The
TPCF is then measured using these reconstructed positions. In this case, the
number density is unchanged but $b\to b-1$: typically, the BAO feature is
sharper, but the amplitude is smaller Nikakhtar et al. (2022). Our methodology
allows an estimate of the precision of the LP distance scale in this
reconstructed signal as follows. Since $\bar{n}(b-1)^{2}P_{\rm Lin}$ is
smaller than $\bar{n}b^{2}P_{\rm Lin}$, the dotted curve in Fig. 1 suggests
that more modes will be needed before we converge to a plateau as in Fig. 5.
Although we do not show it here, we have checked that this is indeed the case.
In addition, the precision of the reconstructed feature is slightly less
precise, less constraining, than the original measurement. I.e.,
reconstruction yields no significant gain in precision (of course, it does
reduce the bias in the mean value, bringing the LP closer to its linear theory
value). To realize the potential for increased precision shown in Fig. 6, one
must combine reconstructed fields, as discussed in Nikakhtar et al. (2024).

We conclude that our methodology is an efficient way of determining accurate
uncertainties on the Linear Point estimate of the BAO distance scale. In
particular, since our analysis suggests that, for reasonable/realistic values
of the shot-noise, the relevant eigenvalues and eigenvectors are entirely
determined by the ratio of $P^{2}(k)$ and the survey volume $V_{s}$, they
scale as $(b\sigma_{8})^{2}/\sqrt{V_{s}}$. Our curves assumed
$(b,\sigma_{8},V_{s})=(1.97,0.84,5~{}h^{-3}{\rm Gpc}^{3})$, so for other
values, the fractional error is given by scaling the numbers in our Fig. 5 by
$(b_{\rm eff}/2.23)^{2}(\sigma_{8}/0.84)^{2}\sqrt{5/V_{s}}$. This is not quite
right, because it ignores the fact that the smearing also depends on
$\sigma_{8}$, but it is a useful first guess.

##  III Dependence of eigenvectors on binning

The previous section considered the structure of the covariance matrix of a
binned estimate of the TPCF. In that analysis, the original bin-size (and
shape) was fixed. How is the analysis modified if we change the binning?

In what follows, we first show that, in general, the eigenvectors of the
binned covariance matrix are not simply binned versions of the original
eigenvectors. Nevertheless, the first few eigenvectors are unchanged by the
variations in binning, permitted by the requirement that one be able to detect
the BAO feature in the first place.

###  III.1 Analytic analysis

Let $\mathbf{x}$ denote our list of bins, $\mathbf{y}$ the list of measured
bin amplitudes, $\mathbf{C}\equiv\langle\mathbf{yy}^{\rm T}\rangle$ the
covariance matrix of the measurements, and $\lambda_{i}$ and
$\mathbf{\Lambda}_{i}$ its eigenvalues and eigenvectors. With some abuse of
notation, let $\bm{\lambda}$ denote the diagonal matrix with the eigenvalues
along the diagonal. If $\mathbf{V}$ is a square matrix whose columns are the
eigenvectors $\mathbf{\Lambda}$, then

$\mathbf{C}=\mathbf{V}\bm{\lambda}\mathbf{V}^{-1}=\mathbf{V}\bm{\lambda}\mathbf{V}^{\rm
T},$ (16)

where the final equality follows because $\mathbf{C}$ is real and symmetric.

Suppose we bin so that $\mathbf{y}_{B}\equiv\mathbf{By}$, with $\mathbf{B}$ a
square matrix, each side having the same dimension as $\mathbf{y}$. E.g.,

$\mathbf{B}=\frac{1}{3}\begin{pmatrix}1&1&0&0&\cdots&0&0&0\\\
1&1&1&0&\cdots&0&0&0\\\ 0&1&1&1&\cdots&0&0&0\\\
\vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots&\vdots\\\
0&0&0&0&\cdots&1&1&1\\\ 0&0&0&0&\cdots&0&1&1\\\ \end{pmatrix}$ (17)

would correspond to averaging the $y$ values in the bins on either side of
each $x$-bin. Note that, when this is done, one typically works with a sparser
set of $\mathbf{x}$ values, so as to not ‘double-count.’ The analysis below is
more transparent when $\mathbf{x}_{\rm B}$ and $\mathbf{x}$, and hence
$\mathbf{y}_{\rm B}$ and $\mathbf{y}$, have the same length.

If $\mathbf{C}_{\rm B}$ denotes the covariance matrix of the binned
measurements, then

$\displaystyle\mathbf{C}_{\rm B}$
$\displaystyle\equiv\langle\mathbf{By\,(By)}^{\rm
T}\rangle=\mathbf{B}\mathbf{C}\mathbf{B}^{\rm T}$
$\displaystyle=\mathbf{B}\,(\mathbf{V}\bm{\lambda}\mathbf{V}^{-1})\,\mathbf{B}^{\rm
T}=(\mathbf{BV})\,\bm{\lambda}\,(\mathbf{V}^{\rm T}\mathbf{B}^{\rm T})$
$\displaystyle=(\mathbf{BV})\,\bm{\lambda}\,(\mathbf{BV})^{\rm T}.$ (18)

If $(\mathbf{BV})^{\rm T}=(\mathbf{BV})^{-1}$, the expression above would be
the eigenvalue decomposition of $\mathbf{C}_{\rm B}$, making it appear that
the eigenvalues of $\mathbf{C}_{\rm B}$ are the same as those of $\mathbf{C}$,
and the eigenvectors are simply those of $\mathbf{C}$, binned using
$\mathbf{B}$. At face value, this is sensible: if the eigenvectors were smooth
on scales smaller than the ‘bin width’ then they will be unchanged by –
essentially invariant to – the binning.

However, notice that

$\mathbf{BV}\,(\mathbf{BV})^{\rm T}=\mathbf{B}\,(\mathbf{VV}^{\rm
T})\,\mathbf{B}^{\rm T}=\mathbf{B}\,\mathbf{B}^{\rm T}$ (19)

is not diagonal (i.e., although $\mathbf{V}$ is an orthonormal basis,
$\mathbf{V}^{\rm T}=\mathbf{V}^{-1}$, the same is not true for $\mathbf{BV}$).
This means that we should not think of $\mathbf{BV}$ as being the
eigenvectors. If we use $\mathbf{V}_{\rm B}$ to denote the eigenvectors of
$\mathbf{C}_{\rm B}$,

$\mathbf{C}_{\rm B}=\mathbf{V}_{\rm B}\bm{\lambda}_{\rm B}\mathbf{V}_{\rm
B}^{-1},$ (20)

then it is natural to ask: How different are the vectors which make up
$\mathbf{V}_{\rm B}$ from those of $\mathbf{BV}$?

###  III.2 Numerical analysis

Heuristically, we expect that if the binning remains smaller than the typical
size of features in the eigenvectors, then they will be unchanged by binning.
This should be particularly true for the primary ‘cosmic variance’ dominated
eigenvectors; the shot-noise dominated eigenvectors oscillate more, but we
argued that they are not interesting anyway. Therefore, we expect the
estimates of the BAO distance scale and their uncertainties should not depend
on how the TPCF was binned, provided this binning is not wider than the BAO
feature itself. (If the bins are too wide, they will not provide a good
description of the BAO feature anyway.)

Figure 7: Dependence of eigenvalues on bin-size and spacing. Wider bins have
smaller total variance, so we have normalized each eigenvalue to show its
fractional contribution to the total. The first 4 normalized eigenvalues are
very similar between all 3 curves, while the first 4 values for non-
overlapping $6~{}h^{-1}$Mpc bins are essentially the same as the overlapping
$6~{}h^{-1}$Mpc bins. Figure 8: First 9 eigenvectors
$\sqrt{\lambda_{i}}\Lambda_{i}(s)$ at $z=0.5057$, for bin-sizes
$2~{}h^{-1}$Mpc and $6~{}h^{-1}$Mpc with the same bin centers. Larger symbols
show the eigenvectors for the $6~{}h^{-1}$Mpc bins that do not overlap (so
there are $3\times$ fewer bin centers, different from the overlapping bins).

Figs. 7 and 8 show the result of two explicit tests. The first increases the
bin-size, but keeps the bin centers, and hence the number of bins, the same.
(As a result, neighboring bins are more correlated, but this just means that
$\mathbf{C}_{\rm B}$, which has the same dimension as in the previous section,
is less diagonal than the original $\mathbf{C}$.) This corresponds to the
$\epsilon\neq 1$ case mentioned previously (Section II.2). Fig. 7 shows the
eigenvalues when we increase $\Delta r$ by factors of 2 and 3 ($\epsilon=1/2$
and $\epsilon=1/3$), respectively. The first few eigenvalues, which dominate
the total variance, are indistinguishable from the original ones, but the more
shot-noise dominated modes are affected. In particular, for wider bins the
shot-noise is smaller, so these shot-noise dominated modes contribute less to
the total variance. To remove the fact that the total variance is reduced, we
normalized each set of eigenvalues by their total: this shows explicitly that,
as the bin-size increases, the shot-noise dominated modes contribute a smaller
fraction of the total variance.

The larger symbols show the eigenvalues when $\Delta r\to 3\Delta r$ but the
bins do not overlap ($\epsilon=1$). In this case, there are $3\times$ fewer
eigenvalues, so the total variance is obviously different. Nevertheless, the
fractional contribution of the first 10 modes to the variance is similar to
that for the overlapping bins (of the same width). Clearly, for the cosmic-
variance dominated modes which contribute most to the total variance, the
binning does not matter.

Fig. 8 shows that this is also true for the eigenvectors. The big symbols show
the case in which $\Delta r\to 3\Delta r$ but now the bins do not overlap (so
$\epsilon=1$ even for this larger bin-size; for clarity, we do not show the
intermediate case where $\Delta r\to 2\Delta r$). Again, the eigenvectors
which dominate the total variance (the first $\sim 5$) are indistinguishable
from the original ones. This is slightly non-trivial since now $C_{\rm B}$ is
$10\times 10$ rather than $30\times 30$, but the leading eigenvectors are
unchanged. Hence, the uncertainties on distance scale estimates provided in
the previous section will be unchanged: they do not depend on the binning, at
least for fiducial values of the shot-noise.

More generally, if the convolution kernel which defines the binning does not
erase features in the original (cosmic variance dominated) eigenvectors, then
these eigenvectors will not depend on the exact bin shape. E.g., this will
certainly be true if the off-diagonal entries in Eq. 17 are less than unity.
Similarly, counting pairs in, e.g., Gaussian-like bins rather than in
rectangles will not change our conclusions.

##  IV Discussion and Conclusions

We presented an eigen-decomposition of the Gauss-Poisson approximation to the
covariance matrix of the two-point correlation function (Eq. 3) and assessed
the importance of the power spectrum-dominated modes that trace cosmic
variance as opposed to the modes which are dominated by shot-noise. For a
fiducial cosmology and noise-levels that are consistent with current and next-
generation surveys, the cosmic variance eigen-modes account for most of the
total variance of the TPCF (Fig. 1). They are also smoother than the shot-
noise dominated modes (Fig. 2), so they are insensitive to the modest changes
in binning that are allowed if one wishes to resolve the BAO feature in the
TPCF (Figs. 7 and 8).

We argued that, as a result, the cosmic variance eigen-modes alone should
provide a good description of the correlated residuals which result from
fitting smooth functional forms to the measured TPCF. We provided a simple
(Eq. 12) but accurate (Fig. 5) approximation for the uncertainty on the Linear
Point estimate of the BAO distance scale which explains why the uncertainty is
greater in the evolved field than in linear theory; allows one to quantify the
gains from working with the reconstructed signal; and does not depend on the
functional form fitted to the TPCF or the binning used. It also provides
insight into why the LP is more robust than the inflection point, and why both
are more precise distance indicators than the peak or dip scales. Perhaps most
importantly, our approximation allows one to quantify the precision of the BAO
distance scale estimate without having to generate a large ensemble of mock
catalogs. Therefore, it should be useful for estimating the gains in precision
which come from making measurements in the reconstructed field (which are
often quoted), after marginalizing over the unknown cosmological model (which
is often ignored).

The plateaus that are apparent in Figs. 5 and 6 suggest that it should be
possible to write down a prescription for determining the optimal number of
eigen-modes which should be used in cosmological analyses. In future work, we
hope to combine our eigen-mode analysis with the Bayesian framework of
Paranjape and Sheth (2022).

###### Acknowledgements.

JL was supported by DOE grant DE-FOA-0002424 and NSF grant AST-2108094. FN
gratefully acknowledges support from the Yale Center for Astronomy and
Astrophysics Prize Postdoctoral Fellowship. AP and RKS are grateful to the
ICTP, Trieste for its hospitality in summer 2024, and RKS is grateful to the
EAIFR, Kigali for its hospitality when this work was completed. The research
of AP is supported by the Associates Scheme of ICTP, Trieste.

## References

  * Eisenstein et al. (2005) D. J. Eisenstein, I. Zehavi, D. W. Hogg, R. Scoccimarro, M. R. Blanton, R. C. Nichol, R. Scranton, H.-J. Seo, M. Tegmark, Z. Zheng, et al., The Astrophysical Journal 633, 560 (2005).
  * Anselmi et al. (2018a) S. Anselmi, G. D. Starkman, P.-S. Corasaniti, R. K. Sheth, and I. Zehavi, Physical Review Letters 121, 021302 (2018a).
  * Alam et al. (2021) S. Alam, M. Aubert, S. Avila, C. Balland, J. E. Bautista, M. A. Bershady, D. Bizyaev, M. R. Blanton, A. S. Bolton, J. Bovy, et al., Physical Review D 103, 083533 (2021).
  * Abbott et al. (2024) T. Abbott, M. Adamow, M. Aguena, S. Allam, O. Alves, A. Amon, F. Andrade-Oliveira, J. Asorey, S. Avila, D. Bacon, et al., arXiv preprint arXiv:2402.10696 (2024).
  * Adame et al. (2024a) A. Adame, J. Aguilar, S. Ahlen, S. Alam, D. Alexander, M. Alvarez, O. Alves, A. Anand, U. Andrade, E. Armengaud, et al., arXiv preprint arXiv:2404.03002 (2024a).
  * Smith et al. (2008a) R. E. Smith, R. Scoccimarro, and R. K. Sheth, Physical Review D 77, 043525 (2008a).
  * Grieb et al. (2016) J. N. Grieb, A. G. Sánchez, S. Salazar-Albornoz, and C. Dalla Vecchia, Monthly Notices of the Royal Astronomical Society 457, 1577 (2016).
  * Anselmi et al. (2016) S. Anselmi, G. D. Starkman, and R. K. Sheth, Monthly Notices of the Royal Astronomical Society 455, 2474 (2016).
  * Anselmi et al. (2019) S. Anselmi, P.-S. Corasaniti, A. G. Sanchez, G. D. Starkman, R. K. Sheth, and I. Zehavi, Physical Review D 99, 123515 (2019).
  * Parimbelli et al. (2021) G. Parimbelli, S. Anselmi, M. Viel, C. Carbone, F. Villaescusa-Navarro, P. Corasaniti, Y. Rasera, R. Sheth, G. Starkman, and I. Zehavi, Journal of Cosmology and Astroparticle Physics 2021, 009 (2021).
  * O’Dwyer et al. (2020) M. O’Dwyer, S. Anselmi, G. D. Starkman, P.-S. Corasaniti, R. K. Sheth, and I. Zehavi, Physical Review D 101, 083517 (2020).
  * Nikakhtar et al. (2021a) F. Nikakhtar, R. K. Sheth, and I. Zehavi, Physical Review D 104, 043530 (2021a).
  * Lee et al. (2024) J. J. Lee, B. Fiorini, F. Nikakhtar, and R. K. Sheth, arXiv e-prints arXiv:2406.09379 (2024), eprint 2406.09379.
  * Paranjape and Sheth (2022) A. Paranjape and R. K. Sheth, Mon. Not. R. Astron. Soc. 517, 4696 (2022), eprint 2209.00668.
  * He et al. (2023) M. He, C. Zhao, and H. Shan, Monthly Notices of the Royal Astronomical Society 525, 1746 (2023).
  * Vargas-Magaña et al. (2014) M. Vargas-Magaña, S. Ho, X. Xu, A. G. Sánchez, R. O’Connell, D. J. Eisenstein, A. J. Cuesta, W. J. Percival, A. J. Ross, E. Aubourg, et al., Mon. Not. R. Astron. Soc. 445, 2 (2014).
  * Anselmi et al. (2018b) S. Anselmi, P.-S. Corasaniti, G. D. Starkman, R. K. Sheth, and I. Zehavi, Physical Review D 98, 023527 (2018b).
  * Nikakhtar et al. (2021b) F. Nikakhtar, R. K. Sheth, and I. Zehavi, Physical Review D 104, 063504 (2021b).
  * Alam et al. (2015) S. Alam, F. D. Albareti, C. A. Prieto, F. Anders, S. F. Anderson, T. Anderton, B. H. Andrews, E. Armengaud, É. Aubourg, S. Bailey, et al., The Astrophysical Journal Supplement Series 219, 12 (2015).
  * Adame et al. (2024b) A. Adame, J. Aguilar, S. Ahlen, S. Alam, D. Alexander, M. Alvarez, O. Alves, A. Anand, U. Andrade, E. Armengaud, et al., arXiv preprint arXiv:2404.03000 (2024b).
  * Kaiser (1987) N. Kaiser, Monthly Notices of the Royal Astronomical Society 227, 1 (1987).
  * Crocce and Scoccimarro (2008) M. Crocce and R. Scoccimarro, Phys. Rev. D 77, 023533 (2008), eprint 0704.2783.
  * Paranjape and Sheth (2023) A. Paranjape and R. K. Sheth, arXiv preprint arXiv:2304.09198 (2023).
  * Smith et al. (2008b) R. E. Smith, R. Scoccimarro, and R. K. Sheth, Physical Review D 77, 043525 (2008b).
  * Nikakhtar et al. (2022) F. Nikakhtar, R. K. Sheth, B. Lévy, and R. Mohayaee, Phys. Rev. Lett. 129, 251101 (2022), eprint 2203.01868.
  * Nikakhtar et al. (2023) F. Nikakhtar, N. Padmanabhan, B. Lévy, R. K. Sheth, and R. Mohayaee, Phys. Rev. D 108, 083534 (2023), eprint 2307.03671.
  * Nikakhtar et al. (2024) F. Nikakhtar, R. K. Sheth, N. Padmanabhan, B. Lévy, and R. Mohayaee, Phys. Rev. D 109, 123512 (2024), eprint 2403.11951.