# Baryon Acoustic Oscillations analyses with Density-Split Statistics

Tengpeng Xu Key Laboratory for Computational Astrophysics, National
Astronomical Observatories, Chinese Academy of Sciences
Beijing 100101, China  College of Astronomy and Space Sciences, University of
Chinese Academy of Sciences
Beijing, 100049, China  Yan-Chuan Cai Institute for Astronomy, University of
Edinburgh, Royal Observatory
Blackford Hill, Edinburgh, EH9 3HJ, UK  Yun Chen Key Laboratory for
Computational Astrophysics, National Astronomical Observatories, Chinese
Academy of Sciences
Beijing 100101, China  College of Astronomy and Space Sciences, University of
Chinese Academy of Sciences
Beijing, 100049, China  Mark Neyrinck Department of Physics and Astronomy,
University of Denver, Denver, CO 80210, USA  Blue Marble Space Institute of
Science, Seattle, WA 98104, USA  Liang Gao Key Laboratory for Computational
Astrophysics, National Astronomical Observatories, Chinese Academy of Sciences
Beijing 100101, China  Qiao Wang Key Laboratory for Computational
Astrophysics, National Astronomical Observatories, Chinese Academy of Sciences
Beijing 100101, China  College of Astronomy and Space Sciences, University of
Chinese Academy of Sciences
Beijing, 100049, China Yan-Chuan Cai, Yun Chen<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Accurate modeling for the evolution of the Baryon Acoustic Oscillations (BAO)
is essential for using it as a standard ruler to probe cosmology. We explore
the non-linearity of the BAO in different environments using the density-split
statistics and compare them to the case of the conventional two-point
correlation function (2PCF). We detect density-dependent shifts for the
position of the BAO with respect to its linear version using halos from N-body
simulations. Around low/high-densities, the scale of the BAO expands/contracts
due to non-linear peculiar velocities. As the simulation evolves from redshift
1 to 0, the difference in the magnitude of the shifts between high- and low-
density regions increases from the sub-percent to the percent level. In
contrast, the scale of the BAO does not evolve in the total 2PCF in the same
redshift range. The width of the BAO around high density regions increases as
the universe evolves, similar to the known broadening of the BAO in the 2PCF
due to non-linear evolution. In contrast, the width is smaller and stable for
low density regions. We discuss possible implications for the reconstructions
of the BAO in light of our results.

Large-scale structure of the universe(902) — Baryon acoustic oscillations(138)
— N-body simulations(1083)

††software:  densitysplit (Paillas et al., 2021), pycorr (Sinha & Garrison,
2019, 2020), Nbodykit (Hand et al., 2018), emcee (Foreman-Mackey et al.,
2013), getdist (Lewis, 2019)

##  1 Introduction

Baryon Acoustic Oscillations (BAO) represents the sound waves in the
primordial baryon-photon plasma before recombination due to the interactions
between gravity, provided by the total matter, and the pressure, provided by
baryons. The largest scale at which sound waves can propagate by the time of
recombination defines the sound horizon. After recombination, the scale of the
sound horizon is frozen into the baryon density fluctuations (and hence the
matter density fluctuations) in comoving units, denoted as $r_{\rm d}\sim 147$
Mpc$\sim 100$ Mpc/$h$. $r_{\rm d}$ represents the sound horizon at the drag
epoch ($z_{\rm d}$), when baryons are released from the Compton drag of the
photons, determined by the ratio of baryon to matter densities (Eisenstein &
Hu, 1998; Hu & Sugiyama, 1996). The sound horizon has been measured precisely
by WMAP and Planck through the temperature fluctuations of the cosmic
microwave background (e.g. Bennett et al., 2013; Planck Collaboration et al.,
2020). These acoustic oscillations, imprinted in the matter density with a
comoving scale of $r_{\rm d}$, result in a peak in the two-point correlation
function of matter, serving as one of the most precise standard rulers in
cosmology for measuring cosmological distances. In the last few decades,
generations of surveys such as 2dFGRS, 6dFGS, SDSS, BOSS, eBOSS and DESI had
used it to set constraints on cosmological parameters (Cole et al., 2005;
Beutler et al., 2011; Eisenstein et al., 2005; Ross et al., 2015; Alam et al.,
2017, 2021; Zarrouk et al., 2018; DESI Collaboration et al., 2024).

However, the position of the BAO is not fixed after recombination in comoving
coordinates, which we assume throughout. Non-linear growth of structure is
expected to cause the location of the BAO to shift, typically at the sub-
percent level; the width of the BAO is broadened relative to its initial state
(e.g. Crocce & Scoccimarro, 2008; Seo et al., 2008, 2010; Orban & Weinberg,
2011). The broadening of the peak of the BAO can be thought of as the
consequence of peculiar velocities around high and low density regions pulling
the BAO in opposite directions, usually in a non-symmetric way. This weakens
the significance of the BAO peak and complicates the modeling. Significant
effort has been made to reconstruct the BAO to reinstall the information from
its initial conditions. A common practice is to perform reconstructions by
using the observed late-time positions of galaxies, estimate their velocities
and move the galaxies back to their initial conditions (e.g. Eisenstein et
al., 2007; Noh et al., 2009; Padmanabhan et al., 2009; Seo et al., 2010;
White, 2015; Schmittfull et al., 2015; Wang et al., 2017; Yu et al., 2017; Mao
et al., 2021; von Hausegger et al., 2022; Nikakhtar et al., 2022, 2023; Chan
et al., 2024). Indeed, reconstructions applied to cubic-box simulations were
shown to substantially reduce the broadening and nonlinear shifts of the BAO,
but this may not always be the case in real observations with complicated
survey masks (e.g. Anderson et al., 2012), or when the shot-noise level of the
galaxy sample is high (White, 2010; Wang & Pen, 2019; Seo et al., 2022). In
principle, if we can track the non-linear evolution of the BAO precisely,
forward modeling can also be applied to extract cosmological information (e.g.
Babić et al., 2022). In addition, introducing environmental-dependence weights
for the reconstruction by up-weighting under-dense environments has been
demonstrated to help sharpening the peak of the BAO (Achitouv & Blake, 2015).
It has also been shown that reconstructions sharpen the BAO around voids more
than for the 2PCF of mock galaxies (Zhao et al., 2020). All these indicate
potential benefits for understanding the environmental dependence for the non-
linear evolution of the BAO, which is the focus of our study.

In light of the above, we take a slightly different approach to analyze the
BAO. We will identify density fluctuations in the late-time Universe around
which the scale of the BAO may have changed relative to its initial position,
and analyze them separately. This was first proposed by Neyrinck et al. (2018)
using a similar concept, which they call sliced correlation functions.

In this paper, we will take a step forward to perform BAO analyses with halos
from the Quijote simulations (Villaescusa-Navarro et al., 2020). We will use
the density-split clustering (DS), or conditional correlation functions for
the analysis (Paillas et al., 2021, 2023). We will split the two-point
correlation function into a series of cross-correlation functions to be
conditioned on local densities of different depth. By doing this, we aim to
identify the possible expansion and contraction of the BAO scale around
different density environments, with the perspective of tracking the
broadening of the BAO. In this method, the 2PCF can be seen as a composition
of the whole spectrum of cross-correlation functions. In addition, recent
studies with DS-clustering have suggested significant gain of cosmological
information over two-point correlation functions (Paillas et al., 2023, 2024).
We would like to explore if BAO analyses with DS-clustering have any
significant impact in the error budget compared to two-point statistics. We
will address both the shifts and the broadening of the BAO together by
performing error analyses.

The paper is organized as follows. In Section 2 we introduce the Quijote
simulations, the halo catalogs used in this work, and describing the
measurements of the density-split clustering. We present our models and their
parameters for the correlation functions in Section 3. Our BAO analyses are
presented in Section 4. We provide discussions and the main conclusions in
Section 5.

##  2 DENSITY-SPLIT clustering measurements from SIMULATIONS

###  2.1 The Quijote simulation

The Quijote project (Villaescusa-Navarro et al., 2020) comprises a series of
N-body simulations aimed at assessing the information content in cosmological
observables, spanning a wide range of values around its fiducial cosmology,
which is consistent with the latest Planck constraints. The parameters of the
fiducial cosmology are set as follows: $\Omega_{\rm m}=0.3175$, $\Omega_{\rm
b}=0.049$, $h=0.6711$, $n_{s}=0.9624$, $\sigma_{8}=0.834$, $M_{\nu}(\rm
eV)=0$, and $w=-1$. There are 15,000 realizations of the fiducial cosmology
encompassed in the Quijote simulations, providing a substantial dataset for
the computation of covariance matrices. Dark matter particles and halos
catalogs obtained with a Friends-of-Friends algorithm (Davis et al., 1985) are
provided. In our work, we utilized the dark matter halo catalogs with a
minimal halo mass of $1.31\times 10^{13}~{}M_{\odot}/h$ from 1000 realizations
of the simulations. Table. 1 summarizes some properties of the catalogs. We
focus on the data at redshifts $z=0$, 0.5, and 1 where most galaxy redshift
surveys can cover.

snapshots ID | 8000 | 8001 | 8002 | 8003 | 8004 | … | 8999 | average number | standard deviation
---|---|---|---|---|---|---|---|---|---
$z=1.0$ | 195633 | 195604 | 195236 | 195747 | 195096 | … | 195805 | 195575.9 | 362.9
$z=0.5$ | 308580 | 310607 | 309485 | 309694 | 309832 | … | 309929 | 309426.1 | 476.6
$z=0.0$ | 406728 | 407631 | 406226 | 406035 | 407212 | … | 406990 | 406764.2 | 642.8
Table 1: The samples of dark matter halos from the Quijote simulation used for
our analyses. They have the minimal halo mass $M_{\rm min}=1.31\times
10^{13}~{}M_{\odot}/h$ at $z=$ 0, 0.5 and 1.

###  2.2 Measurements of clustering

We follow the algorithm presented in Paillas et al. (2021, 2023) to split the
halo field into density bins. We apply spherical top-hat smoothing with a
radius of $R=20h^{-1}\rm Mpc$ to the halo number density field, and measure
the smoothed density at random locations (called DS centers, for centers of
top-hat smoothing spheres) of the simulation box. Note that a DS center does
not usually correspond to the location of a halo. The number of randoms in
each box are set to be $M$ times the number of halos, $N_{\rm random}=M\times
N_{\rm halos}$, where $M$ is the number of density bins we are going to use
for the split. This is to insure fair comparisons with the 2PCF (Paillas et
al., 2023). This is equivalent to counting the number of halos around the
random positions to measure the halo number density fluctuations within $R$:

$\Delta(R)=\frac{\rho(<R)}{\bar{\rho}}-1,$ (1)

where $\bar{\rho}$ is the mean halo density of the entire simulated box. The
probability density function (PDF) of the density contrasts for a single
realization at $z=0$ is shown in Fig. 1.

Figure 1: The normalized PDF of the halo number density contrast smoothed with
a spherical top-hat window of the radius $R=20~{}$Mpc/$h$. Regions of
increasing local densities, corresponding to ${\rm DS}_{1}$ to ${\rm DS}_{5}$
are shown in different colors. The auto-correlations of these DS centers, and
their cross-correlations with halos are presented in Fig. 2. These are the
basis for the density-split clustering analyses described in the paper. The
measurement was conducted for one box of the Quijote simulation at $z=0$.
Figure 2: Measurements from the Quijote simulations for the cross-correlation
functions between halos and density-split centers (Halo $\times$ DS, left
panels) and the auto-correlation functions of density-split centers (DS
$\times$ DS, right panels). The top, middle and bottom panels show the results
from redshifts $z=$ 0, 0.5 and 1, respectively. Solid lines are the averages
from 1000 realizations of Quijote simulations for the fiducial cosmology, with
their errors shown in shaded regions. Different colors indicate the auto- and
cross-correlations for different ${\rm DS}_{i}$ (i=1, 2, 3, 4, 5).  Figure 3:
Left: A 2D slice of halo number counts from one realization of the Quijote
simulations. The thickness of the slice is 50 Mpc/$h$, with each pixel of the
size of 5 Mpc/$h$ $\times$ 5 Mpc/$h$. Middle: Number counts of DS1 centers
(the most under-dense regions) selected from filtering the 3D halo number
density field with a spherical top-hat filter of the radius of $20$Mpc/$h$
(see Section 2.2 and Fig. 1 for more details). Right: Similar to the middle
panel but showing the number counts of DS5 centers (the most over-dense
regions). It can be observed that the distribution of DS1 mostly corresponds
to the void regions in the halo distribution, whereas the DS5 distribution is
largely complementary to that of DS1, occupying high-density regions.

Next, we divide the DS centers equally into $M$-bins and by default, five
quintiles, naming them from $\rm DS_{1}$ to $\rm DS_{5}$. Hence $\rm DS_{1}$
corresponds to under-density regions, while $\rm DS_{5}$ corresponds to over-
density regions. We also repeat our analyses with $M=10$ with results
presented in Appendix B.

Fig. 3 display the distribution of halos, DS1 and DS5 centers in a two-
dimensional slice with a thickness of $50$ Mpc$/h$ at redshift $z=0$
(corresponding to the three subplots from left to right, respectively). We can
visually observe that the distribution of DS1 mostly corresponds to the void
regions in the halo distribution, whereas the DS5 distribution is largely
complementary to that of DS1.

With these DS centers, we utilize the pycorr software package (Sinha &
Garrison, 2019, 2020) to calculate the cross-correlation functions between
$\rm DS_{1}$-$\rm DS_{5}$ centers and dark matter halos for the 1000
Realizations at $z=0,0.5,1.0$. We also compute the auto-correlation functions
of $\rm DS_{1}$-$\rm DS_{5}$ centers. In the following text, we denote the
cross-correlation functions as “Halo $\times$ DS” or $\xi_{i}^{\rm qh}$, and
the auto-correlation functions of ${\rm DS}_{i}$ as “DS $\times$ DS” or
$\xi_{i}^{\rm qq}$, where “q” and “h” are abbreviations for “quintile” and
“halo”, respectively. The measurements of these correlation functions are
shown in Fig. 2.

##  3 THEORETICAL MODELS OF CORRELATION FUNCTIONS

For the model, we first utilize the Nbodykit package (Hand et al., 2018) to
generate the linear power spectrum and the corresponding correlation functions
at $z=0,0.5$ and 1.0 based on the same fiducial cosmology of the Quijote
simulations. We use the sideband algorithm (Kirkby et al., 2013) to decompose
the correlation function into two parts: one that exclusively contains the BAO
peak ($\xi_{\rm peak}(r)$) and another that represents the smooth portion
without the BAO peak ($\xi_{\rm smooth}(r)$). We then Fourier transform these
components to obtain the corresponding power spectra, $P_{\rm peak}(k)$, and
$P_{\rm smooth}(k)$. Subsequently, we can calculate the quasi-linear power
spectrum as:

$P_{\rm QL}(k;\sigma)=P_{\rm smooth}(k)+\exp(-\frac{k^{2}\sigma^{2}}{2})P_{\rm
peak}(k),$ (2)

where $\sigma$ is the BAO peak’s non-linear damping term (referred to as the
BAO width parameter in the following text). Therefore, we can express the
tracer biased power spectrum as given by

$\hat{P}(k;b_{1},b_{2},\sigma)=b_{1}b_{2}P_{\rm QL}(k;\sigma),$ (3)

where $b_{1}$ and $b_{2}$ represent the biases of the two tracers of interest
e.g., $b_{\rm halo}$ and $b_{\rm DSi}$ in the context of this study. We obtain
the corresponding two-point correlation functions
$\hat{\xi}(r;b_{1},b_{2},\sigma)$ through Fourier transformation of the model.
Finally, we introduce the scale dilation parameter $\alpha$ for the
distortions to the correlations due to assuming a cosmology which is different
from the fiducial one (Alcock & Paczynski, 1979; Ballinger et al., 1996). In
general, the scale dilation parameters for the line-of-sight distance and the
transverse distance are:

$\displaystyle\alpha_{\parallel}$ $\displaystyle=\frac{[D_{\rm H}(z)/r_{\rm
d}]}{[D_{\rm H}(z)/r_{\rm d}]_{\rm fid}}$ (4) $\displaystyle\alpha_{\perp}$
$\displaystyle=\frac{[D_{\rm M}(z)/r_{\rm d}]}{[D_{\rm M}(z)/r_{\rm d}]_{\rm
fid}},$ (5)

where $D_{\rm H}$ and $D_{\rm M}$ are the Hubble distance and the comoving
angular diameter distance at $z$ respectively. These parameters is expected to
provide the main cosmological constraints using the BAO. In this work, we
focus on the possible shifts of the BAO positions, and so we do not include
redshift space distortions in our analyses. Therefore, the correlation
functions are isotropic, so do the AP distortions, meaning that
$\alpha_{\parallel}=\alpha_{\perp}$, where we denote as $\alpha$ in the
following text. The theoretical model of the correlation function is given as:

$\xi_{\rm
th}(r;b_{1},b_{2},\alpha,\sigma)=\hat{\xi}(\alpha*r;b_{1},b_{2},\sigma).$ (6)

Following the same notations introduced in Section 2.2, ${\xi}_{i}^{\rm
qh}(r;b_{{\rm halo|DS}_{i}},b_{{\rm DS}_{i}},\alpha^{\rm qh}_{i},\sigma^{\rm
qh}_{i})$ represents the modeled cross-correlation function between dark
matter halos and ${\rm DS}_{i}$ centers, where $b_{{\rm halo|DS}_{i}},b_{{\rm
DS}_{i}},\alpha^{\rm qh}_{i},\sigma^{\rm qh}_{i}$ are free parameters, and
$i=1,2...5$. The same applies to the auto-correlation function of ${\rm
DS}_{i}$ centers, ${\xi}_{i}^{\rm qq}(r;b_{{\rm DS}_{i}},\alpha^{\rm
qq}_{i},\sigma^{\rm qq}_{i})$, and the auto-correlation function of dark
matter halos, $\xi^{\rm hh}(r;b_{\rm halo},\alpha^{\rm hh},\sigma^{\rm hh})$.
It is noteworthy that, in this study, the halo bias in the cross-correlation
functions between dark matter halos and ${\rm DS}_{i}$ centers (denoted as
$b_{{\rm halo|DS}_{i}}$) and the halo bias in the auto-correlation functions
of dark matter halos (denoted as $b_{\rm halo}$) are considered as distinct
free parameters. This allows us to explore potential disparities between the
halo bias around ${\rm DS}_{i}$ centers and the global halo bias.

The $\alpha$ parameter is a correction factor for the true cosmic distance
relative to the fiducial cosmic distance. It also reflects deviations from the
real position of the BAO relative to the reference cosmology. We treat
$\alpha$ of different correlation functions to be independent free parameters,
with $\alpha^{\rm hh}$ for the halo clustering, $\alpha^{\rm qh}_{i}$ for DS-
Halo cross-correlations, and $\alpha^{\rm qq}_{i}$ for DS-DS auto-
correlations. From the definition, when $\alpha>1$, the measured correlation
function exhibits compression along the $r$ axis, and so the BAO scale is
smaller than in the fiducial cosmology. Conversely, when $\alpha<1$, the
opposite occurs. Only when $\alpha=1$, there is no deviation between the
measured BAO peak position and that of the reference universe.

In our analyses, the fiducial cosmology assumed is also the true cosmology,
and so any deviation of $\alpha$ from 1 must come from non-linear evolution.
We will call it the BAO shift parameter in the rest of the paper. The level of
deviation for $\alpha$ for the 2PCF or power spectrum has been known to be at
the sub-percent level at the late-time Universe (e.g. Crocce & Scoccimarro,
2008; Seo et al., 2008, 2010; Orban & Weinberg, 2011). With density-split
clustering, the BAO scale around different background densities may shift
differently, and the possible deviation of $\alpha$ from 1 may become more
obvious. In this work, we focus on the two most extreme densities – ${\rm
DS}_{5}$ and ${\rm DS}_{1}$, corresponding to the highest and lowest density
fields. Guided by Neyrinck et al. (2018), these largest density perturbations
may cause the scale of the BAO to have the most significant shift. We will
perform fitting with the model presented above with the measurements presented
in section 2 to investigate this.

##  4 MCMC fitting for BAO positions

We employ the Markov Chain Monte Carlo (MCMC) method to constrain parameters
in the correlation function model, utilizing the emcee software package
(Foreman-Mackey et al., 2013). For any of the correlation functions, $\xi^{\rm
hh}$, $\xi^{\rm qh}_{i}$ or $\xi^{\rm qq}_{i}$, the log-likelihood can be
computed by comparing the theoretical data vector $\bm{\xi_{\rm th}}$ with the
measured data vector $\bm{\xi_{\rm sim}}$ obtained from simulations:

$\displaystyle\ln\mathcal{L}$ $\displaystyle=-\frac{1}{2}\chi^{2},$ (7)
$\displaystyle\chi^{2}$ $\displaystyle=(\bm{\xi_{\rm sim}}-\bm{\xi_{\rm
th}})\mathbf{C}^{-1}(\bm{\xi_{\rm sim}}-\bm{\xi_{\rm th}})^{\mathbf{T}}.$ (8)

For the data vector $\bm{\xi_{\rm sim}}$, we take the average over the
measurements from the 1000 simulations. The covariance matrix $\mathbf{C}$ is
also computed using the 1000 simulation boxes of the fiducial cosmology:

$\mathbf{C}=\frac{1}{N_{\rm sim}-1}\sum^{N_{\rm sim}}_{k=1}(\bm{\xi_{{\rm
sim},k}}-\bm{\bar{\xi}_{\rm sim}})^{\mathbf{T}}(\bm{\xi_{{\rm
sim},k}}-\bm{\bar{\xi}_{\rm sim}}),$ (9)

where $N_{\rm sim}=1000$. Since the possible shift of the BAO peak position is
expected to be small, we would like to reduce random errors as far as
possible. We scaled down the covariance matrix by a factor of 1/100, which
corresponds to reducing the error bars of each data point to one-tenth of
their original sizes. This is equivalent to increasing the volume of the
simulated box by a factor of 100, namely, the new box has a side length of
approximately $1~{}{\rm Gpc/}h\times 100^{1/3}\approx 4.6~{}{\rm Gpc/}h$. For
the fitting, we use the range of scales from 80 Mpc/$h$ to 180 Mpc/$h$ to
encompasses the co-moving BAO scale ($r\sim 100~{}$Mpc/$h$), and make linear
binning for the data with the bin-width of 4 Mpc/$h$. We have compiled the
parameter constraint results for all the correlation functions described above
in Table 2, and present the full posterior distributions for all parameters in
the Appendix A. Next, we will proceed with a detailed analysis of each
parameter.

Figure 4: Comparing the correlation coefficients between the DS-DS and Halo-DS
clustering at $z=0.0$, 0.5 and 1.0, where the dotted lines represent DS1 (low
densities) and the solid lines represent DS5 (high densities). These are
derived from the full covariance matrix of [$\bf\xi^{qq},\xi^{qh}$]. The
coefficients for DS5 are close to unity, suggesting that there is little gain
of information when combining DS-DS with Halo-DS. The coefficients for DS1 are
typically smaller than those for DS5, suggesting more room for improvement
when combining DS-DS with Halo-DS for under-dense regions.

###  4.1 BAO constraints from halo-halo correlations

The fitted results for the halo-halo correlation functions (2PCF) at $z=0,0.5$
and 1.0 are shown in Fig. 5 (A). We can see that the model agrees very well
with the measurements from simulations across the range of scales shown. Fig.
6 (A, B, C) presents the constrained for $b_{\rm halo},\alpha^{\rm hh}$ and
$\sigma^{\rm hh}$. We can see that at these three redshifts, the central
values of the alpha parameter are approximately 0.36% greater than 1, and
their differences from 1 all exceed the $1\sigma$ level. The fact that
$\alpha^{\rm hh}>1$ suggests a contraction of the BAO scale relative to its
linear version, consistent with previous studies (e.g. Crocce & Scoccimarro,
2008; Seo et al., 2008, 2010; Orban & Weinberg, 2011). This is expected as
halos with the bias greater than 1 occupied over-dense regions. Overall, the
co-moving BAO scales are consistent with no significant evolution within the
redshift range of 1 to 0.

We can also see that the width of the BAO peak, as indicated by the parameter
$\sigma^{\rm hh}$, changes from about 5.0 to 7.7 Mpc/$h$ from $z=1$ to $z=0$,
a 50% increase. This is expected from non-linear evolution (e.g. Seo et al.,
2008). In addition, the halo bias increases from $\sim 1.4$ to $\sim 2.8$ as
the redshift increases from 0 to 1. This is expected as the density peaks
corresponding to halos of the same low-mass cutoff become relatively rare at
high-$z$. Additionally, the bias parameter $b_{\rm halo}$ appears uncorrelated
with the $\alpha^{\rm hh}$, suggesting that the constraints on the BAO scale
is unaffected by the amplitude of the clustering, as expected. All these
indicate that our fitting pipeline works well for the 2PCF, and that there is
no obvious shifts for the BAO scale in the 2PCF in this redshift range.

Figure 5: The best-fit results of the auto-correlation functions (Halo $\times$ Halo), as well as the joint-fit results of the cross-correlation functions (Halo $\times$ DS) and the auto-correlation functions (DS $\times$ DS) at $z=0.0$, 0.5 and 1.0. Points with error bars represent data obtained from simulated samples. The cross sign (‘x’) denotes Halo $\times$ Halo correlation; solid and hollow circles are used to distinguish DS5 and DS1 for both Halo $\times$ DS and DS $\times$ DS correlation functions. Corresponding curves represent the best-fit models. Blue, yellow, and red lines correspond to $z=0.0$, 0.5 and 1.0, respectively.  Figure 6:  Panels A, B, C display the 2D posterior probability distributions of halo bias ($b_{\rm halo}$), the BAO shift parameter ($\alpha^{\rm hh}$), and the BAO width parameter ($\sigma^{\rm hh}$), obtained from MCMC fitting using halo auto-correlation functions. Large and small dashed contours represent the 68% and 95% confidence intervals, respectively. Panel D shows the 2D posterior probability distributions of the BAO shift parameters ($\alpha^{\rm qh}$ and $\alpha^{\rm qq}$) from joint constraints using DS $\times$ DS and Halo $\times$ DS correlation functions. Only the 68% confidence contours are shown, with ${\rm DS}_{1}$ and ${\rm DS}_{5}$ distinguished by dotted and solid lines, respectively. Best-fit $\alpha^{\rm qh}$-$\alpha^{\rm qq}$ values for ${\rm DS}_{1}$ and ${\rm DS}_{5}$ are marked with hollow-downward triangles and solid-upward triangles, respectively, while the best-fit $\alpha^{\rm hh}$ values are indicated with cross signs (“x”). Overall: Blue, yellow, and red lines correspond to $z=0.0$, 0.5 and 1.0, respectively. The marginalized 68% confidence intervals of the $\alpha^{\rm hh}$ parameter are shown as shaded bands in Panels B, C, and D for comparison. Additionally, both horizontal and vertical dashed black lines represent the unshifted alpha parameter ($\alpha=1$). Figure 7:  The marginalized constraints for the BAO shift parameter $\alpha$ and width parameter $\sigma$ at redshifts 0, 0.5, and 1. Black ‘x’ marks indicate the central values obtained from the halo auto-correlation function, while blue and red circles represent the central values obtained from the joint analysis of DS auto-correlation function and the Halo-DS cross-correlation function, respectively (hollow and solid circles distinguish DS1 and DS5). The 68% confidence intervals for each parameter are displayed as error bars in the plot. $z$ | $0.0$ | $0.5$ | $1.0$
---|---|---|---
$b_{\rm halo}$ | $1.425\pm{0.023}$ | $2.015\pm{0.031}$ | $2.832^{+0.042}_{-0.043}$
$\alpha^{\rm hh}$ | $1.0033^{+0.0025}_{-0.0024}$ | $1.0037\pm{0.0020}$ | $1.0039\pm{0.0019}$
$\sigma^{\rm hh}$ | $7.72\pm{0.29}$ | $6.18^{+0.27}_{-0.28}$ | $5.06^{+0.30}_{-0.31}$
$b_{\rm halo|DS_{1}}$ | $1.435\pm{0.023}$ | $2.041\pm{0.033}$ | $2.851^{+0.048}_{-0.049}$
$b_{{\rm DS}_{1}}$ | $-3.02\pm{0.05}$ | $-3.96^{+0.07}_{-0.06}$ | $-4.83^{+0.09}_{-0.08}$
$\alpha^{\rm qq}_{1}$ | $0.9928\pm{0.0031}$ | $0.9946^{+0.0031}_{-0.0030}$ | $0.9943\pm{0.0035}$
$\alpha^{\rm qh}_{1}$ | $0.9969\pm{0.0023}$ | $0.9990\pm{0.0021}$ | $0.9980\pm{0.0023}$
$\sigma^{\rm qq}_{1}$ | $13.78^{+0.25}_{-0.24}$ | $13.66\pm{0.24}$ | $13.76^{+0.28}_{-0.27}$
$\sigma^{\rm qh}_{1}$ | $11.12\pm{0.21}$ | $10.72^{+0.20}_{-0.19}$ | $10.34^{+0.22}_{-0.21}$
$b_{\rm halo|DS_{5}}$ | $1.423\pm{0.020}$ | $2.016^{+0.028}_{-0.029}$ | $2.817^{+0.043}_{-0.044}$
$b_{{\rm DS}_{5}}$ | $3.66\pm{0.05}$ | $4.69\pm{0.07}$ | $5.68\pm{0.09}$
$\alpha^{\rm qq}_{5}$ | $1.0015^{+0.0039}_{-0.0038}$ | $1.0008\pm{0.0034}$ | $0.9986\pm{0.0034}$
$\alpha^{\rm qh}_{5}$ | $1.0023\pm{0.0032}$ | $1.0017\pm{0.0028}$ | $1.0006^{+0.0028}_{-0.0027}$
$\sigma^{\rm qq}_{5}$ | $17.15^{+0.28}_{-0.26}$ | $16.06^{+0.25}_{-0.24}$ | $15.26^{+0.26}_{-0.25}$
$\sigma^{\rm qh}_{5}$ | $13.39^{+0.25}_{-0.24}$ | $12.22\pm{0.22}$ | $11.44^{+0.24}_{-0.23}$
Table 2: The central values of parameters along with their $1\sigma$ errors
from our BAO analyses fitted with $\xi^{\rm hh}$, $\xi^{\rm qq}$ and $\xi^{\rm
qh}$, after marginalizing over other parameters. Note that $\xi^{\rm qq}$ and
$\xi^{\rm qh}$ of the same DS bins are combined into a single data vector for
the fitting. The corresponding posterior distributions from which these values
are derived are shown in Figs. 8 & 9.

###  4.2 BAO constraints from DS-clustering

For density-split clustering, we combine ${\xi}_{i}^{\rm qq}$ and
${\xi}_{i}^{\rm qh}$ together to obtain joint constraints for the shift
parameter. To account for the correlations between ${\xi}_{i}^{\rm qq}$ and
${\xi}_{i}^{\rm qh}$ at each redshift, we compute their joint covariance
matrix and present the normalized correlation coefficients between
${\xi}_{i}^{\rm qq}$ and ${\xi}_{i}^{\rm qh}$ for both the DS1 and DS5 cases
in Fig. 4. It can be seen from the figure that there is a strong correlation
between ${\xi}_{i}^{\rm qq}$ and ${\xi}_{i}^{\rm qh}$. However, the
correlation coefficients are significantly smaller for DS1 than for DS5. The
cross-correlation functions are relatively more independent from their auto-
correlation functions around under-dense regions than the case around high-
density regions. This suggests that it is more valuable to combine
${\xi}_{i}^{\rm qq}$ and ${\xi}_{i}^{\rm qh}$ for DS1 – under-dense regions,
than for DS5 – over-dense regions.

We can see from Fig. 5 (B, C) that the model provides excellent fits to the
simulation measurements for both the auto-correlation functions (DS $\times$
DS) and the cross-correlation functions (Halo $\times$ DS) for ${\rm DS}_{1}$
and ${\rm DS}_{5}$ centers.

####  4.2.1 Biases

Both the theoretical model of the auto-correlation functions ${\xi}_{i}^{\rm
qq}(r;b_{{\rm DS}_{i}},\alpha^{\rm qq}_{i},\sigma^{\rm qq}_{i})$ and the
cross-correlation functions ${\xi}_{i}^{\rm qh}(r;b_{{\rm
halo|DS}_{i}},b_{{\rm DS}_{i}},\alpha^{\rm qh}_{i},\sigma^{\rm qh}_{i})$ have
the parameter $b_{{\rm DS}_{i}}$. Their combination is able to break the
degeneracy between $b_{\rm halo}$ and $b_{{\rm DS}_{i}}$, allowing independent
measurements for the halo bias in different environmental densities. As shown
in Table 2, the biases of ${\rm DS}_{1}$ differ significantly from those of
${\rm DS}_{5}$, with $b_{{\rm DS}_{1}}\sim-3.0$ and $b_{{\rm DS}_{5}}\sim 3.6$
at $z=0$. These are significantly larger than the biases of the halos,
suggesting that large-scale density fluctuations of these kinds are indeed
rarer than the halos in our sample, but they are also more likely to be
responsible for inducing non-linear evolution for the BAO. Furthermore, the
biases for both ${\rm DS}_{1}$ and ${\rm DS}_{5}$ increase with redshift,
similar to the evolution of the halo bias. Additionally, the halo biases in
$\rm DS_{1}$ and $\rm DS_{5}$ environments are consistent with each other at
the $1\sigma$ level for all the three redshifts. They are also consistent with
the global halo bias coming from the fitting with the 2PCF’s. These results
suggest that to the linear order, the halos of the same mass are statistically
similar in different environments, once the environmental factor is accounted
for by the DS biases. This agrees with the theoretical expectations (Bardeen
et al., 1986; Bond et al., 1991; Martino & Sheth, 2009; Sheth & Tormen, 1999)
and findings from simulations (e.g. Alonso et al., 2015; Goh et al., 2019).

####  4.2.2 BAO scales

Fig. 6(D) presents the 2D posterior distributions of the shift parameters
$\alpha^{\rm qq}_{i}$ and $\alpha^{\rm qh}_{i}$ for ${\rm DS}_{1}$ and ${\rm
DS}_{5}$, along with $\alpha^{\rm hh}$ at redshifts $z=0.0$, 0.5, and 1.0.
Additionally, to facilitate more intuitive comparisons at different redshifts,
we present the central values of the BAO shift parameters $\alpha$ with
1$\sigma$ error bars in Fig. 7. The complete posterior distributions for ${\rm
DS}_{1}$ and ${\rm DS}_{5}$ are provided in Appendix A. Based on the
information provided in Fig. 6, Fig. 7, and Table 2, the following insights
can be intuitively gleaned:

$\bullet$ Both the $\alpha^{\rm qq}$ and the $\alpha^{\rm qh}$ parameters for
${\rm DS}_{5}$ tend to be larger than those for ${\rm DS}_{1}$. This suggest
that there are indeed contractions of the BAO scale around high-density
regions compared to low-density regions.

$\bullet$ As the simulations evolves from redshift 1 to 0, the values of the
$\alpha^{\rm qq}$ and $\alpha^{\rm qh}$ parameters for ${\rm DS}_{5}$
increase, whereas ${\rm DS}_{1}$ exhibits a mild opposite trend. Although the
best-fit $\alpha^{\rm qq}$-$\alpha^{\rm qh}$ values for ${\rm DS}_{1}$ at
$z=0.5$ are slightly higher than those at $z=1$, but they are consistent
within the $1\sigma$ confidence level, which is probably due to the
insignificant evolution of the BAO shifts between redshift 1 and 0.5. The
evolution becomes more obvious when we split the density field into 10 DS bins
(see Fig. 10(c)). In the one-dimensional posterior distributions of the
parameters at redshift 0, the difference between $\alpha^{\rm qh}_{1}$ and
$\alpha^{\rm qh}_{5}$ nearly reaches $1\sigma$, while the difference between
$\alpha^{\rm qq}_{1}$ and $\alpha^{\rm qq}_{5}$ exceeds $1\sigma$. The
differences in their best-fit values are approximately 0.5% and 1%,
respectively. The fact that the $\alpha$ parameters around the low and high
density regions evolve in opposite directions confirms the
expansion/contraction of them in these two different environments due to non-
linear evolution.

$\bullet$ The constrained results of parameters associated with the width of
the BAO, $\sigma^{\rm qh}_{i}$ and $\sigma^{\rm qq}_{i}$, are larger than
those of the 2PCF version $\sigma^{\rm hh}$ (bottom panel of Fig. 7). This is
mainly due to the convolution of the top-hat window function when splitting
the densities. In addition, the width of the BAO for $\rm DS_{1}$ is generally
smaller than that for $\rm DS_{5}$ (dotted versus solid lines at the bottom
panel of Fig. 7), which is consistent with the results of (Achitouv & Blake,
2015). This difference is possibly due to the smaller dynamical range of
density for $\rm DS_{1}$ than $\rm DS_{5}$ (Fig. 1). They also exhibit
different evolution between redshift one to zero. Both $\sigma^{\rm qh}_{5}$
and $\sigma^{\rm qq}_{5}$ increases with decreasing redshift, similar to the
case of halo-halo clustering. However, both $\sigma^{\rm qh}_{1}$ and
$\sigma^{\rm qq}_{1}$ remains stable across this redshift range.

The increase for the width of the BAO in the 2PCF as the universe evolve is
expected from non-linear evolution. As halos reside on the peaks of the
density field, they occupy a wide range of densities in the density PDF. The
BAO’s around different local densities evolve at different pace, resulting in
a broadening of its width on average. In DS-clustering, $\rm DS_{5}$ takes the
positive tail of the non-Gaussian PDF of the density field (see Fig. 1), which
also have a relatively wide dynamical range in density. Thus it can be
expected that the width of the BAO around $\rm DS_{5}$ centers follows a
similar evolution trend as that of the 2PCF. For $\rm DS_{1}$, however, their
local densities have a narrow dynamical range (see Fig. 1). The density
perturbations are likely to evolve at a similar pace i.e., BAO’s around under-
dense regions will expand by a similar rate as they evolve. Their widths will
therefore remain stable around these under-dense regions.

In summary, from the prospective of DS-clustering, BAO’s around under- and
over-dense regions evolve at different rates. By analyzing them together, such
as the case of the 2PCF, the situations in low- and high-density regions will
be mixed up. The expansion from under-dense regions and the contraction from
high-density regions will act together to broaden the width of the BAO.

##  5 Discussion and conclusions

In this study, we have explored the possible deviation of the BAO scale in
high-density and low-density regions due to non-linear evolution. We do so by
analyzing the BAO with the density-split statistics using halos found in the
Quijote simulations. The density-split clustering starts from splitting the
halo number density PDF smoothed at a 20 Mpc/$h$-scale into high- and low-
density regions, called DS centers. We then measure the auto-correlations of
these DS centers and their cross-correlations with the halos. We also compare
the results with the standard analysis of the BAO using 2PCF. The main
conclusions drawn from our analyses are as the follows:

For the scales of the BAO, using the the auto-correlation functions of DS
centers and the cross-correlation functions between halos and DS centers, we
observed a contraction (stretching) of the BAO scale in high-density (low-
density) regions relative to its linear version, while the scale of the BAO in
the auto-correlation function of halos remains stable in the same redshift
range. At redshift zero, the difference in the BAO scale between high- and
low-density regions is approximately 1%.

Two-point statistics can be seen as the composition of the whole spectrum of
cross-correlations in DS-clustering. The sub-percent level shifts for the BAO
in 2PCF is the sum of all the effects around high- and low-density regions in
DS-clustering. Due to the opposite direction for the shifts of the BAO scale
around high- and low-densities, a non-linear transformation for the density
field which up-weight low-densities may compensate the contraction effect of
the BAO around high-densities, and therefore reduce the shifts of the BAO in
two-point statistics (McCullagh et al., 2013); a similar finding using the
power spectrum, for general features, was explored by Neyrinck & Yang (2013).

For the width of the BAO, with density-split clustering, we find it to be
larger than its 2PCF counterpart due to a convolution of a smoothing window
function. The width is broadened by $\sim 50\%$ between redshift 1 to 0 for
the 2PCF. Similar evolution is seen for DS-clustering for the high-density
regions. For the highest density bin in DS-clustering, DS5, it encloses the
positive tail of the density PDF, which is more extended than the lowest
density bin. The wide dynamical range of the local density, and the relatively
strong non-linearity around these high-densities may explain the relatively
more rapid evolution for the width of the BAO around them. The width of the
BAO does not appear to evolve for low-density regions, possibly due to having
a similar evolution trend for the BAO around a narrow range of low-density
environments i.e., the evolution of BAO’s around voids may be relatively
synchronized. This may also explains the fact that reconstructions by up-
weighting low-density environments is beneficial in sharpening the BAO feature
(Achitouv & Blake, 2015), which is also seen in the reconstruction with mock
galaxy catalogs (Zhao et al., 2020).

Our results agree qualitatively with the findings of (Neyrinck et al., 2018),
but the magnitudes of the shifts are smaller. This is mainly because we are
looking at a more realistic case where the density field is split into 5
quintiles. The local densities for $\rm DS_{1}$ and $\rm DS_{5}$ are not very
extreme. Their density contrasts $\delta$ differ from the mean by an order of
1. More significant shifts for the BAO scale is expected if one samples the
density contract of $\delta\sim 10$ (Neyrinck et al., 2018). However, at the
smoothing scale of 20 Mpc/$h$, density fluctuations of this amplitude are
rare, and the measurements for the DS-clustering become noisy. With a smaller
smoothing scale, one can increase the amplitude of density fluctuations while
keeping the noise level down. However, having smaller scale fluctuations means
weaker gravitational influence on the BAO, and hence smaller impacts on the
BAO scale.

Nevertheless, with the default smoothing scale of $20$Mpc/$h$, we have
repeated the analyses with 10 DS bins. Results are shown in Appendix B. In
this case, $\rm DS_{1}$ and $\rm DS_{10}$ are more extreme low and high
densities. We have indeed observe a more obvious trend for the expansion of
the BAO scale for $\rm DS_{1}$, but the evolution for the BAO scale for $\rm
DS_{10}$ seems diminishing. The difference between these low and high
densities becomes more extreme, reaching $\sim 2$% at $z=0$. Overall, the
results from using 10 DS-bins remains qualitatively similar to the case of 5
bins.

In conclusion, we have detected clear deviations for the BAO features between
the low- and high- density environments, both in terms of the shift of its
scale and in its broadening. These results may provide theoretical grounds for
better interpretations of the BAO beyond two-point statistics. Accurate
modeling for the environmental dependents of these features, e.g. through
density-split clustering, may allow us to better predict the non-linear
behavior of the BAO, and thus extract more cosmological information.

We caution that with a focus on understanding the physics, our analyses are
conducted in real space, but we anticipate the main conclusions to be similar
in redshift space. If we take the perspective that the scenario in redshift
space is a future universe in one of the three dimensions, the evolution
effects of the BAO we have seen may be even stronger in redshift space.
Indeed, this was seen for the shifts of the BAO in terms of the 2PCF (Seo et
al., 2008).

This work has been supported by the National Natural Science Foundation of
China (Nos. 11988101 and 12033008), the National Key Research and Development
Program of China (Nos. 2022YFA1602903 and 2023YFB3002501), the China Manned
Space Project (No. CMS-CSST-2021-A01), and the K. C. Wong Education
Foundation. YC acknowledges the support of the UK Royal Society through a
University Research Fellowship. WQ acknowledges the supported by the Strategic
Priority Research Program of Chinese Academy of Sciences, Grant No.XDB0500203.
For the purpose of open access, the author has applied a Creative Commons
Attribution (CC BY) licence to any Author Accepted Manuscript version arising
from this submission.

##  Appendix A Complete posterior probability distributions from the density-
split clustering

We display the 1D and 2D probability distributions of parameters of interest
constrained by the auto-correlation and cross-correlation functions of $\rm
DS_{1}$ (Fig. 8) and $\rm DS_{5}$ (Fig. 9) at $z=0$, $0.5$ and $1$ in this
section. We split the halo density field into 5 bins, the default case for the
paper. Results presented in Table 2 are derived from these figures.

Figure 8: The 1D and 2D probability distributions of each pair of parameters
constrained by the auto-correlation and cross-correlation functions of $\rm
DS_{1}$ at $z=$0, 0.5 and 1. The contours correspond to 68% and 95% confidence
levels. The density field is split into 5 density bins – the default choice of
the paper. Results presented in Table 2 are derived from these. Figure 9: The
1D and 2D probability distributions of each pair of parameters constrained by
the auto-correlation and cross-correlation functions of $\rm DS_{5}$ at $z=$0,
0.5 and 1. The contours correspond to 68% and 95% confidence levels. The
density field is split into 5 density bins– the default choice of the paper.
Results presented in Table 2 are derived from these.

##  Appendix B posterior probability distributions from the density-split
clustering with 10 bins

We repeat our analyses by splitting the halo density field into 10 density
bins, with results shown in Figures 10(a), 10(b), 10(c) and 10(d). These can
be compared with Figures 4, 5, 6 and 7 for the default case of having 5
density bins presented in the main part of the paper. The central values for
the parameters with errors are presented in Table 3, to be compared with Table
2.

(a) For the 10 density bin case, comparing the correlation coefficients
between the DS-DS and Halo-DS clustering at $z=0.0$, 0.5 and 1.0, where the
dotted lines represent DS1 and the solid lines represent DS10. The results are
similar to those with 5 density bins presented in Fig. 4, except that the
magnitudes of the correlation coefficients are all smaller, and there is no
significant evolution in the correlation coefficients for high-density regions
as the redshift changes from 0 to 1.

(b) For the case of 10 density bins, the best-fit results of the auto-
correlation functions (Halo $\times$ Halo), as well as the joint fit results
of the cross-correlation functions (Halo $\times$ DS) and the auto-correlation
functions (DS $\times$ DS) at $z=0.0$, 0.5 and 1.0. Similar to the 5 density
bins case (Fig. 5), the correlation function models we used can well fit the
measurements obtained from the simulation data.

(c) For the 10 density bin case, the 2D probability distributions of the BAO
shift parameters ($\alpha^{\rm qh}$ and $\alpha^{\rm qq}$) from joint
constraints using DS $\times$ DS and Halo $\times$ DS correlation functions,
with the details similar to Fig. 6. As the redshift evolves from 1 to 0, the
contraction of the BAO scale in high-density regions is insignificant, whereas
the expansion in low-density regions is pronounced. At $z=0$, the differences
between high- and low-density regions is larger than the 5-bin case shown in
Fig. 6.

(d) For the 10 density bin case, the marginalized constraints for the BAO
shift parameter $\alpha$ and width parameter $\sigma$ at redshifts 0, 0.5, and
1. The evolution of the $\alpha$ parameter with redshift is similar to the
case of 5 density bins shown in Fig. 7, except that the evolution in high-
density regions becomes less significant while the evolution in low-density
regions becomes more pronounced. The evolution of the $\sigma$ parameter is
consistent with the case of 5 density bins, indicating that the width of the
BAO peak in low-density regions is indeed stable in this redshift range.

Figure 10: Main results of density-split analysis with 10 density bins. $z$ | $0.0$ | $0.5$ | $1.0$
---|---|---|---
$b_{\rm halo}$ | $1.425\pm{0.023}$ | $2.015\pm{0.031}$ | $2.832\pm{0.043}$
$\alpha^{\rm hh}$ | $1.0034^{+0.0025}_{-0.0024}$ | $1.0036\pm{0.0020}$ | $1.0039\pm{0.0019}$
$\sigma^{\rm hh}$ | $7.72\pm{0.29}$ | $6.18^{+0.27}_{-0.28}$ | $5.06^{+0.29}_{-0.31}$
$b_{\rm halo|DS_{1}}$ | $1.426^{+0.026}_{-0.027}$ | $2.019\pm{0.036}$ | $2.839^{+0.051}_{-0.052}$
$b_{{\rm DS}_{1}}$ | $-3.56\pm{0.07}$ | $-4.68\pm{0.09}$ | $-5.64\pm{0.11}$
$\alpha^{\rm qq}_{1}$ | $0.9876\pm{0.0036}$ | $0.9897^{+0.0035}_{-0.0034}$ | $0.9952\pm{0.0041}$
$\alpha^{\rm qh}_{1}$ | $0.9952\pm{0.0023}$ | $0.9969\pm{0.0021}$ | $0.9987\pm{0.0023}$
$\sigma^{\rm qq}_{1}$ | $12.85^{+0.30}_{-0.29}$ | $12.90^{+0.30}_{-0.29}$ | $13.50^{+0.32}_{-0.31}$
$\sigma^{\rm qh}_{1}$ | $10.68^{+0.22}_{-0.21}$ | $10.26^{+0.21}_{-0.20}$ | $10.13\pm{0.22}$
$b_{\rm halo|DS_{10}}$ | $1.423\pm{0.021}$ | $2.001^{+0.029}_{-0.030}$ | $2.826^{+0.042}_{-0.043}$
$b_{{\rm DS}_{10}}$ | $4.83\pm{0.07}$ | $6.10^{+0.09}_{-0.10}$ | $7.45^{+0.12}_{-0.13}$
$\alpha^{\rm qq}_{10}$ | $1.0029\pm{0.0053}$ | $1.0022^{+0.0046}_{-0.0045}$ | $1.0020\pm{0.0045}$
$\alpha^{\rm qh}_{10}$ | $1.0039^{+0.0040}_{-0.0039}$ | $1.0031^{+0.0033}_{-0.0032}$ | $1.0037\pm{0.0032}$
$\sigma^{\rm qq}_{10}$ | $18.45^{+0.39}_{-0.36}$ | $16.90^{+0.33}_{-0.31}$ | $16.21^{+0.33}_{-0.31}$
$\sigma^{\rm qh}_{10}$ | $14.44^{+0.29}_{-0.28}$ | $12.90^{+0.26}_{-0.25}$ | $12.18^{+0.26}_{-0.25}$
Table 3: The central values of parameters along with their $1\sigma$ errors
from our BAO analyses using $\xi^{\rm hh}$, $\xi^{\rm qq}$ and $\xi^{\rm qh}$
for the case of $M=10$, after marginalizing over other parameters. Note that
$\xi^{\rm qq}$ and $\xi^{\rm qh}$ of the same DS bins are combined into a
single data vector for the fitting. The results in this table are to be
compared with the results for the case of 5 density bins shown in Table. 2.

## References

  * Achitouv & Blake (2015) Achitouv, I., & Blake, C. 2015, Phys. Rev. D, 92, 083523, doi: 10.1103/PhysRevD.92.083523
  * Alam et al. (2017) Alam, S., Ata, M., Bailey, S., et al. 2017, MNRAS, 470, 2617, doi: 10.1093/mnras/stx721
  * Alam et al. (2021) Alam, S., Aubert, M., Avila, S., et al. 2021, Phys. Rev. D, 103, 083533, doi: 10.1103/PhysRevD.103.083533
  * Alcock & Paczynski (1979) Alcock, C., & Paczynski, B. 1979, Nature, 281, 358, doi: 10.1038/281358a0
  * Alonso et al. (2015) Alonso, D., Eardley, E., & Peacock, J. A. 2015, MNRAS, 447, 2683, doi: 10.1093/mnras/stu2632
  * Anderson et al. (2012) Anderson, L., Aubourg, E., Bailey, S., et al. 2012, MNRAS, 427, 3435, doi: 10.1111/j.1365-2966.2012.22066.x
  * Babić et al. (2022) Babić, I., Schmidt, F., & Tucci, B. 2022, J. Cosmology Astropart. Phys, 2022, 007, doi: 10.1088/1475-7516/2022/08/007
  * Ballinger et al. (1996) Ballinger, W. E., Peacock, J. A., & Heavens, A. F. 1996, MNRAS, 282, 877, doi: 10.1093/mnras/282.3.877
  * Bardeen et al. (1986) Bardeen, J. M., Bond, J. R., Kaiser, N., & Szalay, A. S. 1986, ApJ, 304, 15, doi: 10.1086/164143
  * Bennett et al. (2013) Bennett, C. L., Larson, D., Weiland, J. L., et al. 2013, ApJS, 208, 20, doi: 10.1088/0067-0049/208/2/20
  * Beutler et al. (2011) Beutler, F., Blake, C., Colless, M., et al. 2011, MNRAS, 416, 3017, doi: 10.1111/j.1365-2966.2011.19250.x
  * Bond et al. (1991) Bond, J. R., Cole, S., Efstathiou, G., & Kaiser, N. 1991, ApJ, 379, 440, doi: 10.1086/170520
  * Chan et al. (2024) Chan, K. C., Lu, G., & Wang, X. 2024, MNRAS, 529, 1667, doi: 10.1093/mnras/stae638
  * Cole et al. (2005) Cole, S., Percival, W. J., Peacock, J. A., et al. 2005, MNRAS, 362, 505, doi: 10.1111/j.1365-2966.2005.09318.x
  * Crocce & Scoccimarro (2008) Crocce, M., & Scoccimarro, R. 2008, Phys. Rev. D, 77, 023533, doi: 10.1103/PhysRevD.77.023533
  * Davis et al. (1985) Davis, M., Efstathiou, G., Frenk, C. S., & White, S. D. M. 1985, ApJ, 292, 371, doi: 10.1086/163168
  * DESI Collaboration et al. (2024) DESI Collaboration, Adame, A. G., Aguilar, J., et al. 2024, arXiv e-prints, arXiv:2404.03002, doi: 10.48550/arXiv.2404.03002
  * Eisenstein & Hu (1998) Eisenstein, D. J., & Hu, W. 1998, ApJ, 496, 605, doi: 10.1086/305424
  * Eisenstein et al. (2007) Eisenstein, D. J., Seo, H.-J., Sirko, E., & Spergel, D. N. 2007, ApJ, 664, 675, doi: 10.1086/518712
  * Eisenstein et al. (2005) Eisenstein, D. J., Zehavi, I., Hogg, D. W., et al. 2005, ApJ, 633, 560, doi: 10.1086/466512
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306, doi: 10.1086/670067
  * Goh et al. (2019) Goh, T., Primack, J., Lee, C. T., et al. 2019, MNRAS, 483, 2101, doi: 10.1093/mnras/sty3153
  * Hand et al. (2018) Hand, N., Feng, Y., Beutler, F., et al. 2018, AJ, 156, 160, doi: 10.3847/1538-3881/aadae0
  * Hu & Sugiyama (1996) Hu, W., & Sugiyama, N. 1996, ApJ, 471, 542, doi: 10.1086/177989
  * Kirkby et al. (2013) Kirkby, D., Margala, D., Slosar, A., et al. 2013, J. Cosmology Astropart. Phys, 2013, 024, doi: 10.1088/1475-7516/2013/03/024
  * Lewis (2019) Lewis, A. 2019, arXiv e-prints, arXiv:1910.13970, doi: 10.48550/arXiv.1910.13970
  * Mao et al. (2021) Mao, T.-X., Wang, J., Li, B., et al. 2021, MNRAS, 501, 1499, doi: 10.1093/mnras/staa3741
  * Martino & Sheth (2009) Martino, M. C., & Sheth, R. K. 2009, MNRAS, 394, 2109, doi: 10.1111/j.1365-2966.2009.14467.x
  * McCullagh et al. (2013) McCullagh, N., Neyrinck, M. C., Szapudi, I., & Szalay, A. S. 2013, ApJ, 763, L14, doi: 10.1088/2041-8205/763/1/L14
  * Neyrinck et al. (2018) Neyrinck, M. C., Szapudi, I., McCullagh, N., et al. 2018, MNRAS, 478, 2495, doi: 10.1093/mnras/sty1074
  * Neyrinck & Yang (2013) Neyrinck, M. C., & Yang, L. F. 2013, MNRAS, 433, 1628, doi: 10.1093/mnras/stt949
  * Nikakhtar et al. (2023) Nikakhtar, F., Padmanabhan, N., Lévy, B., Sheth, R. K., & Mohayaee, R. 2023, Phys. Rev. D, 108, 083534, doi: 10.1103/PhysRevD.108.083534
  * Nikakhtar et al. (2022) Nikakhtar, F., Sheth, R. K., Lévy, B., & Mohayaee, R. 2022, Phys. Rev. Lett., 129, 251101, doi: 10.1103/PhysRevLett.129.251101
  * Noh et al. (2009) Noh, Y., White, M., & Padmanabhan, N. 2009, Phys. Rev. D, 80, 123501, doi: 10.1103/PhysRevD.80.123501
  * Orban & Weinberg (2011) Orban, C., & Weinberg, D. H. 2011, Phys. Rev. D, 84, 063501, doi: 10.1103/PhysRevD.84.063501
  * Padmanabhan et al. (2009) Padmanabhan, N., White, M., & Cohn, J. D. 2009, Phys. Rev. D, 79, 063523, doi: 10.1103/PhysRevD.79.063523
  * Paillas et al. (2021) Paillas, E., Cai, Y.-C., Padilla, N., & Sánchez, A. G. 2021, MNRAS, 505, 5731, doi: 10.1093/mnras/stab1654
  * Paillas et al. (2023) Paillas, E., Cuesta-Lazaro, C., Zarrouk, P., et al. 2023, MNRAS, 522, 606, doi: 10.1093/mnras/stad1017
  * Paillas et al. (2024) Paillas, E., Ding, Z., Chen, X., et al. 2024, arXiv e-prints, arXiv:2404.03005, doi: 10.48550/arXiv.2404.03005
  * Planck Collaboration et al. (2020) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, A&A, 641, A6, doi: 10.1051/0004-6361/201833910
  * Ross et al. (2015) Ross, A. J., Samushia, L., Howlett, C., et al. 2015, MNRAS, 449, 835, doi: 10.1093/mnras/stv154
  * Schmittfull et al. (2015) Schmittfull, M., Feng, Y., Beutler, F., Sherwin, B., & Chu, M. Y. 2015, Phys. Rev. D, 92, 123522, doi: 10.1103/PhysRevD.92.123522
  * Seo et al. (2022) Seo, H.-J., Ota, A., Schmittfull, M., Saito, S., & Beutler, F. 2022, MNRAS, 511, 1557, doi: 10.1093/mnras/stac082
  * Seo et al. (2008) Seo, H.-J., Siegel, E. R., Eisenstein, D. J., & White, M. 2008, ApJ, 686, 13, doi: 10.1086/589921
  * Seo et al. (2010) Seo, H.-J., Eckel, J., Eisenstein, D. J., et al. 2010, ApJ, 720, 1650, doi: 10.1088/0004-637X/720/2/1650
  * Sheth & Tormen (1999) Sheth, R. K., & Tormen, G. 1999, MNRAS, 308, 119, doi: 10.1046/j.1365-8711.1999.02692.x
  * Sinha & Garrison (2019) Sinha, M., & Garrison, L. 2019, in Software Challenges to Exascale Computing, ed. A. Majumdar & R. Arora (Singapore: Springer Singapore), 3–20.  https://doi.org/10.1007/978-981-13-7729-7_1
  * Sinha & Garrison (2020) Sinha, M., & Garrison, L. H. 2020, MNRAS, 491, 3022, doi: 10.1093/mnras/stz3157
  * Villaescusa-Navarro et al. (2020) Villaescusa-Navarro, F., Hahn, C., Massara, E., et al. 2020, ApJS, 250, 2, doi: 10.3847/1538-4365/ab9d82
  * von Hausegger et al. (2022) von Hausegger, S., Lévy, B., & Mohayaee, R. 2022, Phys. Rev. Lett., 128, 201302, doi: 10.1103/PhysRevLett.128.201302
  * Wang & Pen (2019) Wang, X., & Pen, U.-L. 2019, ApJ, 870, 116, doi: 10.3847/1538-4357/aaf231
  * Wang et al. (2017) Wang, X., Yu, H.-R., Zhu, H.-M., et al. 2017, ApJ, 841, L29, doi: 10.3847/2041-8213/aa738c
  * White (2010) White, M. 2010, arXiv e-prints, arXiv:1004.0250, doi: 10.48550/arXiv.1004.0250
  * White (2015) —. 2015, MNRAS, 450, 3822, doi: 10.1093/mnras/stv842
  * Yu et al. (2017) Yu, Y., Zhu, H.-M., & Pen, U.-L. 2017, ApJ, 847, 110, doi: 10.3847/1538-4357/aa89e7
  * Zarrouk et al. (2018) Zarrouk, P., Burtin, E., Gil-Marín, H., et al. 2018, MNRAS, 477, 1639, doi: 10.1093/mnras/sty506
  * Zhao et al. (2020) Zhao, C., Chuang, C.-H., Kitaura, F.-S., et al. 2020, MNRAS, 491, 4554, doi: 10.1093/mnras/stz3339


u)+(P_{L}-u\otimes u),$

implying that ${\rm
tr}(\Sigma_{1}^{-1}\Sigma_{0})-d=\frac{\lambda}{\lambda_{1}}-1$ and ${\rm
det}(\Sigma_{1}^{-1}\Sigma_{0})=\frac{\lambda}{\lambda_{1}}.$ Hence

$\displaystyle
K(N(0,\Sigma_{0})\|N(0,\Sigma_{1}))=\frac{\lambda}{\lambda_{1}}-1-\log\frac{\lambda}{\lambda_{1}}.$

Since

$\displaystyle\frac{\lambda}{\lambda_{1}}-1=\frac{\gamma_{1}a}{\gamma_{1}a+c\gamma_{1}\frac{a}{\sqrt{n}}}-1=\frac{1}{1+\frac{c}{\sqrt{n}}}-1=-\frac{\frac{c}{\sqrt{n}}}{1+\frac{c}{\sqrt{n}}},$

we get

$\displaystyle
K(N(0,\Sigma_{0})\|N(0,\Sigma_{1}))=-\frac{\frac{c}{\sqrt{n}}}{1+\frac{c}{\sqrt{n}}}-\log\biggl{(}1-\frac{\frac{c}{\sqrt{n}}}{1+\frac{c}{\sqrt{n}}}\biggr{)}\asymp\frac{c^{2}}{n},$

provided that $\frac{c}{\sqrt{n}}<1/4.$ This implies that

$\displaystyle K(N(0,\Sigma_{0})^{\otimes n}\|N(0,\Sigma_{1})^{\otimes
n})=nK(N(0,\Sigma_{0})\|N(0,\Sigma_{1}))\asymp c^{2}.$

Thus, if $c^{2}$ is small enough, then by [26], Theorem 2.1, Theorem 2.2(iii),
we can conclude that

$\displaystyle{\mathbb{E}}_{\Sigma_{0}}\phi(X_{1},\dots,X_{n})\bigvee{\mathbb{E}}_{\Sigma_{1}}(1-\phi(X_{1},\dots,X_{n}))\gtrsim
1.$

Combining this with (4.3) and (4.3) we get that

$\displaystyle\inf_{T}\sup_{\Sigma\in{\mathcal{S}}(a,r),\|\Sigma-\Sigma_{0}\|<\delta}{\mathbb{E}}_{\Sigma}(T(X_{1},\dots,X_{n})-f(\Sigma))^{2}\gtrsim\frac{\gamma_{1}^{2}a^{2}}{n},$

which implies the claim.

∎

The proof of the next proposition is much more involved.

######  Proposition 4.2.

Let $s>0.$ Suppose $\delta$ satisfies the condition

$\displaystyle
c_{1}\sqrt{\gamma_{1}\gamma_{2}}a\sqrt{\frac{r}{n}}\leq\delta\leq
c_{2}(\gamma_{1}-\gamma_{2})a\wedge 1$

for sufficiently large constant $c_{1}>0$ and sufficiently small constant
$c_{2}>0$ and let $U=B(\Sigma_{0},\delta).$ Then

$\displaystyle\sup_{\|f\|_{C^{s}(U)}\leq
1}\sup_{\Sigma\in{\mathcal{S}}(a,r),\|\Sigma-\Sigma_{0}\|<\delta}\Bigl{\|}T(X_{1},\dots,X_{n})-f(\Sigma)\Bigr{\|}_{L_{2}({\mathbb{P}}_{\Sigma})}\gtrsim\gamma_{1}^{s/2}\gamma_{2}^{s/2}(a\wedge
1)^{s}\Bigl{(}\sqrt{\frac{r}{n}}\Bigr{)}^{s}.$

proof. It will be enough to prove the bound for a sufficiently large $r,$ say,
$r\geq r_{0},$ where $r_{0}$ is a numerical constant. Indeed, then the bound
also holds for $r\leq r_{0}$ with a constant depending on $r_{0}.$

Recall that $d=[r].$ In what follows, we will use certain “well separated”
finite sets of unit vectors $\theta_{\omega}\in L$ parametrized by vertices
$\omega$ of binary cube $\\{-1,1\\}^{d}.$ Recall that $u\in L$ is a unit
vector and $\Sigma_{0}=\Sigma_{u,\lambda,\mu}.$ For $\omega\in\\{-1,1\\}^{d},$
define

$\displaystyle
t_{\omega}:=\varepsilon\frac{\omega}{\sqrt{d}}+\sqrt{1-\varepsilon^{2}}u$

(recall that vectors are replaced by their coordinate representations). Then

$\displaystyle\|t_{\omega}\|^{2}=1+\frac{2\varepsilon\sqrt{1-\varepsilon^{2}}}{\sqrt{d}}\langle
u,\omega\rangle.$ (4.32)

Let

$\displaystyle
h(\omega,\omega^{\prime}):=\sum_{j=1}^{d}I(\omega_{j}\neq\omega_{j}^{\prime}),\omega,\omega^{\prime}\in\\{-1,1\\}^{d}$

be the Hamming distance in the binary cube $\\{-1,1\\}^{d}.$ We will use the
following simple lemma (which is a modification of Varshamov-Gilbert bound).

######  Lemma 4.3.

There exists a subset $B\subset\\{-1,1\\}^{d}$ such that

  1. (i)

${\rm card}(B)\geq\frac{e^{d/8}}{2};$

  2. (ii)

$\langle u,\omega\rangle<2,\omega\in B;$

  3. (iii)

$h(\omega,\omega^{\prime})\geq d/4,\ \omega,\omega^{\prime}\in
B,\omega\neq\omega^{\prime}.$

proof. Assuming that ${\mathbb{P}}$ is a uniform probability distribution on
the binary cube $\\{-1,1\\}^{d},$ $\langle u,\omega\rangle$ is a subgaussian
r.v. with parameter $1,$ so, we have

$\displaystyle{\mathbb{P}}\\{\langle u,\omega\rangle\geq t\\}\leq
e^{-t^{2}/2}.$

Therefore,

$\displaystyle 2^{-d}{\rm card}\\{\omega\in\\{-1,1\\}^{d}:|\langle
u,\omega\rangle|\geq t\\}\leq 2e^{-t^{2}/2}.$

For $t=2,$ we have $2e^{-t^{2}/2}<1/2.$ This implies that

$\displaystyle{\rm card}\\{\omega\in\\{-1,1\\}^{d}:|\langle
u,\omega\rangle|<2\\}>2^{d-1}.$

Note that the cardinality of any ball of radius $m$ with respect to the
Hamming distance in the $d$-dimensional binary cube is equal to

$\displaystyle{d\choose{\leq m}}=\sum_{j=0}^{m}{d\choose j}.$

Let $B\subset\\{\omega\in\\{-1,1\\}^{d}:|\langle u,\omega\rangle|<2\\}$ be an
$m$-separated subset w.r.t. the Hamming distance of maximal cardinality. Then
$\\{\omega\in\\{-1,1\\}^{d}:|\langle u,\omega\rangle|<2\\}$ is covered by the
union of balls of radius $m$ with centers in $B,$ implying that ${\rm
card}(B){d\choose{\leq m}}\geq 2^{d-1},$ or

$\displaystyle{\rm card}(B)\frac{{d\choose{\leq m}}}{2^{d}}\geq 1/2.$

If $\nu\sim B(d,1/2)$ denotes the binomial r.v. with parameters $m$ and $1/2,$
then the last inequality can be interpreted as

$\displaystyle{\rm card}(B){\mathbb{P}}\\{\nu\leq m\\}\geq 1/2.$

We will choose $m=d/4.$ Then, by Hoeffding’s inequality,
${\mathbb{P}}\\{\nu\leq d/4\\}\leq e^{-d/8}.$ Thus, for the corresponding set
$B,$ we have ${\rm card}(B)\geq\frac{e^{d/8}}{2}.$ ∎

Denote

$\displaystyle\theta_{\omega}:=\frac{t_{\omega}}{\|t_{\omega}\|},\ \omega\in
B.$

######  Lemma 4.4.

Assume that $\frac{\varepsilon}{\sqrt{d}}\leq 1/8.$ Then, for all
$\omega,\omega^{\prime}\in B,$

$\displaystyle\frac{\varepsilon}{\sqrt{d}}(\sqrt{h(\omega,\omega^{\prime})}-4)\leq\|\theta_{\omega}-\theta_{\omega^{\prime}}\|\leq\frac{8\varepsilon}{\sqrt{d}}\sqrt{h(\omega,\omega^{\prime})}.$

Moreover, if $d$ is sufficiently large (namely, $d\geq 4\times 8^{2}$), then,
for all $\omega,\omega^{\prime}\in B,\omega\neq\omega^{\prime},$

$\displaystyle\frac{\varepsilon}{2\sqrt{d}}\sqrt{h(\omega,\omega^{\prime})}\leq\|\theta_{\omega}-\theta_{\omega^{\prime}}\|\leq\frac{8\varepsilon}{\sqrt{d}}\sqrt{h(\omega,\omega^{\prime})},$

which implies that

$\displaystyle\frac{\varepsilon}{4}\leq\|\theta_{\omega}-\theta_{\omega^{\prime}}\|\leq
8\varepsilon,\ \omega,\omega^{\prime}\in B,\omega\neq\omega^{\prime}.$

proof. Note that, in view of (4.32) for all $\omega\in B,$ we have

$\displaystyle 1-\frac{4\varepsilon}{\sqrt{d}}\leq\|t_{\omega}\|^{2}\leq
1+\frac{4\varepsilon}{\sqrt{d}},$

which also implies that

$\displaystyle 1-\frac{4\varepsilon}{\sqrt{d}}\leq\|t_{\omega}\|\leq
1+\frac{4\varepsilon}{\sqrt{d}},$ (4.33)

and, under the assumption that $\varepsilon/\sqrt{d}\leq 1/8,$ we have
$1/2\leq\|t_{\omega}\|\leq 3/2.$

Since
$t_{\omega}-t_{\omega^{\prime}}=\varepsilon\frac{\omega-\omega^{\prime}}{\sqrt{d}},$
we have

$\displaystyle\|t_{\omega}-t_{\omega^{\prime}}\|=\frac{2\varepsilon}{\sqrt{d}}\sqrt{h(\omega,\omega^{\prime})}.$
(4.34)

Note also that

$\displaystyle\|\theta_{\omega}-\theta_{\omega^{\prime}}\|\leq\Bigl{\|}\frac{t_{\omega}}{\|t_{\omega}\|}-\frac{t_{\omega^{\prime}}}{\|t_{\omega}\|}\Bigr{\|}+\Bigl{\|}\frac{t_{\omega^{\prime}}}{\|t_{\omega}\|}-\frac{t_{\omega^{\prime}}}{\|t_{\omega^{\prime}}\|}\Bigr{\|}$
$\displaystyle\leq\frac{\|t_{\omega}-t_{\omega^{\prime}}\|}{\|t_{\omega}\|}+\|t_{\omega^{\prime}}\|\Bigl{|}\frac{1}{\|t_{\omega}\|}-\frac{1}{\|t_{\omega^{\prime}}\|}\Bigr{|}\leq\frac{2\|t_{\omega}-t_{\omega^{\prime}}\|}{\|t_{\omega}\|}\leq
4\|t_{\omega}-t_{\omega^{\prime}}\|.$

On the other hand, it follows from (4.33) that
$|\|t_{\omega}\|-\|t_{\omega^{\prime}}\||\leq\frac{8\varepsilon}{\sqrt{d}},$
and we have

$\displaystyle\|t_{\omega}-t_{\omega^{\prime}}\|\leq\|t_{\omega}\|\|\theta_{\omega}-\theta_{\omega^{\prime}}\|+\|\theta_{\omega^{\prime}}\||\|t_{\omega}\|-\|t_{\omega^{\prime}}\||\leq
2\|\theta_{\omega}-\theta_{\omega^{\prime}}\|+\frac{8\varepsilon}{\sqrt{d}}.$

Thus, under the assumption $\frac{\varepsilon}{\sqrt{d}}\leq 1/8,$ we get that
for all $\omega,\omega^{\prime}\in B,$

$\displaystyle\frac{1}{2}\|t_{\omega}-t_{\omega^{\prime}}\|-\frac{4\varepsilon}{\sqrt{d}}\leq\|\theta_{\omega}-\theta_{\omega^{\prime}}\|\leq
4\|t_{\omega}-t_{\omega^{\prime}}\|.$

Combining this with (4.34) yields the first inequality of the lemma. For the
second inequality, note that, if $d\geq 4\times 8^{2},$ then, in view of bound
(iii) of Lemma 4.3, $\sqrt{h(\omega,\omega^{\prime})}\geq 8.$

∎

We need the following simple lemma.

######  Lemma 4.5.

Let $u_{1},u_{2}\in L$ with $\|u_{1}\|=1,\|u_{2}\|=1$ and let $\lambda>\mu>0.$
The following formula holds for the Kullback-Leibler distance between
$N(0,\Sigma_{u_{1},\lambda,\mu})$ and $N(0,\Sigma_{u_{2},\lambda,\mu}):$

$\displaystyle
K(N(0,\Sigma_{u_{1},\lambda,\mu})\|N(0,\Sigma_{u_{2},\lambda,\mu}))=\frac{(\lambda-\mu)(\mu^{-1}-\lambda^{-1})}{2}\|u_{1}\otimes
u_{1}-u_{2}\otimes u_{2}\|_{2}^{2}.$

proof. Indeed, taking into account that $\Sigma_{u_{1},\lambda,\mu}$ and
$\Sigma_{u_{2},\lambda,\mu}$ have the same spectrum and using well known
formulas for the KL-distance between two multivariate normal distributions, we
have

$\displaystyle
K(N(0,\Sigma_{u_{1},\lambda,\mu})\|N(0,\Sigma_{u_{2},\lambda,\mu}))={\rm
tr}(\Sigma_{u_{1},\lambda,\mu}^{-1}(\Sigma_{u_{2},\lambda,\mu}-\Sigma_{u_{1},\lambda,\mu})),$

where $\Sigma_{u_{1},\lambda,\mu}^{-1}$ is the inverse of operator
$\Sigma_{u_{1},\lambda,\mu}:L\mapsto L.$ Let $P_{1}=u_{1}\otimes u_{1}$
$P_{2}=u_{2}\otimes u_{2},$ $P=P_{L}.$ Note that
$\Sigma_{u_{2},\lambda,\mu}-\Sigma_{u_{1},\lambda,\mu}=(\lambda-\mu)(P_{2}-P_{1})$
and

$\displaystyle\Sigma_{u_{1},\lambda,\mu}^{-1}=\frac{1}{\lambda}P_{1}+\frac{1}{\mu}(P-P_{1})=\Bigl{(}\frac{1}{\lambda}-\frac{1}{\mu}\Bigr{)}P_{1}+\frac{1}{\mu}P.$

Therefore,

$\displaystyle{\rm
tr}(\Sigma_{u_{1},\lambda,\mu}^{-1}(\Sigma_{u_{2},\lambda,\mu}-\Sigma_{u_{1},\lambda,\mu}))=(\lambda-\mu){\rm
tr}\Bigl{(}\Bigl{(}\Bigl{(}\frac{1}{\lambda}-\frac{1}{\mu}\Bigr{)}P_{1}+\frac{1}{\mu}P\Bigr{)}(P_{2}-P_{1})\Bigr{)}$
$\displaystyle=(\lambda-\mu)(\lambda^{-1}-\mu^{-1}){\rm
tr}(P_{1}P_{2}-P_{1}^{2})+(\lambda-\mu)\mu^{-1}{\rm tr}(P(P_{2}-P_{1}))$
$\displaystyle=(\lambda-\mu)(\mu^{-1}-\lambda^{-1})(1-\langle
P_{1},P_{2}\rangle)=\frac{(\lambda-\mu)(\mu^{-1}-\lambda^{-1})}{2}\|P_{1}-P_{2}\|_{2}^{2},$

implying the claim.

∎

Let $X_{1},\dots,X_{n}$ be i.i.d. $N(0,\Sigma_{\theta,\lambda,\mu}),$ where
$\lambda>\mu>0$ are known numbers and $\theta\in L,\|\theta\|=1$ is an unknown
parameter to be estimated based on observations $X_{1},\dots,X_{n}.$ Further
assume that $\theta\in\Theta_{\varepsilon}:=\\{\theta_{\omega}:\omega\in
B\\}.$ The following fact and its proof are straightforward modifications of
well known minimax bounds on the error rate for principal eigenvector in
spiked covariance model.

######  Lemma 4.6.

Suppose that, for a sufficiently small constant $c>0,$

$\displaystyle\varepsilon^{2}\leq\frac{c}{(\lambda-\mu)(\mu^{-1}-\lambda^{-1})}\frac{d}{n}.$
(4.35)

Then

$\displaystyle\inf_{\hat{\theta}}\max_{\theta\in\Theta_{\varepsilon}}{\mathbb{E}}_{\theta}\|\hat{\theta}-\theta\|^{2}\gtrsim\varepsilon^{2}.$

proof. First note that, by Lemma 4.5 for all
$\theta,\theta^{\prime}\in\Theta_{\varepsilon},$

$\displaystyle K\Bigl{(}N(0,\Sigma_{\theta,\lambda,\mu})^{\otimes
n}\|N(0,\Sigma_{\theta^{\prime},\lambda,\mu})^{\otimes n}\Bigr{)}$
$\displaystyle=nK\Bigl{(}N(0,\Sigma_{\theta,\lambda,\mu})\|N(0,\Sigma_{\theta^{\prime},\lambda,\mu})\Bigr{)}$
$\displaystyle=\frac{(\lambda-\mu)(\mu^{-1}-\lambda^{-1})}{2}n\|\theta\otimes\theta-\theta^{\prime}\otimes\theta^{\prime}\|_{2}^{2}$
$\displaystyle=(\lambda-\mu)(\mu^{-1}-\lambda^{-1})n(1-\langle\theta,\theta^{\prime}\rangle^{2}).$

On the other hand, by the last bound of Lemma 4.4, for all
$\theta,\theta^{\prime}\in\Theta_{\varepsilon},$

$\displaystyle\langle\theta,\theta^{\prime}\rangle=\frac{2-\|\theta-\theta^{\prime}\|^{2}}{2}\geq
1-\frac{8^{2}\varepsilon^{2}}{2},$

implying that

$\displaystyle 1-\langle\theta,\theta^{\prime}\rangle^{2}\leq
1-\Bigl{(}1-\frac{8^{2}\varepsilon^{2}}{2}\Bigr{)}^{2}=\frac{8^{2}\varepsilon^{2}}{2}\Bigl{(}2-\frac{8^{2}\varepsilon^{2}}{2}\Bigr{)}\leq
8^{2}\varepsilon^{2}.$

Therefore, for all $\theta,\theta^{\prime}\in\Theta_{\varepsilon},$

$\displaystyle K\Bigl{(}N(0,\Sigma_{\theta,\lambda,\mu})^{\otimes
n}\|N(0,\Sigma_{\theta^{\prime},\lambda,\mu})^{\otimes n}\Bigr{)}\leq
8^{2}(\lambda-\mu)(\mu^{-1}-\lambda^{-1})n\varepsilon^{2}.$

Note also that, by Lemma 4.3 (i) under the assumptions $d\geq 4\times 8^{2}$
and (4.35) with a sufficiently small constant $c,$ we have

$\displaystyle\frac{1}{10}\log{\rm card}(\Theta_{\varepsilon})$
$\displaystyle=\frac{1}{10}\log{\rm card}(B)\geq\frac{1}{10}(d/8-\log 2)\geq
8^{2}(\lambda-\mu)(\mu^{-1}-\lambda^{-1})n\varepsilon^{2}$
$\displaystyle\geq\max_{\theta,\theta^{\prime}\in\Theta_{\varepsilon}}K\Bigl{(}N(0,\Sigma_{\theta,\lambda,\mu})^{\otimes
n}\|N(0,\Sigma_{\theta^{\prime},\lambda,\mu})^{\otimes n}\Bigr{)}.$

The claim now follows by an application of a standard Fano type inequality
(see, e.g., [26], Theorem 2.5).

∎

######  Lemma 4.7.

If $d\geq 8^{2}$ and $\frac{\varepsilon}{\sqrt{d}}\leq 1/8,$ then, for all
$\omega\in B,$

$\displaystyle\|\Sigma_{\theta_{\omega},\lambda,\mu}-\Sigma_{u,\lambda,\mu}\|\leq
5\sqrt{2}|\lambda-\mu|\varepsilon.$

proof. Note that
$\|t_{\omega}-u\|\leq\varepsilon+|1-\sqrt{1-\varepsilon^{2}}|\leq
2\varepsilon$ and, using bound (4.33), it is easy to get that, under the
assumptions $d\geq 8^{2}$ and $\frac{\varepsilon}{\sqrt{d}}\leq 1/8,$ for all
$\omega\in B,$

$\displaystyle\|\theta_{\omega}-u\|\leq\frac{\|t_{\omega}-u\|}{\|t_{\omega}\|}+\frac{|\|t_{\omega}\|-1|}{\|t_{\omega}\|}\leq
4\varepsilon+\frac{8\varepsilon}{\sqrt{d}}\leq 5\varepsilon.$

In addition, for all $\omega\in B,$

$\displaystyle\|\Sigma_{\theta_{\omega},\lambda,\mu}-\Sigma_{u,\lambda,\mu}\|$
$\displaystyle=|\lambda-\mu|\|\theta_{\omega}\otimes\theta_{\omega}-u\otimes
u\|\leq|\lambda-\mu|\|\theta_{\omega}\otimes\theta_{\omega}-u\otimes u\|_{2}$
$\displaystyle=|\lambda-\mu|\sqrt{2-2\langle\theta_{\omega},u\rangle^{2}}\leq|\lambda-\mu|\sqrt{4(1-\langle\theta_{\omega},u\rangle)}=\sqrt{2}|\lambda-\mu|\|\theta_{\omega}-u\|$
$\displaystyle\leq 5\sqrt{2}|\lambda-\mu|\varepsilon.$

∎

Let $\Sigma_{0}:=\Sigma_{u,\lambda,\mu}$ and let $\delta<(\lambda-\mu)/2.$
Then, for any covariance operator $\Sigma$ such that
$\|\Sigma-\Sigma_{0}\|<\delta,$ the top eigenvalue of
$\lambda(\Sigma)=\|\Sigma\|$ of $\Sigma$ is simple and
$|\lambda(\Sigma)-\lambda|<\delta.$ Moreover, if $\theta(\Sigma)$ is a unit
eigenvector of $\Sigma$ corresponding to $\lambda(\Sigma),$ then, by standard
perturbation bounds (see, e.g., [10]), bound (2.7) of Lemma 1)

$\displaystyle\|\theta(\Sigma)\otimes\theta(\Sigma)-u\otimes
u\|\leq\frac{4\delta}{\lambda-\mu}.$ (4.36)

Since eigenvector $\theta(\Sigma)$ is defined up to its sign, we will assume
that $\langle\theta(\Sigma),u\rangle\geq 0.$ We will need the following lemma.

######  Lemma 4.8.

Let $\bar{\delta}:=(\lambda-\mu)/8\wedge 1$ and $U:=B(\Sigma_{0},2\delta)$
with $\delta<\bar{\delta}/2.$ Then, for all $s>0,$ the mapping
$\Sigma\mapsto\theta(\Sigma)$ belongs to the Hölder space $C^{s}(U)$ with
$\|\theta(\cdot)\|_{C^{s}(U)}\lesssim_{s}\bar{\delta}^{-s}.$

proof. Let $g:{\mathbb{R}}\mapsto{\mathbb{R}}$ be a function supported in
$(\lambda-\bar{\delta},\lambda+\bar{\delta})$ such that $0\leq g(x)\leq
1,x\in{\mathbb{R}},$
$g(x)=1,x\in(\lambda-\bar{\delta}/2,\lambda+\bar{\delta}/2)$ and
$\|g\|_{B_{\infty,1}^{s}}\lesssim\bar{\delta}^{-s}$ for all $s>0.$ For
instance, one can take a $C^{\infty}$ function $\psi:{\mathbb{R}}\mapsto[0,1]$
with ${\rm supp}(\psi)\subset[-1,1]$ and such that $\psi(x)=1,x\in[-1/2,1/2]$
and set $g(x):=\psi(x/\bar{\delta}),x\in{\mathbb{R}}$ (for the proof of the
bound on the Besov norms, see [27], Section 3.4.1). Note that
$g(\Sigma)=\theta(\Sigma)\otimes\theta(\Sigma)$ and operator function
$A\mapsto g(A)$ 111in the space of bounded self-adjoint operators in
${\mathbb{H}}$ equipped with the operator norm belongs to the Hölder space
$C^{s}$ with $\|g\|_{C^{s}}\lesssim_{s}\bar{\delta}^{-s}$ for all $s>0$ (see
[7], Section 2, Corollary 2). In addition, bound (4.36) implies that
$g(U)\subset B(P,1/2),$ where $P=u\otimes u.$

Consider now the following mapping $A\mapsto F(A)$ from the ball $B(P,1/2)$
into ${\mathbb{H}}:$

$\displaystyle F(A):=\frac{Au}{\sqrt{\langle Au,u\rangle}}.$

Note that, for all $A\in B(P,1/2),$ we have $\langle Au,u\rangle>1/2,$ so, the
above mapping is well defined. Moreover, since $A\mapsto Au$ and
$A\mapsto\langle Au,u\rangle$ are both linear mappings with operator norms
bounded by $1$ and $\langle Au,u\rangle>1/2$ for $A\in B(P,1/2),$ it is easy
to check that the mapping $A\mapsto F(A)$ belongs to the Hölder space
$C^{s}(B(P,1/2))$ and $\|F\|_{C^{s}(B(P,1/2))}\lesssim_{s}1.$ Finally, note
that for any unit vector $v\in{\mathbb{H}}$ with $v\otimes v\in B(P,1/2)$ and
$\langle v,u\rangle\geq 0,$ we have $F(v\otimes v)=v.$ Therefore,
$\theta(\Sigma)=F(g(\Sigma)),\Sigma\in B(\Sigma_{0},2\delta)$ and, using Faà
di Bruno type calculus (see, e.g., [9], Section 3.1), it is easy to bound the
Hölder norms of this superposition, yielding the bound
$\|\theta(\cdot)\|_{C^{s}(U)}=\|F\circ
g\|_{C^{s}(U)}\lesssim_{s}\bar{\delta}^{-s}.$

∎

Our next goal is to construct a set of functionals $f_{k}\in
C^{s}(U),k=1,\dots,d,$ defined as $f_{k}(\Sigma)=h_{k}(\theta(\Sigma)),$ where
$h_{k}\in C^{s}({\mathbb{H}}).$ Each functional $h_{k}$ will be a sum of
“bumps” with disjoint supports around well separated vectors
$\theta_{\omega},\omega\in B.$ Namely, let $\varphi:{\mathbb{R}}\mapsto[0,1]$
be a $C^{\infty}$ function with support in $[-1,1]$ and with $\varphi(0)>0.$
Suppose that $\|\varphi\|_{C^{s}({\mathbb{R}})}\lesssim 1.$ Let
$\phi(u):=\varphi(\|u\|^{2}),u\in{\mathbb{H}}.$ Then
$\|\phi\|_{C^{s}({\mathbb{H}})}\lesssim 1.$ Define

$\displaystyle h_{k}(\theta):=\sum_{\omega\in
B}\omega_{k}\varepsilon^{s}\phi\Bigl{(}\frac{\theta-\theta_{\omega}}{c\varepsilon}\Bigr{)},\theta\in{\mathbb{H}},k=1,\dots,d,$

where $c\in(0,1/16)$ is a small enough constant. The following facts are
straightforward:

  1. (i)

each function
$\theta\mapsto\phi\Bigl{(}\frac{\theta-\theta_{\omega}}{c\varepsilon}\Bigr{)}$
is supported in a ball of radius $c\varepsilon$ centered at $\theta_{\omega}.$

  2. (ii)

Functions
$\phi\Bigl{(}\frac{\theta-\theta_{\omega}}{c\varepsilon}\Bigr{)},\omega\in B$
have disjoint supports. Moreover, the supports of any two of those functions
are separated by distance $\geq\varepsilon/8.$

  3. (iii)

For all $\omega\in B,$

$\displaystyle\Bigl{\|}\varepsilon^{s}\phi\Bigl{(}\frac{\cdot-\theta_{\omega}}{c\varepsilon}\Bigr{)}\Bigr{\|}_{C^{s}({\mathbb{H}})}\lesssim
1.$

The next proposition easily follows from the above properties.

######  Proposition 4.3.

  1. (i)

For all $k=1,\dots,d,$

$\displaystyle\|h_{k}\|_{C^{s}({\mathbb{H}})}\lesssim 1.$

  2. (ii)

For all $k=1,\dots,d$ and all $\omega\in B,$

$\displaystyle h_{k}(\theta_{\omega})=\omega_{k}\varepsilon^{s}\varphi(0).$

The second claim of Proposition 4.3 means that the values of functionals
$h_{k},k=1,\dots,d$ at the vectors $\theta_{\omega},\omega\in B$ could be used
to recover the vectors $\omega\in B$ and, hence, the vectors $\theta_{\omega}$
themselves. This property plays a crucial role in our argument (its idea goes
back to [22]).

As in the statement of the theorem, we fix $\lambda=\gamma_{1}a$ and
$\mu=\gamma_{2}a$ for some numerical constants $\gamma_{1}>\gamma_{2}>0$ and
denote $\Sigma_{\theta}:=\Sigma_{\theta,\lambda,\mu}.$ Recall that
$\Sigma_{0}=\Sigma_{u}$ and also recall the statement and the notations of
Lemma 4.8. Define functionals

$\displaystyle
f_{k}(\Sigma):=\bar{\delta}^{s}h_{k}(\theta(\Sigma)),k=1,\dots,d,\Sigma\in
U=B(\Sigma_{0},2\delta).$

Then $\|f_{k}\|_{C^{s}(U)}\lesssim 1$ and, moreover, by using
$c^{\prime}\varphi$ with a small enough $c^{\prime}>0$ instead of $\varphi,$
we can achieve $\|f_{k}\|_{C^{s}(U)}\leq 1.$ Note also that
$f_{k}(\Sigma_{\theta_{\omega}})=\bar{\delta}^{s}h_{k}(\theta_{\omega})=\varphi(0)\bar{\delta}^{s}\varepsilon^{s}\omega_{k},k=1,\dots,d.$
The last relationship clearly means that the values
$f_{k}(\Sigma_{\theta_{\omega}}),k=1,\dots,d$ could be used to recover vectors
$\omega\in B$ and $\theta_{\omega}.$ Moreover, define the following distance

$\displaystyle\tau(\omega,\omega^{\prime}):=\biggl{(}\frac{1}{d}\sum_{k=1}^{d}(f_{k}(\Sigma_{\omega})-f_{k}(\Sigma_{\omega^{\prime}}))^{2}\biggr{)}^{1/2},\
\omega,\omega^{\prime}\in\\{-1,1\\}^{d}.$

Then

$\displaystyle\tau(\omega,\omega^{\prime})=\frac{\varphi(0)\bar{\delta}^{s}\varepsilon^{s}}{\sqrt{d}}\sqrt{h(\omega,\omega^{\prime})},\
\omega,\omega^{\prime}\in\\{-1,1\\}^{d}.$

By the second bound of Lemma 4.4, we easily get the following statement:

######  Lemma 4.9.

Suppose $d\geq 4\times 8^{2}$ and $\frac{\varepsilon}{\sqrt{d}}\leq 1/8.$
Then, for all $\omega,\omega^{\prime}\in B,$

$\displaystyle\frac{1}{8}\varphi(0)\bar{\delta}^{s}\varepsilon^{s-1}\|\theta_{\omega}-\theta_{\omega^{\prime}}\|\leq\tau(\omega,\omega^{\prime})\leq
2\varphi(0)\bar{\delta}^{s}\varepsilon^{s-1}\|\theta_{\omega}-\theta_{\omega^{\prime}}\|.$

Consider now arbitrary estimators $T_{k}(X_{1},\dots,X_{n})$ of functionals
$f_{k}(\Sigma).$ For $k=1,\dots,d,$ define $\tilde{\omega}_{k}:={\rm
sign}(T_{k}(X_{1},\dots,X_{n})),$
$\tilde{\omega}:=(\tilde{\omega}_{1},\dots,\tilde{\omega}_{d})$ and

$\displaystyle\tilde{T}_{k}(X_{1},\dots,X_{n}):=\varphi(0)\bar{\delta}^{s}\varepsilon^{s}\tilde{\omega}_{k}=f_{k}(\Sigma_{\theta_{\tilde{\omega}}}).$

Then, for all $\omega\in B$ and $k=1,\dots,d,$

$\displaystyle{\mathbb{E}}_{\Sigma_{\theta_{\omega}}}(\tilde{T}_{k}(X_{1},\dots,X_{n})-f_{k}(\Sigma_{\theta_{\omega}}))^{2}\leq
4{\mathbb{E}}_{\Sigma_{\theta_{\omega}}}(T_{k}(X_{1},\dots,X_{n})-f_{k}(\Sigma_{\theta_{\omega}}))^{2}.$
(4.37)

Indeed, recall that
$f_{k}(\Sigma_{\theta_{\omega}})=\varphi(0)\bar{\delta}^{s}\varepsilon^{s}\omega_{k}$
and, without loss of generality, assume that $\omega_{k}=+1.$ Then
$\tilde{T}_{k}(X_{1},\dots,X_{n})-f_{k}(\Sigma_{\theta_{\omega}})=0$ when
$T_{k}(X_{1},\dots,X_{n})\geq 0$ and otherwise

$\displaystyle|\tilde{T}_{k}(X_{1},\dots,X_{n})-f_{k}(\Sigma_{\theta_{\omega}})|\leq
2|T_{k}(X_{1},\dots,X_{n})-f_{k}(\Sigma_{\theta_{\omega}})|,$

which implies the claim.

Define also $\hat{\omega}:={\rm argmin}_{\omega\in B}h(\tilde{\omega},\omega)$
and $\hat{\theta}=\theta_{\hat{\omega}}.$ By Lemma 4.9, for all $\omega\in B,$

$\displaystyle\|\hat{\theta}-\theta_{\omega}\|$
$\displaystyle\leq\frac{8}{\varphi(0)\bar{\delta}^{s}\varepsilon^{s-1}}\tau(\hat{\omega},\omega)=\frac{8}{\varphi(0)\bar{\delta}^{s}\varepsilon^{s-1}}\frac{\varphi(0)\bar{\delta}^{s}\varepsilon^{s}}{\sqrt{d}}\sqrt{h(\hat{\omega},\omega)}$
$\displaystyle\leq\frac{8}{\varphi(0)\bar{\delta}^{s}\varepsilon^{s-1}}\frac{\varphi(0)\bar{\delta}^{s}\varepsilon^{s}}{\sqrt{d}}\sqrt{h(\tilde{\omega},\omega)+h(\tilde{\omega},\hat{\omega})}\leq\frac{8}{\varphi(0)\bar{\delta}^{s}\varepsilon^{s-1}}\frac{\varphi(0)\bar{\delta}^{s}\varepsilon^{s}}{\sqrt{d}}\sqrt{2h(\tilde{\omega},\omega)}$
$\displaystyle=\frac{8\sqrt{2}}{\varphi(0)\bar{\delta}^{s}\varepsilon^{s-1}}\tau(\tilde{\omega},\omega).$

Therefore,

$\displaystyle{\mathbb{E}}_{\theta_{\omega}}\|\hat{\theta}-\theta_{\omega}\|^{2}$
$\displaystyle\leq\frac{2\times
8^{2}}{\varphi^{2}(0)\bar{\delta}^{2s}\varepsilon^{2(s-1)}}{\mathbb{E}}_{\theta_{\omega}}\tau^{2}(\tilde{\omega},\omega)=\frac{2\times
8^{2}}{\varphi^{2}(0)\bar{\delta}^{2s}\varepsilon^{2(s-1)}}\frac{1}{d}\sum_{k=1}^{d}{\mathbb{E}}_{\theta_{\omega}}\Bigl{(}f_{k}(\Sigma_{\theta_{\tilde{\omega}}})-f_{k}(\Sigma_{\theta_{\omega}})\Bigr{)}^{2}$
$\displaystyle=\frac{2\times
8^{2}}{\varphi^{2}(0)\bar{\delta}^{2s}\varepsilon^{2(s-1)}}\frac{1}{d}\sum_{k=1}^{d}{\mathbb{E}}_{\Sigma_{\theta_{\omega}}}\Bigl{(}\tilde{T}_{k}(X_{1},\dots,X_{n})-f_{k}(\Sigma_{\theta_{\omega}})\Bigr{)}^{2}$
$\displaystyle\leq\frac{8^{3}}{\varphi^{2}(0)\bar{\delta}^{2s}\varepsilon^{2(s-1)}}\frac{1}{d}\sum_{k=1}^{d}{\mathbb{E}}_{\Sigma_{\theta_{\omega}}}\Bigl{(}T_{k}(X_{1},\dots,X_{n})-f_{k}(\Sigma_{\theta_{\omega}})\Bigr{)}^{2},$

where we also used bound (4.37). This implies that

$\displaystyle\max_{\omega\in
B}{\mathbb{E}}_{\theta_{\omega}}\|\hat{\theta}-\theta_{\omega}\|^{2}\leq\frac{8^{3}}{\varphi^{2}(0)\bar{\delta}^{2s}\varepsilon^{2(s-1)}}\frac{1}{d}\sum_{k=1}^{d}\max_{\omega\in
B}{\mathbb{E}}_{\Sigma_{\theta_{\omega}}}\Bigl{(}T_{k}(X_{1},\dots,X_{n})-f_{k}(\Sigma_{\theta_{\omega}})\Bigr{)}^{2}.$

Assuming that

$\displaystyle\max_{1\leq k\leq d}\inf_{T}\max_{\omega\in
B}{\mathbb{E}}_{\Sigma_{\theta_{\omega}}}\Bigl{(}T(X_{1},\dots,X_{n})-f_{k}(\Sigma_{\theta_{\omega}})\Bigr{)}^{2}<\tau^{2},$

one can find, for each $k=1,\dots,d,$ an estimator $T_{k}(X_{1},\dots,X_{n})$
such that

$\displaystyle\max_{\omega\in
B}{\mathbb{E}}_{\Sigma_{\theta_{\omega}}}\Bigl{(}T_{k}(X_{1},\dots,X_{n})-f_{k}(\Sigma_{\theta_{\omega}})\Bigr{)}^{2}<\tau^{2},$

which implies the existence of $\hat{\theta}$ such that

$\displaystyle\max_{\omega\in
B}{\mathbb{E}}_{\theta_{\omega}}\|\hat{\theta}-\theta_{\omega}\|^{2}\leq\frac{8^{3}\tau^{2}}{\varphi^{2}(0)\bar{\delta}^{2s}\varepsilon^{2(s-1)}}.$

Note that, for $\lambda=\gamma_{1}a,\mu=\gamma_{2}a,$
$(\lambda-\mu)(\mu^{-1}-\lambda^{-1})=(\gamma_{1}-\gamma_{2})(\gamma_{2}^{-1}-\gamma_{1}^{-1}),$
so condition (4.35) of Lemma 4.6 becomes

$\displaystyle\varepsilon\leq\frac{c}{\sqrt{(\gamma_{1}-\gamma_{2})(\gamma_{2}^{-1}-\gamma_{1}^{-1})}}\sqrt{\frac{d}{n}}:=\varepsilon(\gamma_{1},\gamma_{2})$
(4.38)

with a sufficiently small constant $c>0.$ In view of Lemma 4.6, it follows
that

$\displaystyle\frac{\tau^{2}}{\varphi^{2}(0)\bar{\delta}^{2s}\varepsilon^{2(s-1)}}\gtrsim\varepsilon^{2}$

for an arbitrary $\varepsilon$ satisfying the condition (4.38). We set
$\varepsilon:=\varepsilon(\gamma_{1},\gamma_{2})$ and also observe that
$\bar{\delta}=\frac{(\gamma_{1}-\gamma_{2})a}{8}\wedge
1\geq\frac{\gamma_{1}-\gamma_{2}}{8}(a\wedge 1),$ which allows us to conclude
that

$\displaystyle\tau^{2}\gtrsim\bar{\delta}^{2s}\varepsilon^{2s}\gtrsim\gamma_{1}^{s}\gamma_{2}^{s}(a\wedge
1)^{2s}\Bigl{(}\frac{d}{n}\Bigr{)}^{s}.$

In other words,

$\displaystyle\max_{1\leq k\leq d}\inf_{T}\max_{\omega\in
B}{\mathbb{E}}_{\Sigma_{\theta_{\omega}}}\Bigl{(}T(X_{1},\dots,X_{n})-f_{k}(\Sigma_{\theta_{\omega}})\Bigr{)}^{2}\gtrsim\gamma_{1}^{s}\gamma_{2}^{s}(a\wedge
1)^{2s}\Bigl{(}\frac{d}{n}\Bigr{)}^{s}.$

Finally, by Lemma 4.7 for all $\omega\in B,$
$\|\Sigma_{\theta_{\omega}}-\Sigma_{0}\|\leq
5\sqrt{2}(\gamma_{1}-\gamma_{2})a\varepsilon.$ For
$\varepsilon=\varepsilon(\gamma_{1},\gamma_{2}),$ we get
$\|\Sigma_{\theta_{\omega}}-\Sigma_{0}\|<\delta$ provided that
$\delta>5\sqrt{2}c\sqrt{\gamma_{1}\gamma_{2}}a\sqrt{\frac{d}{n}}.$ Taking also
into account that $\|f_{k}\|_{C^{s}(U)}\leq 1,k=1,\dots,d,$ we can conclude
that

$\displaystyle\sup_{\|f\|_{C^{s}(U)}\leq 1}\sup_{{\rm
Im}(\Sigma)=L,\|\Sigma-\Sigma_{0}\|<\delta}{\mathbb{E}}_{\Sigma}\Bigl{(}T(X_{1},\dots,X_{n})-f(\Sigma)\Bigr{)}^{2}\gtrsim\gamma_{1}^{s}\gamma_{2}^{s}(a\wedge
1)^{2s}\Bigl{(}\frac{d}{n}\Bigr{)}^{s}.$

for a sufficiently small constant $c^{\prime}>0.$ It remains to combine the
last bound with (4.3) and to recall that $d=[r],$ to get

$\displaystyle\sup_{\|f\|_{C^{s}(U)}\leq
1}\sup_{\Sigma\in{\mathcal{S}}(a,r),\|\Sigma-\Sigma_{0}\|<\delta}\Bigl{\|}T(X_{1},\dots,X_{n})-f(\Sigma)\Bigr{\|}_{L_{2}({\mathbb{P}}_{\Sigma})}\gtrsim\gamma_{1}^{s/2}\gamma_{2}^{s/2}(a\wedge
1)^{s}\Bigl{(}\sqrt{\frac{r}{n}}\Bigr{)}^{s}.$

∎

The claim of Theorem 1.9 immediately follows by combining propositions 4.1 and
4.2.

## References

  * [1] P. Bickel and Y. Ritov.  Estimating integrated square density derivatives: sharp best order of convergence estimates.  Sankhya, 1988, 50, 381–393.
  * [2] L. Birgé and P. Massart.  Estimation of integral functionals of a density.  Annals of Statistics, 1995, 23, 11-29.
  * [3] P. Hall and M.A. Martin.  On Bootstrap Resampling and Iteration.  Biometrika, 1988, 75, 4, 661–671.
  * [4] I. A. Ibragimov and R.Z. Khasminskii.  Statistical Estimation: Asymptotic Theory.  Springer-Verlag, New York, 1981.
  * [5] I.A. Ibragimov, A.S. Nemirovski and R.Z. Khasminskii.  Some problems of nonparametric estimation in Gaussian white noise.  Theory of Probab. and Appl., 1987, 31, 391–406.
  * [6] J. Jiao and Y. Han.  Bias correction with Jackknife, Bootstrap and Taylor Series. IEEE Trans. on Information Theory, 2020, 66, 7, 4392–4418.
  * [7] V. Koltchinskii.  Asymptotically Efficient Estimation of Smooth Functionals of Covariance Operators.  J. European Mathematical Society, 2021, 23, 3, 765–843.
  * [8] V. Koltchinskii.  Asymptotic Efficiency in High-Dimensional Covariance Estimation.  Proc. ICM 2018, Rio de Janeiro, 2018, vol. 3, 2891–2912.
  * [9] V. Koltchinskii. Estimation of smooth functionals in high-dimensional models: bootstrap chains and Gaussian approximation. 2021\. arXiv:2011.03789.
  * [10] V. Koltchinskii and K. Lounici.  Asymptotics and concentration bounds for bilinear forms of spectral projectors of sample covariance. Ann. Inst. H. Poincaré Probab. Statist., 2016, 52, 4, 1976–2013.
  * [11] V. Koltchinskii and K. Lounici.  Concentration inequalities and moment bounds for sample covariance operators.  Bernoulli, 2017, 23, 1, 110–133.
  * [12] V. Koltchinskii and K. Lounici.  Normal approximation and concentration of spectral projectors of sample covariance.  Annals of Statistics, 2017, 45, 1, 121–157.
  * [13] V. Koltchinskii, M. Löffler and R. Nickl.  Efficient Estimation of Linear Functionals of Principal Components.  Annals of Statistics, 2020, 48, 1, 464–490.
  * [14] V. Koltchinskii and M. Wahl. Functional estimation in log-concave location families. 2021\. arXiv:2108.00263.
  * [15] V. Koltchinskii and M. Zhilova.  Efficient estimation of smooth functionals in Gaussian shift models.  Ann. Inst. H. Poincaré - Probab. et Statist., 2021, 57, 1, 351–386.
  * [16] V. Koltchinskii and M. Zhilova.  Estimation of Smooth Functionals in Normal Models: Bias Reduction and Asymptotic Efficiency, Annals of Statistics, 2021, 49, 5, 2577–2610.
  * [17] V. Koltchinskii and M. Zhilova. Estimation of smooth functionals of location parameter in Gaussian and Poincaré random shift models. Sankhya A, 2021, 83, 569–596.
  * [18] S. Kwapien and B. Szymanski. Some remarks on Gaussian measures in Banach spaces.  Probability and Mathematical Statistics, 1980, 1, 1, 59–65.
  * [19] B. Levit.  On the efficiency of a class of non-parametric estimates.  Theory of Prob. and applications, 1975, 20(4), 723–740.
  * [20] B. Levit.  Asymptotically efficient estimation of nonlinear functionals.  Probl. Peredachi Inf. (Problems of Information Transmission), 1978, 14(3), 65–72.
  * [21] O. Lepski, A. Nemirovski and V. Spokoiny.  On estimation of the $L_{r}$ norm of a regression function.  Probab. Theory Relat. Fields, 1999, 113, 221–253.
  * [22] A. Nemirovski.  On necessary conditions for the efficient estimation of functionals of a nonparametric signal which is observed in white noise.  Theory of Probab. and Appl., 1990, 35, 94–103.
  * [23] A. Nemirovski.  Topics in Non-parametric Statistics.  Ecole d’Ete de Probabilités de Saint-Flour.  Lecture Notes in Mathematics, v. 1738, Springer, New York, 2000.
  * [24] E. Rio. Upper bounds for minimal distances in the central limit theorem. Ann. Inst. Henry Poincarè Probab. Stat., 2009, 45, 3, 802–817.
  * [25] J. Robins, L. Li, E. Tchetgen and A. van der Vaart.  Higher order influence functions and minimax estimation of nonlinear functionals.  IMS Collections Probability and Statistics: Essays in Honor of David. A. Freedman, 2008, vol. 2, 335-421.
  * [26] A.B. Tsybakov.  Introduction to Nonparametric Estimation. Springer, 2009.
  * [27] H. Triebel.  Theory of Function Spaces, Springer, 1983.
  * [28] F. Zhou, P. Li and C.-H. Zhang.  High-Order Statistical Functional Expansion and Its Application To Some Nonsmooth Problems. 2021, arXiv:2112.15591.