Primordial Features
from Linear to Nonlinear Scales

Florian Beutler,1,2 Matteo Biagetti,3 Daniel Green,4

Anže Slosar5 and Benjamin Wallisch6,7,3,4

1 Institute of Cosmology & Gravitation, University of Portsmouth, Portsmouth,
PO1 3FX, UK

2 Lawrence Berkeley National Laboratory, 1 Cyclotron Road, Berkeley, CA 94720,
USA

3 Institute for Theoretical Physics, University of Amsterdam, 1098 XH
Amsterdam, NL

4 Department of Physics, University of California, San Diego, La Jolla, CA
92093, USA

5 Physics Department, Brookhaven National Laboratory, Upton, NY 11973, USA

6 School of Natural Sciences, Institute for Advanced Study, Princeton, NJ
08540, USA

7 DAMTP, University of Cambridge, Cambridge CB3 0WA, UK

Abstract
Sharp features in the primordial power spectrum are a powerful window into the
inflationary epoch. To date, the cosmic microwave background (CMB) has offered
the most sensitive avenue to search for these signatures. In this paper, we
demonstrate the power of large-scale structure observations to surpass the CMB
as a probe of primordial features. We show that the signatures in galaxy
surveys can be separated from the broadband power spectrum and are as robust
to the nonlinear evolution of matter as the standard baryon acoustic
oscillations. As a result, analyses can exploit a significant range of scales
beyond the linear regime available in the datasets. We develop a feature
search for large-scale structure, apply it to the final data release of the
Baryon Oscillation Spectroscopic Survey and find new bounds on oscillatory
features that exceed the sensitivity of Planck for a significant range of
frequencies. Moreover, we forecast that the next generation of galaxy surveys,
such as DESI and Euclid, will be able to improve current constraints by up to
an order of magnitude over an expanded frequency range.

###### Contents

  1. 1 Introduction
  2. 2 Primordial Features and Galaxy Surveys
    1. 2.1 Oscillatory Features in the Primordial Spectrum
    2. 2.2 Robustness of Features to Small-Scale Nonlinearities
    3. 2.3 Damping from Large-Scale Nonlinearities
  3. 3 Feature Search in Large-Scale Structure
    1. 3.1 BOSS DR12 Dataset and Analysis Pipeline
    2. 3.2 Forecasting Methodology for BOSS
    3. 3.3 Validation on Mock Catalogs and in Forecasts
  4. 4 First Large-Scale Structure Constraints
    1. 4.1 Limits on Features from BOSS DR12
    2. 4.2 Comparison with Planck CMB Bounds
    3. 4.3 Future LSS and CMB Constraints
  5. 5 Conclusions
  6. A Nonlinear Damping of Logarithmic Features
    1. A.1 Perturbative Treatment
    2. A.2 Infrared Resummation
  7. B Large-Scale Structure Forecasts
    1. B.1 Forecasting with the Wiggle Spectrum
    2. B.2 Experimental Specifications
    3. B.3 Additional Tests of the Pipeline
    4. B.4 Forecasts for Future LSS Surveys
  8. C Further Details on the BOSS Analysis
    1. C.1 From Posteriors to Upper Limits
    2. C.2 Upper Limits or Detections?
  9. D CMB Analysis and Forecasts
    1. D.1 Analysis of Planck Data
    2. D.2 Transfer of Feature Power
    3. D.3 Joint CMB and LSS Analysis
    4. D.4 Forecasts for Future CMB Surveys

##  1 Introduction

Characterizing the nature of inflation is one of the major challenges in
cosmology. While current data are compatible with the simplest incarnation of
inflation, a single weakly-coupled scalar field on a very flat potential, the
space of possibilities for inflation are vast and should ultimately be settled
by data. An appealing aspect of the simplest models is that the observed
(near) scale invariance of the power spectrum of fluctuations is easily
explained by the flatness of this potential [1]. However, attempts to realize
inflation from a more fundamental starting point can lead to much more
complicated models where many interconnected pieces are needed to achieve
scale invariance [2, 3]. This dichotomy between the simplest models and
ultraviolet-complete examples originates from quantum gravity itself:
sufficiently flat potentials can be engineered using symmetries, but quantum
gravity famously abhors them [4, 5]. As a result, a wide variety of models
gives rise to non-trivial deviations from canonical slow-roll and scale
invariance (i.e. features) [3, 6]. Other mechanisms avoid this picture
altogether by invoking non-trivial interactions that can lead to non-Gaussian
$n$-point correlation functions (primordial non-Gaussianity) [7].

Our ability to test these ideas directly with data relies on separating the
primordial signatures of interest from a broad range of processes at late
times. This challenge is particularly acute for constraints on inflation from
large-scale structure (LSS) surveys since the observed objects, i.e. galaxies,
owe their existence to the nonlinear gravitational evolution of matter
fluctuations in the late universe. While both current and future galaxy
surveys have the raw statistical power to compete with other cosmological
probes, such as the cosmic microwave background (CMB), the useful information
is greatly diminished if we restrict our analyses to modes that are
sufficiently linear to use forward modeling in order to isolate the primordial
information. There have been steady improvements in this direction, but
eventually these techniques are expected to be limited by the complexity of
(astro)physics at short distances [8, 9].

An alternate and already successful approach is to look for special
observables that are (at least partially) immune to the complications
presented by LSS data. The best-known example of this type are the baryon
acoustic oscillations (BAO). Although they manifest themselves as an
oscillation in the power spectrum in a range of wavenumbers sensitive to
nonlinearities, it is more usefully understood as a sharp peak in the two-
point correlation function at the size of the sound horizon, which is a scale
that is much larger than the scale where nonlinear evolution dominates [10].
More recently, a constant phase of the baryon acoustic oscillations was shown
to be immune to nonlinear evolution [11]. These are useful examples as they
show that smooth and oscillatory power spectra are not sensitive to the same
nonlinear effects.

Nevertheless, it has proven challenging to constrain the most common
inflationary parameters (e.g. the scalar spectral index $n_{\mathrm{s}}$ or
its running $\alpha_{\mathrm{s}}$) by LSS. Changes to these parameters
typically lead to smooth variations in the power spectrum (as a function of
wavenumber $k$) and are therefore degenerate with other contributions such as
galaxy biasing and baryonic effects. Furthermore, most inflationary
observables get their constraining power from the smallest physical scales
accessible in a given survey where gravitational nonlinearities dominate. With
current data, large-scale structure is most competitive with the CMB as a
probe of inflation for constraints on local primordial non-Gaussianity [12,
13, 14]. In this case, the non-Gaussian signature in the initial conditions
manifests itself in the biasing of galaxies on the largest scales where
nonlinearities are negligible [15] (see [16] for a recent review). Future
surveys will search these large scales with increasing sensitivity and have
the potential to ultimately exceed the CMB [17, 18].

In this paper, we demonstrate that features in the primordial spectra, much
like the standard BAO signal itself, are immune to short-distance nonlinear
processes of the late universe and the effects of large-scale bulk flows can
be captured analytically. We can therefore test these models with the full
statistical power of LSS surveys. Phenomenologically, primordial features are
generally characterized by significant deviations from scale invariance over a
narrow range of scales, usually in the power spectrum. The most canonical
examples, shown in Fig. 1,

Figure 1: Illustration of linear and logarithmic oscillations in the
primordial power spectrum $\mathcal{P}(k)$. The employed frequencies
$\omega_{X}$ are comparable in the range of scales
$k=$0.10.2\text{\,}\mathrm{\mathit{h}}\text{\,}{\mathrm{Mpc}}^{-1}$$ where
BOSS has the largest signal-to-noise ratio. We see that the effective
frequency of the logarithmic oscillations decreases as we go to larger
wavenumbers $k$.

are written as oscillations in either $k$ or $\log k$ [6]:

$P_{\zeta}(k)=P^{s}(k)\,\big{[}1+A_{\mathrm{lin}}\sin(\omega_{\mathrm{lin}}k+\varphi_{\mathrm{lin}})+A_{\mathrm{log}}\sin(\omega_{\mathrm{log}}\log(k/k_{\star})+\varphi_{\mathrm{log}})\big{]}\,,$
(1.1)

where $P_{\zeta}(k)$ is the primordial power spectrum of adiabatic density
fluctuations ($\zeta$) and $P^{s}(k)$ is a smooth function of $k$. For a
linear oscillation ($A_{\mathrm{lin}}$), the insensitivity to nonlinear
effects is identical to the case of the baryon acoustic oscillations: if we
Fourier transform the signal, the linear oscillation is a sharp peak at the
scale $\omega_{\mathrm{lin}}$ in the two-point correlation function. For
sufficiently large $\omega_{\mathrm{lin}}$, one is effectively looking for a
second BAO peak. For logarithmic oscillations ($A_{\mathrm{log}}$), there is
not a simple description in terms of scales, but we will show that local
nonlinear evolution is incapable of producing the same oscillation for
sufficiently large frequencies $\omega_{\mathrm{log}}$.

While the analogy with the BAO signal is useful, it is worth observing that
the primordial features arise directly in the initial conditions of the dark
matter and baryons, and are suppressed only by the amplitude of the
oscillation. By contrast, the baryon acoustic oscillations themselves are a
consequence of the physics of baryons. This means that their impact on LSS
data is suppressed by $\omega_{b}/\omega_{m}$ and the growth of dark matter
density fluctuations prior to recombination. As a result, the amplitude of the
BAO spectrum is suppressed relative to a primordial feature and is roughly
equivalent to a linear feature amplitude of $A_{\mathrm{lin}}=0.05$. The fact
that we do not see an additional oscillation beyond the BAO signal by eye
already suggests that $A_{\mathrm{lin}}\ll 0.1$ without any analysis.

In addition to being protected quantities in LSS, primordial features of this
form are a well-motivated probe of the early universe in their own right [19,
20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 6]. They notably arise
in axion monodromy inflation [34] as a direct consequence of the fundamental
symmetry structure needed to produce large-field inflation. The periodic,
nonperturbative potential generated for axions gives rise to an oscillation in
the inflationary potential and manifests itself as a logarithmic oscillation
in the power spectrum [26]. In addition, there are a number of scenarios where
particles are excited from the vacuum at a specific time through non-adiabatic
evolution and can give rise to linear oscillations. Moreover, most features in
the power spectrum can be efficiently decomposed in this basis of functions
which therefore captures large parts of model space.

The outline of the paper is as follows. In Section 2, we show that features
are robust to small-scale nonlinearities and compute the nonlinear damping
effect due to long-wavelength modes. In Section 3, we introduce our new
analysis to search for these features in LSS data and verify that these
oscillatory signals can be reliably constrained. In Section 4, we apply this
pipeline to BOSS DR12 data and present a new constraint that exceeds Planck
over a significant range of frequencies. Moreover, we forecast the sensitivity
of future observational surveys (cf. e.g. [35, 36, 37, 38, 39, 40, 41, 42, 43,
44, 45, 46] for previous work). We conclude in Section 5. Additional details
on the theoretical calculation, the employed forecasts, and the performed LSS
and CMB analyses are provided in a set of four appendices.

##  2 Primordial Features and Galaxy Surveys

In this section, we determine how a primordial feature will appear in the
nonlinear (low-redshift) universe. We first characterize the signals in the
linear matter power spectrum. Then, we will use a linear-response argument to
show that nonlinear evolution on small scales does not change the amplitude of
the feature in the nonlinear power spectrum. Finally, we will use infrared
resummation to determine the nonlinear damping of the features from large-
scale modes.

###  2.1 Oscillatory Features in the Primordial Spectrum

The physics of inflation determines the primordial power spectrum
$P_{\zeta}(k)$. In the simplest versions of inflation, this power spectrum
arises from the freeze-out of quantum mechanical fluctuations when the
physical wavelength reaches the Hubble radius, $k=aH(t_{\star})$, where $k$ is
the constant comoving wavenumber, $a(t)$ is the scale factor and $H(t)$ is the
Hubble parameter. (This equation defines the freeze-out time $t_{\star}$.) As
a result, the amplitude of fluctuations for a comoving wavenumber $k$ is
determined by the physics of inflation around the time $t_{\star}$.
Furthermore, the near scale invariance of the resulting spectrum is a
consequence of the weak time dependence in the evolution of the perturbations
during inflation.

Scale-dependent features in the primordial spectrum therefore arise from
strongly time-dependent physics during inflation.111We focus on the
inflationary origin of features, but note that they might also appear in
alternatives to inflation [47]. This may be due to sharp features in the
underlying potential for a scalar field, or special locations in the field
space of the inflaton where other particles become light and can be excited
from the vacuum. Despite the wide range of possibilities, the signatures do
not significantly depend on the details of the model since the nature of the
time dependence controls the deviations from scale invariance.

We will assume that the smooth spectrum $P^{s}(k)$ of (1.1) is the almost
scale-invariant power spectrum of curvature perturbations in vanilla models of
inflation,

$P_{\zeta,0}(k)=\frac{2\pi^{2}}{k^{3}}\mathcal{P}_{\zeta,0}(k)=\frac{2\pi^{2}A_{\mathrm{s}}}{k^{3}}\left(\frac{k}{k_{\star}}\right)^{\\!n_{\mathrm{s}}-1}\,,$
(2.1)

where $A_{\mathrm{s}}$ and $n_{\mathrm{s}}$ are the scalar amplitude and
spectral index at the pivot scale $k_{\star}$, which we generally take to be
$k_{\star}=$0.05\text{\,}{\mathrm{Mpc}}^{-1}$$. We then write the full power
spectrum $P_{\zeta}(k)$ including the contribution from features, $\delta
P_{\zeta}(k)$, as

$P_{\zeta}(k)=P_{\zeta,0}(k)\left[1+\delta P_{\zeta}(k)\right].$ (2.2)

As suggested in (1.1), we will consider oscillatory features with linearly-
spaced oscillations, $\delta P_{\zeta}^{\mathrm{lin}}$, and logarithmically-
spaced oscillations, $\delta P_{\zeta}^{\mathrm{log}}$. We phenomenologically
parameterize the former, which we refer to as linear features, as
follows:222Note that the feature frequency is sometimes defined in the
literature as $\tilde{\omega}_{\mathrm{lin}}=\omega_{\mathrm{lin}}/2$ or with
respect to a pivot scale $k_{\star}$. More generally, we highlight that
$\omega_{\mathrm{lin}}$ is a frequency in Fourier space which corresponds to a
physical scale in real space.

$\displaystyle\delta P_{\zeta}^{\mathrm{lin}}(k)$
$\displaystyle=A_{\mathrm{lin}}^{\mathrm{sin}}\sin\\!\left(\omega_{\mathrm{lin}}k\right)+A_{\mathrm{lin}}^{\mathrm{cos}}\cos\\!\left(\omega_{\mathrm{lin}}k\right)$
(2.3)
$\displaystyle=A_{\mathrm{lin}}\sin\\!\left(\omega_{\mathrm{lin}}k+\varphi_{\mathrm{lin}}\right),$

with the feature frequency $\omega_{\mathrm{lin}}$, and the amplitudes of the
sine and cosine contributions $A_{\mathrm{lin}}^{\mathrm{sin}}$ and
$A_{\mathrm{lin}}^{\mathrm{cos}}$, respectively, or the overall feature
amplitude $A_{\mathrm{lin}}$ and corresponding phase $\varphi_{\mathrm{lin}}$.
The so-called logarithmic features are similarly defined as

$\displaystyle\delta P_{\zeta}^{\mathrm{log}}(k)$
$\displaystyle=A_{\mathrm{log}}^{\mathrm{sin}}\sin\\!\left[\omega_{\mathrm{log}}\log(k/k_{\star})\right]+A_{\mathrm{log}}^{\mathrm{cos}}\cos\\!\left[\omega_{\mathrm{log}}\log(k/k_{\star})\right]$
(2.4)
$\displaystyle=A_{\mathrm{log}}\sin\\!\left[\omega_{\mathrm{log}}\log(k/k_{\star})+\varphi_{\mathrm{log}}\right],$

with the feature frequency $\omega_{\mathrm{log}}$, and the amplitudes of the
sine and cosine contributions $A_{\mathrm{log}}^{\mathrm{sin}}$ and
$A_{\mathrm{log}}^{\mathrm{cos}}$, respectively, or the overall feature
amplitude $A_{\mathrm{log}}$ and corresponding phase $\varphi_{\mathrm{log}}$.
We refer to the parameterization in terms of two amplitudes as ‘amplitude
parameterization’ and the one in terms of the overall amplitude and a phase as
‘phase parametrization’. We note that it has been customary in the literature
to define the linear feature frequency $\omega_{\mathrm{lin}}$ as a
dimensionful quantity in units of $\mathrm{Mpc}$, whereas the logarithmic
feature frequency $\omega_{\mathrm{log}}$ is dimensionless. In addition, we
remark that the feature amplitudes give the contribution relative to the
standard power spectrum $P_{\zeta,0}$ and the feature contribution $\delta
P_{\zeta}(k)$ oscillates around zero by construction.

The information in the primordial power spectrum is transferred to the matter
power spectrum, as illustrated in Fig. 2,

Figure 2: Imprint of primordial features in several large-scale structure
observables at redshift $z=0$: the matter power spectrum $P(k)$ (top) and the
relative wiggle spectrum $P^{\mathrm{w}}(k)/P^{\mathrm{nw}}(k)$ (middle) in
Fourier space, and the (rescaled) two-point correlation function $\xi(r)$ in
real space (bottom). We compare a featureless model (gray) to scenarios
involving a linear feature (left) and a logarithmic feature (right) with
$A^{\mathrm{sin}}_{X}=0.05$ and $A^{\mathrm{cos}}_{X}=0$,
$X=\mathrm{lin},\mathrm{log}$. In addition to the predictions in linear theory
(dashed), we also show the observables including nonlinear corrections (solid)
from a theoretical calculation of the damping.

through the usual linear evolution from initial conditions,

$P(k,z)=k^{4}\,T(k)^{2}\,D(z)^{2}\,P_{\zeta}(k)\,,$ (2.5)

where $D(z)$ is the linear growth rate and $T(k)$ the transfer function. For
simplicity, we omit the redshift dependence of the matter power spectrum,
$P(k)=P(k,z)$, from now on. For large enough feature frequencies, these
oscillations can be distinguished from the broadband shape of the power
spectrum, similar to the baryon acoustic oscillations. It is therefore natural
to constrain the feature models as contributions to the BAO spectrum, which is
why we split the power spectrum into a smooth (‘no-wiggle’) part and an
oscillatory (‘wiggle’) part,

$P(k)\equiv P^{\mathrm{nw}}(k)+P^{\mathrm{w}}(k)\,.$ (2.6)

Since primordial features with large enough frequencies are only contained in
the second term, it is useful to further decompose the wiggle spectrum as
follows:

$P^{\mathrm{w}}(k)\equiv
P^{\mathrm{w}}_{\mathrm{BAO}}(k)+P^{\mathrm{w}}_{X}(k)+P^{\mathrm{w}}_{\mathrm{BAO}}(k)\,\delta
P^{X}_{\zeta}(k)\,,$ (2.7)

where $P^{\mathrm{w}}_{\mathrm{BAO}}(k)$ is the standard BAO spectrum in a
featureless $\Lambda$CDM cosmology, the auto-spectrum of possible primordial
features is

$P^{\mathrm{w}}_{X}(k)=P^{\mathrm{nw}}(k)\delta P_{\zeta}^{X}(k)\,,$ (2.8)

with $X=\mathrm{lin},\mathrm{log}$, and the third term is the BAO-feature
cross-correlation power spectrum. Since the BAO signal itself is only a small
(five-percent) contribution to $P(k)$, we will be able to generally neglect
this cross-spectrum term in our theoretical considerations. Given that (2.6)
is the linear matter power spectrum, we will show in the following two
sections that small-scale modes do not affect $P^{\mathrm{w}}(k)$ for large
enough feature frequencies (§2.2), but that each of its oscillatory components
in (2.7) is affected by gravitational large-scale nonlinearities and
exponentially damped (§2.3).

###  2.2 Robustness of Features to Small-Scale Nonlinearities

The smallest scales in an LSS survey carry most of the statistical power, but
are also the most prone to nonlinear corrections, including galaxy bias and
baryonic effects. For this reason, a typical analysis might cut at
$k=$0.10.2\text{\,}\mathrm{\mathit{h}}\text{\,}{\mathrm{Mpc}}^{-1}$$ to avoid
the complications of modeling and marginalizing over these effects. In the
case of the baryon acoustic oscillations, it has long been known that they are
robust to those effects which change the power spectrum only by a smooth
window function and that it is possible to aggressively marginalize over these
smooth corrections without losing any information. In the following, we show
that high-frequency oscillatory features in the power spectrum are protected
from small-scale nonlinearities in precisely the same fashion as the BAO
signal. We first give an intuitive argument that this is indeed the case which
we then confirm more rigorously.

####  2.2.1 Intuitive argument

The power spectrum of linear oscillations $P^{\mathrm{w}}_{\mathrm{lin}}(k)$
as defined in (2.3) and (2.8) is the same as the approximate form of the BAO
signal, where the linear feature frequency $\omega_{\mathrm{lin}}$ corresponds
to the sound horizon $r_{s}$ and the feature phase $\varphi_{\mathrm{lin}}$ is
the phase shift due to free-streaming neutrinos [11, 48, 49].333The real form
of the BAO spectrum has an additional $k$-dependence in the phase and
amplitude [48] which are however unimportant for this analogy. It is a long-
established fact that the baryon acoustic oscillations are essentially immune
from nonlinearities on small scales because the sound horizon is a large-
distance scale. The same will be true for any linear oscillation with
$\omega_{\mathrm{lin}}$ larger than the BAO scale $r_{s}$. Smaller frequencies
will remain immune to nonlinearities provided that they are sufficiently large
compared to the nonlinear scale. Conservatively, we will use
$\omega_{\mathrm{lin}}\gtrsim$75\text{\,}\mathrm{Mpc}$$ as a reference value.
As illustrated in Fig. 3,

Figure 3: Separation of zeros (or, equivalently, peaks and troughs) in the
spectrum of logarithmically-spaced oscillations (2.4) as a function of
wavenumber $k$ for a range of frequencies $\omega_{\mathrm{log}}$. For
comparison, the separations in the standard BAO spectrum and for linear
oscillations with $\omega_{\mathrm{lin}}=$75\text{\,}\mathrm{Mpc}$$ and
$\omega_{\mathrm{lin}}=\omega_{\mathrm{Ny}}\approx$929\text{\,}\mathrm{Mpc}$$,
which is the Nyquist frequency in our BOSS analysis below, are also shown. The
dotted and dashed gray lines indicate the nonlinear scale
$k_{\mathrm{nl}}\approx$0.15\text{\,}\mathrm{\mathit{h}}\text{\,}{\mathrm{Mpc}}^{-1}$$
and $k_{\omega}=k$, respectively. The former marks the approximate wavenumber
where nonlinearities are expected to become large, $k\gtrsim k_{\mathrm{nl}}$,
while these nonlinearities do not alter features with $k_{\omega}\lesssim k$.
The zero separations are denoted by $k_{\omega}/2$ since they are equivalent
to half of a feature oscillation period.

this frequency is consistent with our more rigorous constraint on the scales
where small-scale nonlinearities are not important (cf. §2.2.2). More
recently, it has been shown that the phase $\varphi_{\mathrm{lin}}$ is also
protected [11] which implies that this argument will hold for both the sine
and cosine contributions.

A priori, it is less clear that logarithmic oscillations will share the same
nice properties. Unlike the linear frequency $\omega_{\mathrm{lin}}$, the
frequency $\omega_{\mathrm{log}}$ of the logarithmically-spaced oscillations
is a dimensionless parameter and, therefore, does not refer to any fixed
physical scale in either Fourier or configuration space. As a result, the
signal must appear on all scales in both descriptions and, consequently, is
not clearly distinct from nonlinear effects. Fortunately, in practice, our
information only comes from a limited range of scales which is given by
$k\sim$0.10.3\text{\,}\mathrm{\mathit{h}}\text{\,}{\mathrm{Mpc}}^{-1}$$ for
BOSS. We can always (Fourier) decompose any feature into a sum over linear
oscillations. For sufficiently large $\omega_{\mathrm{log}}$, the signal in
this range of wavenumbers is reliably reproduced keeping only the linear
oscillations with frequencies large enough to be protected by our previous
argument. We illustrate this argument in Fig. 3, where we show the peak-trough
separation for the logarithmically-oscillating power spectrum
$P^{\mathrm{w}}_{\mathrm{log}}(k)$ defined in (2.4) and (2.8) as a function of
$k$ in our range of interest. For $\omega_{\mathrm{log}}\gtrsim 20$, the peak-
trough separation is at the level of the BAO spectrum, while nonlinearities
are expected to affect the power spectrum at smaller scales. For the purposes
of an analysis of BOSS data,444In future surveys covering a larger range of
$k$, more work may be required to account for the fact that the effective
frequency evolves with $k$ and may not be protected over the full range of
scales. Moreover, these surveys will generally require more theoretical
control since they will be sensitive to smaller signals and, therefore,
subleading effects. the logarithmic features therefore do not present any
significant complications.

####  2.2.2 Rigorous argument

We can rigorously derive these results using the techniques developed in [11].
Since the amplitude is known to be small (percent level at best), we can find
the nonlinear density fluctuations as the linear response in the feature
amplitude ($A_{\mathrm{lin}}$ or $A_{\mathrm{log}}$). On general grounds, this
response must take the form [11]

$\delta^{\mathrm{w}}({\vec{x}},\tau)=\int\\!\mathrm{d}^{3}x^{\prime}\,G({\color[rgb]{0.92,0.,0.}{\vec{x}}},{\color[rgb]{0.25,0.41,0.88}{\vec{x}}-{\vec{x}}\hskip
1.0pt^{\prime}};\tau)\,\delta_{\mathrm{in}}^{\mathrm{w}}({\vec{x}}\hskip
1.0pt^{\prime})\,,$ (2.9)

where
$G({\color[rgb]{0.92,0.,0.}{\vec{x}}},{\color[rgb]{0.25,0.41,0.88}{\vec{x}}-{\vec{x}}\hskip
1.0pt^{\prime}};\tau)$ is the response function and
$\delta_{\mathrm{in}}^{\mathrm{w}}({\vec{x}}^{\prime})$ is the contribution to
the initial density contrast at linear order in $A_{X}$,
$X=\mathrm{lin},\mathrm{log}$. The response function is independent of the
oscillatory signal and has previously been studied for more general
applications [50, 51, 52]. The first index in red,
${\color[rgb]{0.92,0.,0.}{\vec{x}}}$, arises from the underlying inhomogeneity
of the universe and would be absent in a translation-invariant system. The
second index in blue, ${\color[rgb]{0.25,0.41,0.88}{\vec{x}}-{\vec{x}}\hskip
1.0pt^{\prime}}$, characterizes the propagation of information from one point
to another.

Unlike for the BAO spectrum, we are not concerned that nonlinear evolution
will make a small change to the frequency or phase. Since nonlinear evolution
is not expected to create such a frequency out of nothing, the primary concern
is that nonlinear evolution will make a large change to the amplitude by some
incalculable amount, making it impossible to relate bounds on the nonlinear
spectrum to the amplitude in the initial conditions. To be dangerous, this
change must specifically alter the amplitude of a high-frequency oscillation
relative to the amplitude of the underlying smooth nonlinear matter power
spectrum. To proceed, we Fourier transform (2.9) to arrive at

$\delta^{\mathrm{w}}({\vec{k}},\tau)=\int\\!\frac{\mathrm{d}^{3}q}{(2\pi)^{3}}\,G({\color[rgb]{0.92,0.,0.}{\vec{k}}-{\vec{q}}},{\color[rgb]{0.25,0.41,0.88}{\vec{q}}\hskip
2.0pt};\tau)\,\delta_{\mathrm{in}}^{\mathrm{w}}({\vec{q}})\,.$ (2.10)

Note that it is the wavevector
${\color[rgb]{0.92,0.,0.}{\vec{k}}-{\vec{q}}}\equiv{\color[rgb]{0.92,0.,0.}{\vec{p}}}$
that characterizes the scale of the inhomogeneities. Following [11], we will
define a scale $k_{\omega}\ll k$ as the approximate period of oscillations in
the power spectrum (e.g. $k_{\omega}=2\pi/\omega_{\mathrm{lin}}$). We will
separate the domains $p>k_{\omega}$ and $p<k_{\omega}$ to distinguish the
“small-scale” and “large-scale” inhomogeneities, respectively.

Using an argument from [11], it is easy to see that small-scale
inhomogeneities do not contribute an oscillatory signal in the nonlinear
matter power spectrum.555We note that we would have
$G({\color[rgb]{0.92,0.,0.}{\vec{p}}},{\color[rgb]{0.25,0.41,0.88}{\vec{q}}\hskip
2.0pt};\tau)=\delta({\color[rgb]{0.92,0.,0.}{\vec{p}}}\hskip
1.0pt)\tilde{G}({\color[rgb]{0.25,0.41,0.88}q})$ if the matter distribution
was in the linear regime. As a result, the absence of an oscillatory signal
from $p>k_{\omega}$ in the nonlinear case is not a suppression of the
oscillation in the initial conditions, but the absence of a correction.
Suppose now that $p>k_{\omega}$ did contribute a high-frequency oscillatory
signal. This would imply that we should see a large change in
$\delta^{\mathrm{w}}(k)$ if we shifted
${\vec{k}}\to{\vec{k}}+{\vec{\alpha}}k_{\omega}$, where
$|\vec{\alpha}|\sim\mathcal{O}(1)$. However, recall that
$G({\color[rgb]{0.92,0.,0.}{\vec{p}}},{\color[rgb]{0.25,0.41,0.88}{\vec{q}}\hskip
2.0pt};\tau)$ is determined from the nonlinear density field in a universe
without the oscillatory signals and should therefore be a smooth function of
${\vec{p}}$ and ${\vec{q}}$. As a result, we can Taylor expand
$G({\color[rgb]{0.92,0.,0.}{\vec{k}}+{\vec{\alpha}}k_{\omega}-{\vec{q}}},{\color[rgb]{0.25,0.41,0.88}{\vec{q}}\hskip
2.0pt};\tau)$ in ${\vec{\alpha}}$ to find

$\displaystyle\delta^{\mathrm{w}}_{p>k_{\omega}}({\vec{k}}+\alpha{\vec{k}},\tau)$
$\displaystyle=\int\\!\frac{\mathrm{d}^{3}q}{(2\pi)^{3}}\,G({\color[rgb]{0.92,0.,0.}{\vec{k}}+{\vec{\alpha}}k_{\omega}-{\vec{q}}},{\color[rgb]{0.25,0.41,0.88}{\vec{q}}\hskip
2.0pt};\tau)\,\delta_{\mathrm{in}}^{\mathrm{w}}({\vec{q}})$
$\displaystyle\approx\int\\!\frac{\mathrm{d}^{3}q}{(2\pi)^{3}}\left[G({\color[rgb]{0.92,0.,0.}{\vec{k}}-{\vec{q}}},{\color[rgb]{0.25,0.41,0.88}{\vec{q}}\hskip
2.0pt};\tau)+k_{\omega}{\vec{\alpha}}\cdot\vec{\nabla}_{{\vec{k}}}G({\color[rgb]{0.92,0.,0.}{\vec{k}}-{\vec{q}}},{\color[rgb]{0.25,0.41,0.88}{\vec{q}}\hskip
2.0pt};\tau)\right]\delta_{\mathrm{in}}^{\mathrm{w}}({\vec{q}})$ (2.11)
$\displaystyle\approx\delta^{\mathrm{w}}_{p>k_{\omega}}({\vec{k}},\tau)+\mathcal{O}\left(\frac{\alpha
k_{\omega}}{k}\right),$

where we used $\nabla_{{\vec{k}}}\sim k^{-1}$ because $G$ is a smooth function
of ${\color[rgb]{0.92,0.,0.}{\vec{k}}-{\vec{q}}}$. From this argument, we see
that $\delta^{\mathrm{w}}_{p>k_{\omega}}({\vec{k}},\tau)$ is a smooth function
of $k$. This means that it does not change rapidly over one period of the
initial oscillations and, therefore, does not contribute to the oscillatory
signal. In other words, small-scale nonlinearities do not change the amplitude
of features in the power spectrum provided that they have a large enough
frequency.

To conclude the discussion of small-scale nonlinearities, we note that this
result is independent of the precise shape of the initial feature. It only
requires that the feature, $\delta^{\mathrm{w}}(k)$, changes by order one over
a very small range of ${\vec{k}}$. We do not require that it is sinusoidal or
that it is associated with a large physical scale in configuration space.
Furthermore, the result that $k_{\omega}/k\ll 1$ is not altered by small-scale
nonlinearities is the same condition which requires that the oscillation is
distinct from a smooth polynomial. The power spectrum on a scale $k$ is
“smooth” if $\frac{\partial\log P}{\partial\log k}\lesssim 1$. This implies
that our features will be “sharp” if

$\frac{\partial\log P^{\mathrm{w}}_{X}(k)}{\partial\log k}\gg 1\quad\to\quad
k\omega_{\mathrm{lin}}\gg 1\,,\quad\omega_{\mathrm{log}}\gg 1\,.$ (2.12)

These conditions are illustrated in the left panel of Fig. 1, where the linear
oscillations at small wavenumbers are smooth, while the logarithmic feature is
sharp on all scales. Of course, implicit in this discussion is that $k$ is a
scale where nonlinearities are important. In our universe, nonlinear effects
are strongly $k$-dependent and, therefore, primarily affect modes near the
nonlinear scale, $k\approx k_{\mathrm{nl}}$.

###  2.3 Damping from Large-Scale Nonlinearities

In the previous section, we established that short-wavelength inhomogeneities
cannot alter the appearance of an oscillatory signal in the initial
conditions. We now turn to the long-wavelength modes. In general, it is hard
to compute the consequences of nonlinearities on the matter power spectrum
from first principles. Nevertheless, it has been shown that the nonlinear
effects of large-scale modes on the BAO spectrum can be computed and resummed
in perturbation theory, resulting in a damping of the amplitude and shape of
the standard BAO signal (cf. e.g. [10, 53, 54, 55, 56, 57, 58, 59]). In the
following, we generalize this calculation to a generic linear feature and
further extend it to the case of logarithmically-spaced oscillations.

####  2.3.1 Perturbative Treatment

Our aim is to compute the damping of a generic oscillatory feature due to
long-wavelength modes. It is well known that a simple perturbative treatment
is not enough to capture the full damping effect in the case of the standard
oscillatory features in the matter power spectrum, the BAO signal [53, 54, 55,
56, 57, 58, 59]. As in that case, it is however also useful to start with a
perturbative treatment for generic features and then include nonperturbative
effects in the calculation.

The full one-loop power spectrum is given by

$P_{1\text{-loop}}(k)=\int\\!\frac{\mathrm{d}^{3}q}{(2\pi)^{3}}\left[6F_{3}({\vec{q}},-{\vec{q}},{\vec{k}}\hskip
1.0pt)P(k)+2F_{2}^{2}({\vec{q}},{\vec{k}}-{\vec{q}}\hskip
1.0pt)P(|{\vec{k}}-{\vec{q}}\hskip 1.0pt|)P(q)\right],$ (2.13)

where $F_{n}$ are the usual perturbation theory kernels (see [60] for a
review). For a generic oscillatory component $P^{\mathrm{w}}(k)$ of the matter
power spectrum, the effects of long modes $q$, with $q<\Lambda\equiv\epsilon
k$,666In this case, long modes $q$ are defined as those modes with wavelengths
much longer than the typical width $\sigma$ of the feature, $q<\Lambda\ll
2\pi/\sigma$. This implies
$\Lambda_{\mathrm{BAO}}\lesssim$0.6\text{\,}\mathrm{\mathit{h}}\text{\,}{\mathrm{Mpc}}^{-1}$$
for the standard BAO signal. However, in practice, we want to predict the
matter power spectrum at any wavenumber $k$, which implies that the
prescription for long modes should also satisfy $q\ll k$. Since the latter
requirement is stronger than the former for the range of scales under
consideration in this work, we choose the separation scale to be
$\Lambda=\epsilon k$, with $\epsilon\ll 1$ [55]. can be captured by

$P^{\mathrm{w}}_{1\text{-loop}}(k)=\frac{1}{2}\int^{\Lambda}\\!\frac{\mathrm{d}^{3}q}{(2\pi)^{3}}\,\frac{({\vec{q}}\cdot{\vec{k}}\hskip
1.0pt)^{2}}{q^{4}}\,P^{\mathrm{nw}}(q)\left[P^{\mathrm{w}}(|{\vec{k}}+{\vec{q}}\hskip
1.0pt|)+P^{\mathrm{w}}(|{\vec{k}}-{\vec{q}}\hskip
1.0pt|)-2P^{\mathrm{w}}(k)\right].$ (2.14)

Since features break the scale invariance of the matter power spectrum, we
cannot keep only the first few orders in the Taylor expansion
$P^{\mathrm{w}}(|{\vec{k}}+{\vec{q}}\hskip
1.0pt|)=P^{\mathrm{w}}(k)+{\vec{q}}\cdot\nabla_{{\vec{k}}}P^{\mathrm{w}}(k)+\mathcal{O}(q^{2}/k^{2})$,
but have to resum the entire series into the exponential

$P^{\mathrm{w}}(|{\vec{k}}+{\vec{q}}\hskip
1.0pt|)=\mathrm{e}^{{\vec{q}}\cdot\nabla_{{\vec{k}}}}P^{\mathrm{w}}(k)\,.$
(2.15)

We can therefore rewrite (2.14) as

$P^{\mathrm{w}}_{1\text{-loop}}(k)=\int^{\Lambda}\\!\frac{\mathrm{d}^{3}q}{(2\pi)^{3}}\,\frac{({\vec{q}}\cdot{\vec{k}}\hskip
1.0pt)^{2}}{q^{4}}\,P^{\mathrm{nw}}(q)\left[\cosh\hskip-1.0pt\left({\vec{q}}\cdot\nabla_{{\vec{k}}}\right)-1\right]P^{\mathrm{w}}(k)\,.$
(2.16)

We stress that we did not assume any particular form of $P^{\mathrm{w}}(k)$ in
this expression. It is only based on (i) calculating the contribution of long
modes to the one-loop power spectrum (2.13) and (ii) the matter power spectrum
having an oscillatory component of any kind. Since the calculation proceeds by
applying the operator
$\cosh\hskip-1.0pt\left({\vec{q}}\cdot\nabla_{{\vec{k}}}\right)$ to the wiggle
power spectrum $P^{\mathrm{w}}(k)$, we now consider the two oscillatory
feature models separately.

##### Linear features.

As a consequence of the baryon acoustic oscillations in the early universe,
there is an enhanced probability to find pairs of galaxies at a separation
given by the size of the sound horizon at the drag epoch,
$r_{s}\approx$150\text{\,}\mathrm{Mpc}$$. We therefore find a peak in the
galaxy two-point correlation function and linearly-spaced oscillations in
Fourier space, with the location and frequency given by $r_{s}$, respectively.
An enhanced probability of finding galaxy pairs at another distance scale
$\omega_{\mathrm{lin}}$ would produce the same signatures. This is why we can
compute the damping of such feature oscillations in exactly the same way as
for the BAO spectrum. We will review them here for a generic scale
$\omega_{\mathrm{lin}}\gtrsim$75\text{\,}\mathrm{Mpc}$$ in order to not be
affected by small-scale nonlinearities.

Since the BAO signal is itself a small contribution to the overall matter
power spectrum, we neglect its contribution and simply use
$P^{\mathrm{w}}(k)\equiv\leavevmode\nobreak\ P^{\mathrm{w}}_{\mathrm{lin}}(k)$
as defined in (2.8). Applying $2n$ gradients to
$P^{\mathrm{w}}_{\mathrm{lin}}(k)$ results in

$q^{i_{1}}\dots
q^{i_{2n}}\nabla_{k_{i_{1}}}\dots\nabla_{k_{i_{2n}}}P^{\mathrm{w}}_{\mathrm{lin}}(k)\approx(-1)^{n}q^{i_{1}}\dots
q^{i_{2n}}\hat{k}_{i_{1}}\dots\hat{k}_{i_{2n}}\,\omega_{\mathrm{lin}}^{2n}\,P^{\mathrm{w}}_{\mathrm{lin}}(k)\,,$
(2.17)

where we used the series expansion of the hyperbolic cosine and neglected
small corrections that arise from acting with the derivative operators on the
smooth envelope $P^{\mathrm{nw}}(k)$. Plugging this result into (2.16), we get

$P^{\mathrm{w}}_{1\text{-loop}}(k)=-k^{2}\,\Sigma^{2}_{\mathrm{lin}}\,P^{\mathrm{w}}_{\mathrm{lin}}(k)\,,$
(2.18)

where we defined

$\Sigma^{2}_{\mathrm{lin}}\equiv\frac{1}{6\pi^{2}}\int_{0}^{\Lambda}\\!\mathrm{d}q\,P^{\mathrm{nw}}(q)\left[1-j_{0}\left(q\,\omega_{\mathrm{lin}}\right)+2j_{2}\left(q\,\omega_{\mathrm{lin}}\right)\right],$
(2.19)

with the spherical Bessel function of the first kind $j_{n}(x)$. We note that
$\Sigma_{\mathrm{lin}}$ is independent of the wavenumber $k$ unless an
implicit dependence is introduced by taking $\Lambda=\Lambda(k)$ as mentioned
above. However, the damping scale $\Sigma_{\mathrm{lin}}$ of course depends on
redshift via the broadband power spectrum,
$P^{\mathrm{nw}}(q)=P^{\mathrm{nw}}(q,z)$.

##### Logarithmic features.

In contrast to linear features, logarithmically-spaced oscillations in Fourier
space do not have a simple intuitive interpretation in real space. However, we
can proceed with the calculation of the damping without additional caveats
because the derivation of (2.16) is valid for any oscillations in the matter
power spectrum. In Appendix A, we show how the operator
$\cosh\hskip-1.0pt\left({\vec{q}}\cdot\nabla_{{\vec{k}}}\right)$ acts on the
wiggle component of the linear matter power spectrum in the presence of
logarithmic features, which we introduced in (2.4) and (2.8). In consequence,
the one-loop wiggle spectrum can be computed to be

$P^{\mathrm{w}}_{1\text{-loop}}(k)=-k^{2}\left[\Sigma^{2}_{\mathrm{log}}(k)P^{\mathrm{w}}_{\mathrm{log}}(k)+\hat{\Sigma}^{2}_{\mathrm{log}}(k)\,\frac{P^{\mathrm{nw}}(k)}{\omega_{\mathrm{log}}}\frac{\mathrm{d}\delta
P_{\zeta}^{\mathrm{log}}(k)}{\mathrm{d}\log k}\right],$ (2.20)

where we introduced

$\displaystyle\Sigma^{2}_{\mathrm{log}}(k)$
$\displaystyle\equiv\frac{1}{4\pi^{2}}\int_{0}^{\Lambda}\\!\mathrm{d}q\,P^{\mathrm{nw}}(q)\int_{-1}^{1}\\!\mathrm{d}\mu\,\mu^{2}\left\\{1-\cos\\!\left[\omega_{\mathrm{log}}\log\\!\left(1-\frac{q\,\mu}{k}\right)\right]\right\\},$
(2.21) $\displaystyle\hat{\Sigma}^{2}_{\mathrm{log}}(k)$
$\displaystyle\equiv-\frac{1}{4\pi^{2}}\int_{0}^{\Lambda}\\!\mathrm{d}q\,P^{\mathrm{nw}}(q)\int_{-1}^{1}\\!\mathrm{d}\mu\,\mu^{2}\,\sin\\!\left[\omega_{\mathrm{log}}\log\\!\left(1-\frac{q\,\mu}{k}\right)\right],$
(2.22)

with $\mu=\hat{q}\cdot\hat{k}$. In contrast to the damping scales of the BAO
signal, $\Sigma_{\mathrm{BAO}}$, and of linear features,
$\Sigma_{\mathrm{lin}}$, which are constant in $k$, the damping factors
$\Sigma_{\mathrm{log}}$ and $\hat{\Sigma}_{\mathrm{log}}$ are scale dependent
(but similarly redshift dependent). Moreover, the one-loop wiggle power
spectrum is also no longer directly proportional to the oscillatory power
spectrum at linear order. We however note that (2.21) and (2.22) are only
valid if $q\ll k$. In this limit and for large enough values of
$\omega_{\mathrm{log}}$, these expressions can be simplified into a form
similar to (2.19) for linear features. Since we can always decompose a
logarithmic feature in a basis of linear oscillations, this also conforms with
our expectation to recover this result in the appropriate limit. For now, we
however choose to keep the calculation general and will discuss these limits
in detail below.

The crucial aspect of both one-loop results (2.18) and (2.20) is that they
correct the linear power spectrum by $\mathcal{O}(1)$-terms for a wide range
of parameter space and wavenumbers, exactly as in the case of the standard BAO
spectrum. In other words, we have
$P^{\mathrm{w}}_{1\text{-loop}}(k)\approx\mathcal{O}(1)\,P^{\mathrm{w}}_{\text{tree-
level}}(k)$ for $k\in$]0.1\text{\,}\mathrm{0}\,\mathrm{3}${}$ because
$k^{2}\,\Sigma^{2}\approx\mathcal{O}(1)$ in this range. This indicates that
the perturbative treatment is insufficient. Fortunately, it is possible to
compute the leading-order correction of long modes to the wiggle power
spectrum at all orders in perturbation theory.

####  2.3.2 Infrared Resummation

The infrared (IR) resummation of the large-scale bulk flows that damp the BAO
signal has been studied in various ways [53, 54, 55, 56, 57, 58, 59]. Here, we
follow the approach of [56], in which the class of loop diagrams that are most
IR-enhanced are first identified and then resummed into the nonperturbative
effect, the well-known exponential BAO damping. Their $L$-loop diagram is
given by

$P^{\mathrm{w}}_{L\text{-loop},\,\mathrm{LO}}(k)=\frac{1}{L!}\prod_{i=1}^{L}\left[\frac{1}{2}\int^{\Lambda}\\![\mathrm{d}q_{i}]\,P^{\mathrm{nw}}(q_{i})\mathcal{D}_{{\vec{q}}_{i}}\mathcal{D}_{-{\vec{q}}_{i}}\right]P^{\mathrm{w}}(k)\,,$
(2.23)

where the subscript ‘LO’ indicates that the leading-order IR-enhanced loops
are taken into account. Furthermore, we introduced the notation
$[\mathrm{d}q]=\mathrm{d}^{3}q$ and defined

$\mathcal{D}_{{\vec{q}}_{i}}P^{\mathrm{w}}(k)=\frac{{\vec{q}}_{i}\cdot{\vec{k}}}{q_{i}^{2}}\left(P^{\mathrm{w}}(|{\vec{k}}+{\vec{q}}\hskip
1.0pt|)-P^{\mathrm{w}}(k)\right)=\frac{{\vec{q}}_{i}\cdot{\vec{k}}}{q_{i}^{2}}\left.\left(\mathrm{e}^{{\vec{q}}_{i}\cdot\nabla_{{\vec{k}}^{\prime}}}-1\right)P^{\mathrm{w}}(k^{\prime})\right|_{k^{\prime}=k}.$
(2.24)

It is easy to verify that we exactly recover (2.16) for $L=1$. Since the rest
of the calculation depends on the form of $P^{\mathrm{w}}(k)$, we discuss the
linearly- and logarithmically-spaced oscillations again in turn.

##### Linear features.

By employing (2.18), it is straightforward to compute the wiggle power
spectrum of linear features at $L^{\mathrm{th}}$ order,

$P^{\mathrm{w}}_{L\text{-loop},\,\mathrm{LO}}(k)=\frac{\left(-k^{2}\Sigma_{\mathrm{lin}}^{2}\right)^{L}}{L!}\,P^{\mathrm{w}}_{\mathrm{lin}}(k)\,.$
(2.25)

Resumming these terms to all orders, we obtain

$P^{\mathrm{w}}_{\mathrm{IR},\,\mathrm{LO}}(k)=\sum_{L=0}^{\infty}P^{\mathrm{w}}_{L\text{-loop},\,\mathrm{LO}}(k)=\mathrm{e}^{-k^{2}\Sigma_{\mathrm{lin}}^{2}}P^{\mathrm{w}}_{\mathrm{lin}}(k)\,.$
(2.26)

This result is a generalization of the BAO expression with
$r_{s}\rightarrow\omega_{\mathrm{lin}}$ (cf. e.g. [56]). Since the value of
$\Sigma_{\mathrm{lin}}$ however saturates for
$\omega_{\mathrm{lin}}\gtrsim$75\text{\,}\mathrm{Mpc}$$, we can simply use the
BAO damping scale, $\Sigma_{\mathrm{lin}}\approx\Sigma_{\mathrm{BAO}}$. We can
therefore factor out the damping and write the full matter power spectrum as

$P_{m}(k)\approx
P^{\mathrm{nw}}_{m}(k)+\mathrm{e}^{-k^{2}\Sigma_{\mathrm{BAO}}^{2}}\left[P^{\mathrm{w}}_{\mathrm{BAO}}(k)+P^{\mathrm{w}}_{\mathrm{lin}}(k)\right],$
(2.27)

which constitutes a simple generalization of the known result for the standard
BAO signal.

##### Logarithmic features.

It is slightly less trivial to compute the expression of the
$L^{\mathrm{th}}$-order loop for logarithmic features with arbitrary frequency
$\omega_{\mathrm{log}}$. We proceed via induction by computing the first few
orders and then deriving the general formula. In this way, the IR-resummed
wiggle power spectrum is found to be

$P^{\mathrm{w}}_{\mathrm{IR},\,\mathrm{LO}}(k)=\mathrm{e}^{-k^{2}\Sigma^{2}_{\mathrm{log}}(k)}\cos\hskip-1.0pt\left(k^{2}\hat{\Sigma}^{2}_{\mathrm{log}}(k)\right)P^{\mathrm{w}}_{\mathrm{log}}(k)-\mathrm{e}^{-k^{2}\Sigma^{2}_{\mathrm{log}}(k)}\sin\hskip-1.0pt\left(k^{2}\hat{\Sigma}^{2}_{\mathrm{log}}(k)\right)\frac{P^{\mathrm{nw}}(k)}{\omega_{\mathrm{log}}}\frac{\mathrm{d}\delta
P_{\zeta}^{\mathrm{log}}(k)}{\mathrm{d}\log k}\,.$ (2.28)

We refer to Appendix A for further details on this calculation.

While this expression provides the leading nonperturbative damping for
logarithmic features, we find by explicit calculation that generically
$\hat{\Sigma}_{\mathrm{log}}\ll\Sigma_{\mathrm{log}}$. This can be understood
by noticing that the integrands in (2.21) and (2.22) can be expanded in
$\omega_{\mathrm{log}}(q\mu/k)$ and $\omega_{\mathrm{log}}(q\mu/k)^{2}$,
respectively. The fact that the integrals get their largest contributions from
$q\ll k$ explains the hierarchy between $\Sigma_{\mathrm{log}}$ and
$\hat{\Sigma}_{\mathrm{log}}$. As a result, when
$k^{2}\hat{\Sigma}_{\mathrm{log}}$ is large enough to be important, the signal
is already exponentially suppressed. It is therefore a good approximation to
set $\hat{\Sigma}_{\mathrm{log}}\approx 0$ for our analysis choice of
$\omega_{\mathrm{log}}\geq 10$ (see below), as confirmed by Fig. 11 in
Appendix B. We therefore get

$P^{\mathrm{w}}_{\mathrm{IR},\,\mathrm{LO}}(k)\approx\mathrm{e}^{-k^{2}\Sigma^{2}_{\mathrm{log}}(k)}P^{\mathrm{w}}_{\mathrm{log}}(k)\,.$
(2.29)

Moreover, it is straightforward to show that in the limit $q/k\ll 1$,
$\Sigma_{\mathrm{log}}(k)$ approaches the functional form of
$\Sigma_{\mathrm{lin}}$ with the substitution
$\omega_{\mathrm{lin}}\rightarrow\omega_{\mathrm{log}}/k$,777This result has
also been independently derived in [61] using $\omega_{\mathrm{log}}\gg 1$ as
an expansion parameter.

$\Sigma_{\mathrm{log}}(k)\approx\Sigma_{\mathrm{lin}}|_{\omega_{\mathrm{lin}}=\omega_{\mathrm{log}}/k}=\frac{1}{6\pi^{2}}\int_{0}^{\Lambda}\\!\mathrm{d}q\,P^{\mathrm{nw}}(q)\left[1-j_{0}\left(\frac{\omega_{\mathrm{log}}}{k}q\right)+2j_{2}\left(\frac{\omega_{\mathrm{log}}}{k}q\right)\right],$
(2.30)

with $\Sigma_{\mathrm{lin}}$ given by (2.19). Therefore, whenever
$\omega_{\mathrm{log}}/k\gtrsim$75\text{\,}\mathrm{Mpc}$$, we can use the same
approximation as in the linear case,
$\Sigma_{\mathrm{log}}(k)\approx\Sigma_{\mathrm{BAO}}$. As we also show in
Appendix B, this approximation is good enough for the scales and frequencies
of interest in an analysis of BOSS data.888These findings are also confirmed
by analyses of $N$-body simulations performed in [62]. We can therefore
further approximate (2.29) and write the full nonlinear power spectrum as

$P_{m}(k)\approx
P^{\mathrm{nw}}_{m}(k)+\mathrm{e}^{-k^{2}\Sigma_{\mathrm{BAO}}^{2}}\left[P^{\mathrm{w}}_{\mathrm{BAO}}(k)+P^{\mathrm{w}}_{\mathrm{log}}(k)\right],$
(2.31)

i.e. in the same way as for linear features in (2.27).

Before concluding the discussion of the theoretical damping calculation, a
couple of remarks are in order regarding subleading corrections to (2.27) and
(2.31):

  * •

We computed the leading-order IR-resummed power spectrum. It has however been
shown that there are subleading contributions which improve the fit to
$N$-body simulations in the case of the featureless BAO signal [62, 56, 57,
58, 59]. Since we do not employ the theoretically computed results, but fit
the damping scale in our data analysis below, we can neglect these corrections
for both the BAO signal (as in the standard BAO analyses), and the primordial
linear and logarithmic features.

  * •

We have not computed the damping of the mixed BAO-feature term of (2.7),
$P^{\mathrm{w}}_{\mathrm{BAO}}(k)\,\delta P_{\zeta}^{X}$. The size of this
contribution is of order $A_{\mathrm{BAO}}\times A_{X}\approx 0.05\times 0.01$
and, therefore, contributes less than per mil to the matter power spectrum.
For linearly-spaced oscillations, we checked that it is again a good
approximation for $\omega_{\mathrm{lin}}\gtrsim 2\hskip
1.0ptr_{s}\approx$300\text{\,}\mathrm{Mpc}$$ to also factor out the
exponential BAO damping for the mixed term. We expect this to also be the case
for logarithmic features. Consequently, we implement the mixed term in the
nonlinear matter power spectrum as follows:

$P_{m}(k)=P^{\mathrm{nw}}_{m}(k)+\mathrm{e}^{-k^{2}\Sigma_{\mathrm{BAO}}^{2}}\left[P^{\mathrm{w}}_{\mathrm{BAO}}(k)+P^{\mathrm{w}}_{X}(k)+P^{\mathrm{w}}_{\mathrm{BAO}}(k)\,\delta
P_{\zeta}^{X}(k)\right].$ (2.32)

Finally, from now on, we switch to the notation for the BAO damping scale that
has been adopted in data analyses,
$\Sigma^{2}_{\mathrm{BAO}}\to\Sigma_{\mathrm{nl}}^{2}/2$. With this, we are
ready to implement (2.32) and search for primordial feature models in
observational LSS data.

##  3 Feature Search in Large-Scale Structure

In this section, we introduce and establish our search for features in the
BOSS dataset. We propose an analysis in which the amplitude and frequency of
the linear and logarithmic features can be constrained. Moreover, we check the
validity of our approach on mock data and compare the results with the
expected constraining power obtained in forecasts.

###  3.1 BOSS DR12 Dataset and Analysis Pipeline

The approach of our analysis, which we introduce in the following, is very
general and we expect it to apply to a wide range of surveys. Having said
this, some of its aspects are particular to BOSS, such as the validity of some
of the employed approximations, and should therefore be revisited in future
analyses.

Our analysis is based on the BAO pipeline of [63] which utilizes nbodykit
[64]. We use the commonly employed density field reconstruction procedure [65]
to reduce the damping scale caused by gravitational evolution and move
information from higher-order statistics back to the power spectrum [66].999We
showed in §2.3 that features are damped by large-scale modes in the same way
as the BAO signal. In consequence, reconstruction will remove the damping of
features in the same way, as encoded by the post-reconstruction value of the
damping scale $\Sigma_{\mathrm{nl}}$. We then measure the galaxy power
spectrum following the steps described in [63, 67]. The corresponding
covariance matrix is obtained by measuring the power spectrum monopole in 999
mock catalogs (see §3.3 for more details on these mock catalogs). To extract
the BAO (and potential feature) signal, we marginalize over the smooth galaxy
broadband power spectrum,

$P^{\mathrm{nw}}_{g}(k)=B^{2}P^{\mathrm{nw}}(k)F(k,\Sigma_{s})+A(k)\,,$ (3.1)

with five polynomial terms

$A(k)=\frac{a_{1}}{k^{3}}+\frac{a_{2}}{k^{2}}+\frac{a_{3}}{k}+a_{4}+a_{5}k^{2}\,.$
(3.2)

Here, the bias parameter $B$ is used to marginalize over the power spectrum
amplitude, $P^{\mathrm{nw}}(k)$ is the linear no-wiggle power spectrum model
without any BAO signal and

$F(k,\Sigma_{s})=\frac{1}{(1+k^{2}\Sigma_{s}^{2}/2)^{2}}$ (3.3)

is the velocity damping term arising from the nonlinear velocity field.
Finally, the standard BAO signal and the oscillatory features left after the
marginalization described above are modeled as

$P_{g}(k)=P^{\mathrm{nw}}_{g}(k)\left\\{1+\left[O(k/\alpha)+\delta
P_{\zeta}^{X}(k)+O(k/\alpha)\,\delta
P_{\zeta}^{X}(k)\right]\mathrm{e}^{-k^{2}\Sigma^{2}_{\mathrm{nl}}/2}\right\\},$
(3.4)

where $O(k)\equiv P^{\mathrm{w}}_{\mathrm{BAO}}(k)/P^{\mathrm{nw}}(k)$ is the
standard linear BAO spectrum, $\alpha$ is the associated isotropic scaling
parameter and $\Sigma_{\mathrm{nl}}$ is the nonlinear damping scale, which we
keep as a free parameter. While we left the redshift dependence to be
implicit, this model for the galaxy power spectrum is appropriately evaluated
in each redshift bin. The fiducial $\Lambda$CDM cosmology is taken to be the
same as in [63], with matter density $\Omega_{m}=0.31$, physical baryon
density $\omega_{b}=0.022$, amplitude of linear matter fluctuations on
$8\text{\,}\mathrm{\mathit{h}}\text{\,}{\mathrm{Mpc}}^{-1}$ scales
$\sigma_{8}=0.824$, scalar spectral index $n_{\mathrm{s}}=0.96$ and Hubble
constant
$H_{0}=$67.6\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{Mpc}}^{-1}$$.
The linear and logarithmic features are contained in the relative primordial
spectrum $\delta P_{\zeta}^{X}(k)$, with $X=\mathrm{lin},\mathrm{log}$, that
was introduced in (2.3) and (2.4).

While the above model is described in terms of a continuous wavenumber $k$,
cosmological experiments can only access a finite number of modes due to their
finite survey volume. The associated fundamental mode $k_{f}$ is given by the
largest scale included in the dataset and is naturally setting a resolution
limit on the oscillation frequency that we can measure. In practice, this is
realized through the survey window function which introduces couplings between
modes separated by the fundamental mode or less [68, 69, 67, 70]. In addition,
we only measure the power spectrum in discrete bandpowers $P_{i}$, which
average wavenumbers $\vec{k}$ with $k\in[k_{i}-\Delta k/2,k_{i}+\Delta k/2)$,
where the bandwidth $\Delta k$ is a choice of the analysis. As a consequence,
a signal with an (effective) linear frequency above the Nyquist frequency,
$\omega_{\mathrm{Ny}}=\pi/\Delta k$,101010The fundamental Nyquist frequency of
a survey is determined by $k_{f}$, but for $\Delta k\gtrsim k_{f}$, it is the
bandwidth that sets the limiting scale in an analysis. will be aliased and is
therefore out of reach. Figure 4

Figure 4: Illustration of the impact of the finite volume observed in BOSS on
the imprint of primordial features in the linear wiggle spectrum for linearly-
(left) and logarithmically-spaced (right) oscillations with
$A^{\mathrm{sin}}_{X}=0.05$ and $A^{\mathrm{cos}}_{X}=0$. The frequency of the
standard featureless BAO spectrum (dark gray) is small enough to be
essentially unaffected by the effects of the bandpower estimation and the
window function. The larger the frequency of the features (light colors), the
larger the effect of bandpass filtering with $\Delta
k=$0.005\text{\,}\mathrm{\mathit{h}}\text{\,}{\mathrm{Mpc}}^{-1}$$ (medium
colors). Including an approximation of the window function (dark colors; see
Appendix B for details) further smooths the primordial wiggles. In particular
at the Nyquist frequency, which is given by
$\omega_{\mathrm{Ny}}=$929\text{\,}\mathrm{Mpc}$$ for the BOSS DR12 dataset
used in this work, the feature oscillations are completely aliased with only
the featureless BAO spectrum remaining. The light gray lines indicate the
estimated noise curves for the high-redshift bin of BOSS including the
nonlinear exponential damping which differs between the BOSS data (solid) and
mock catalogs (dashed). Finally, we note that the comparison of the left- and
right-hand panels together with these noise curves allows for another way to
estimate the reliable frequency range for the logarithmic features that is
complementary to Fig. 3.

highlights the effects that the finite bandwidth and the window function have
on the power spectrum, and illustrates why they limit the range of frequencies
that are accessible in an LSS analysis.

In the next sections, we apply this analysis pipeline to the Baryon
Oscillation Spectroscopic Survey (BOSS), which was part of SDSS-III [71, 72]
and used the SDSS multi-fibre spectrographs [73, 74] at the
$2.5\text{\,}\mathrm{m}$ Sloan Telescope [75] of the Apache Point Observatory
in New Mexico. We employ the final version of this dataset, known as data
release 12 (DR12) [76], which contains spectroscopic redshifts of
$1\,198\,006$ galaxies. The survey covered
$10\,252\text{\,}\mathrm{d}\mathrm{e}\mathrm{g}^{2}$ of the sky, divided in
two patches called the North Galactic Cap (NGC) and the South Galactic Cap
(SGC), and a redshift range of $0.2-0.75$. Following the main BOSS analysis
[77], we split this redshift range into two (independent) redshift bins given
by $0.2<z<0.5$ (‘low-$z$’) and $0.5<z<0.75$ (‘high-$z$’). While the standard
BOSS analysis uses $\Delta
k=$0.01\text{\,}\mathrm{\mathit{h}}\text{\,}{\mathrm{Mpc}}^{-1}$$ [63], we
employ a bandwidth of $\Delta
k=$0.005\text{\,}\mathrm{\mathit{h}}\text{\,}{\mathrm{Mpc}}^{-1}$$, which is
close to the fundamental mode of BOSS, to maximize the feature frequency range
accessible in this dataset. This limits our analysis to
$\omega_{\mathrm{lin}}\leq\omega_{\mathrm{Ny}}\approx$929\text{\,}\mathrm{Mpc}$$,
but we conservatively take
$\omega_{\mathrm{lin}}\leq$900\text{\,}\mathrm{Mpc}$$.

We analyze this dataset by producing a Markov Chain Monte Carlo (MCMC) with a
modified version of emcee [78] which includes the Gelman & Rubin convergence
criterion [79] with scale parameter $\epsilon<0.04$. Since the inflationary
signal under consideration is isotropic, we focus on the power spectrum
monopole and perform our analysis with one isotropic BAO parameter $\alpha$
per redshift bin. Moreover, we treat the low- and high-redshift bins of BOSS
DR12 independently. Since we use separate broadband marginalization parameters
for the NGC and the SGC, we fit for a total of 18 free parameters per redshift
bin:

$\alpha,\omega_{X},A^{\mathrm{cos}}_{X},A^{\mathrm{sin}}_{X};B^{\mathrm{NGC}},B^{\mathrm{SGC}},\Sigma_{s},\Sigma_{\mathrm{nl}},a^{\mathrm{NGC}}_{n},a^{\mathrm{SGC}}_{n}\,.$
(3.5)

We impose flat priors on all parameters, including the feature frequencies
which are sampled within $]100\text{\,}\mathrm{9}\mathrm{0}\mathrm{0}$ and
$[10,80]$ for linearly- and logarithmically-spaced oscillations, respectively.
These ranges are motivated by our discussion of small- and large-scale
nonlinearities in Section 2 (see also Fig. 3). Given that the primordial
feature parameters $\omega_{X}$, $A^{\mathrm{sin}}_{X}$ and
$A^{\mathrm{cos}}_{X}$ are independent of the redshift bin, we combine them
when inferring bounds on these models while marginalizing over the other
parameters (see Appendix C).

###  3.2 Forecasting Methodology for BOSS

To estimate the expected level of sensitivity, validate and cross-check the
described analysis of BOSS data, we perform two types of forecasts: based on
the Fisher information matrix and based on the likelihood itself. The Fisher
forecasts have to be used with care, but provide useful guidelines over a
large range of possible parameters and experimental configurations since they
are relatively fast to compute. On the other hand, the likelihood-based
forecasts present a more direct picture of the sensitivity and, in particular,
allow to study the effects of noisy data and injected signals. In addition, we
checked a number of theoretical approximations using Fisher forecasts (see
Appendix B).

Our general forecasting methodology and modeling is based on [48], with some
modifications that are detailed in Appendix B. Since we are only interested in
oscillatory features (and not broadband effects), our forecasts directly
employ the relative wiggle spectrum $\mathcal{O}_{\\!g}(k)\equiv
P^{\mathrm{w}}_{g}(k)/P^{\mathrm{nw}}_{g}(k)$ and not the galaxy power
spectrum $P_{g}(k)$ as the observable. After marginalizing over the bias $B$
and the polynomial coefficients $a_{n}$ of (3.2), we expect these two to be
identical. Working with $\mathcal{O}_{\\!g}(k)$ removes much of the degeneracy
with these broadband parameters and makes the forecasts more reliable (in
particular for the Fisher matrix). The Fisher information matrix $F_{ij}$ is
typically defined as the average curvature of the log-likelihood
$\log\\!\mathcal{L}(\vec{\theta})$ around a fiducial point in parameter space
spanned by $\vec{\theta}$. In our BAO forecasts, we will generally use
$\vec{\theta}\equiv\\{\alpha,A^{\mathrm{sin}}_{X},A^{\mathrm{cos}}_{X}\\}$,
with $X=\mathrm{lin},\mathrm{log}$, where $\alpha$ is the standard BAO
parameter, and $A^{\mathrm{sin}}_{X}$ and $A^{\mathrm{cos}}_{X}$ are the
respective feature amplitudes. As in [48], we employ a conservative broadband
marginalization scheme. Since the inverse Fisher matrix is the covariance
matrix for a Gaussian likelihood, the Cramér-Rao bound,
$\sigma(\theta_{i})\geq\sqrt{(F^{-1})_{ii}}$, provides a lower limit on the
marginalized constraints, with equality commonly assumed for Fisher forecasts.

The likelihood-based forecasts are based on the same modeling and have
previously been utilized successfully in [48, 49]. In this type of forecast,
we compute the likelihood function $\mathcal{L}(\vec{\theta})$ on a grid in
parameter space, given a specific fiducial (‘data’) spectrum computed for a
fixed set of parameters $\vec{\theta}_{\mathrm{fid}}$. When specified, the
fiducial model includes a random realization of the noise to mimic the scatter
of experimental data points due to shot noise and cosmic variance. In this
case, we talk about “noisy forecasts” which will be useful in our estimates of
the probability of experimental noise mimicking the presence of oscillatory
features. This is in contrast to the “noiseless forecasts”, for which the
experimental effects are only captured by the covariance matrix, as commonly
employed. (We emphasize that the latter forecasts are not cosmic variance
limited.) Except where noted otherwise, all of the following BOSS forecasts
are based on this likelihood-based approach.

Finally, we note that we produce two sets of forecasts as in [49]: one for the
comparison to the mock catalogs and another to compare to the results from the
actual BOSS data. This is due to the fact that the mock catalogs have a known
problem of overdamping the BAO spectrum which results in an approximately 30%
weaker signal for the traditional BAO analysis [63]. When comparing to mocks,
we use a (post-reconstruction) nonlinear damping scale of
$\Sigma_{\mathrm{nl}}\approx$7\text{\,}\mathrm{\mathit{h}}\text{\,}{\mathrm{Mpc}}^{-1}$$,
while we employ the standard (redshift-dependent) values otherwise (see [48]).
The gray noise curves in Fig. 4, which include the nonlinear damping terms,
indicate this difference and we can anticipate that the bounds on the feature
amplitudes will be stronger on the data than in the mocks.

###  3.3 Validation on Mock Catalogs and in Forecasts

To validate our analysis pipeline, we first perform an analysis on the
MultiDark Patchy mock catalogs [80], which mimic the galaxy clustering
behavior observed in BOSS. These mock data have been produced using
approximate gravity solvers and analytical-statistical biasing models. The
catalogs have been calibrated to an $N$-body-based reference sample extracted
from one of the BigMultiDark simulations [81], which was performed using
GADGET-2 [82] with $$3840$^{3}$ particles in a volume of
$($2.5\text{\,}{\mathrm{\mathit{h}}}^{-1}\text{\,}\mathrm{Gpc}$)^{3}$ assuming
a $\Lambda$CDM cosmology with $\Omega_{m}=0.307115$, $\Omega_{b}=0.048206$,
$\sigma_{8}=0.8288$, $n_{\mathrm{s}}=0.9611$ and
$H_{0}=$67.77\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{Mpc}}^{-1}$$.
The mock catalogs use halo abundance matching to reproduce the observed BOSS
two- and three-point clustering measurements [83]. This technique is applied
as a function of redshift to reproduce the BOSS DR12 redshift evolution. Note
that these are the same mock catalogs that we use to derive the covariance
matrix of our analysis as mentioned above. In addition, we remark that the
mock catalogs do not contain any inflationary features since they assume a
featureless primordial power spectrum (see Appendix B for a check with
injected signals).

We apply the described MCMC analysis pipeline to 100 NGC and SGC Patchy mock
catalogs for the low- and high-redshift bins. Since nearby feature frequencies
are correlated due to the finite range of wavenumbers used in the analysis,
$0.010.3\text{\,}\mathrm{\mathit{h}}\text{\,}{\mathrm{Mpc}}^{-1}$, we bin the
samples of the resulting Markov chains with widths of
$\Delta\omega_{\mathrm{lin}}=$10\text{\,}\mathrm{Mpc}$$ and
$\Delta\omega_{\mathrm{log}}=1.0$. These values have been obtained from the
Markov chains by estimating a scale-independent correlation length of the
feature frequencies $\omega_{X}$. For a number of different purposes, we
condense the Markov chains into the following three statistical quantities
shown in Fig. 5:

Figure 5: Comparison of the mock analysis with likelihood-based forecasts for
linear (top) and logarithmic features (bottom) in the high-redshift bin. We
present the standard deviation $\sigma(\theta_{i})$, mean $\bar{\theta}_{i}$
and significance $\bar{\theta}_{i}/\sigma(\theta_{i})$ for the BAO parameter
$\alpha$ and the amplitudes $A^{\mathrm{sin}}_{X}$ and $A^{\mathrm{cos}}_{X}$
as a function of the frequency $\omega_{X}$. We display these quantities for
the individual mock catalogs together with their mean, showing very good
agreement with the noisy and noiseless forecasts.

the mean value $\bar{\theta}_{i}$, the variance $\sigma^{2}(\theta_{i})$ and
their ratio, the significance $\bar{\theta}_{i}/\sigma(\theta_{i})$, for the
parameters
$\vec{\theta}=\\{\alpha,A^{\mathrm{sin}}_{X},A^{\mathrm{cos}}_{X}\\}$ as a
function of the frequency $\omega_{X}$. Apart from providing a validation of
important parts of our analysis pipeline, the comparison of these quantities
to those obtained in forecasts serves as a check of both the mock analysis
leading to results within the expected sensitivity and the forecasts being
suitable to compare to the data as well as to perform additional checks.
Moreover, the significance provides a metric that helps to quantify any
possible detections of features in the data analysis, for instance.

The results from the analysis of the low-redshift mocks are presented in Fig.
5, with the results from the high-redshift bin being similar. The middle
panels of the mean values clearly show that our pipeline results in an
unbiased estimation of both the BAO parameter $\alpha$ and the feature
amplitudes given that the mocks are generated from a featureless primordial
power spectrum. For linear features, the larger variance and non-zero mean
values around $\omega_{\mathrm{lin}}=$150\text{\,}\mathrm{Mpc}$$ indicate the
expected degeneracy between the primordial features and the standard BAO
spectrum with a sound horizon of $r_{s}\approx$150\text{\,}\mathrm{Mpc}$$.
(Note that the degeneracy is not perfect because the BAO spectrum is not a
perfect sine oscillation, but contains a $k$-dependent amplitude and phase
shift.) As expected, we also observe that the BAO parameter $\alpha$ is
independent of the primordial parameters away from the scale of the sound
horizon reproducing the constraints in the standard BOSS analysis [63]. The
fact that the constraints on the feature amplitudes become (slightly) weaker
with growing frequency can be attributed to the finite-survey effects of
bandpowers and window function discussed above (cf. Appendix B). Finally, the
right column shows that one typically finds a $2\sigma$ fluctuation in some
frequency bins for any given mock catalog, or for a given frequency bin for
some of the 100 mocks. This should not be surprising given the roughly 80
sampled (but partly correlated) frequencies. Furthermore, the fact that we do
not find many greater-than-$3\sigma$ fluctuations is consistent with the
statistical expectations.

While the fluctuations seen in individual mocks are consistent with the
variance inferred from the posterior, we would also like to know if this
variance is consistent with the expectations for this type of survey. For this
purpose, we turn to the likelihood-based forecasts whose results are shown in
direct comparison in Fig. 5. We see that the mean values of
$\sigma(\theta_{i})$ and $\bar{\theta}_{i}$ are in excellent agreement with
the noiseless forecasts across the entire mock catalogs. Furthermore, when a
specific realization of noise is added to these forecasts, one finds the
fluctuations in both the mean and the variance are consistent with the
fluctuations observed in the mocks. (As expected, the mean of the noisy
forecasts approaches the noiseless forecasts in the limit of many
realizations.) This therefore further establishes that the observed
fluctuations in the significances are entirely generated by and consistent
with the experimental noise of cosmic variance and shot noise of BOSS. In
other words, these fluctuations occur because we fit random noise. Together
with the extensive checks using forecasts presented in Appendix B, this
establishes our feature search in the clustering of galaxies and we can now
turn to the BOSS DR12 data.

##  4 First Large-Scale Structure Constraints

In this section, we discuss the constraints on primordial feature models that
we infer from the BOSS DR12 galaxy power spectrum. Furthermore, we compare and
combine these novel bounds with those obtained from current Planck CMB data.
We conclude this section with estimates of the future sensitivity of
cosmological observations of the CMB and LSS.

###  4.1 Limits on Features from BOSS DR12

We now apply our analysis and forecasting pipeline to the BOSS DR12 dataset.
Figure 6

Figure 6: Comparison of the analysis of the BOSS DR12 dataset with likelihood-
based forecasts for linear (top) and logarithmic (bottom) features in the low-
redshift bin with the same quantities as in Fig. 5. We again observe very good
agreement.

shows the posterior distributions derived from the Monte Carlo Markov chains
of the low-redshift bin in terms of the same characterizing statistical
quantities as in Fig. 5 for the mock catalogs. When comparing these results
inferred from the data chains with likelihood-based forecasts, we again find
very good agreement for the low-redshift bin and similar results for the high-
redshift bin. We reiterate that these forecasts differ from those in Fig. 5
especially in the value of the nonlinear damping scale. Since the smaller
damping scale in the data leads to a larger signal-to-noise ratio and
considerably extends the range of wavenumbers contributing to the feature
search, we observe a smaller variance, less scatter in the mean values and a
smaller (but statistically consistent) number of greater-than-$2\sigma$
fluctuations than in the mocks. The fact that the inferred significances in
the third column of Fig. 6 agree well with those found in the noisy forecasts
indicates that we do not have any significant detection of a feature, but
rather that the data analysis is consistent with fitting experimental noise.
We note that the oscillations in $\bar{A}_{\mathrm{log}}^{Y}$, $Y=\sin,\cos$,
that are visible towards smaller $\omega_{\mathrm{log}}$ in the noiseless
forecasts, arise due to interference of the logarithmic feature spectrum with
the BAO spectrum in the range
$k\sim$0.10.2\text{\,}\mathrm{\mathit{h}}\text{\,}{\mathrm{Mpc}}^{-1}$$. The
noisy forecasts show however that this does not impact our BOSS analysis.

Having established the reliability and robustness of our data analysis in the
amplitude parameterization of $(A^{\mathrm{sin}}_{X},A^{\mathrm{cos}}_{X})$,
we want to infer the constraints on the overall feature amplitude $A_{X}$
while marginalizing over the phase $\varphi_{X}$. Since we do not find any
significant detections (see also Appendix C), we are mainly interested in
deriving limits on the presence of primordial features which is why we take
$\varphi_{X}\in[0,2\pi)$ and the amplitude to be positive semi-definite:

$A_{X}=\sqrt{\left(A^{\mathrm{sin}}_{X}\right)^{2}+\left(A^{\mathrm{cos}}_{X}\right)^{2}}\,.$
(4.1)

In this way, we can directly infer the upper limits on $A_{X}$ at 95% c.l.
from the Markov chains of the low- and high-redshift bin, respectively. To
derive constraints from the entire BOSS DR12 data, we combine the two sets of
Markov chains by multiplying the binned posterior distributions. In this
process, we neglect a possible correlation between the BAO parameter $\alpha$
and the feature amplitudes $A^{\mathrm{sin}}_{X}$ and $A^{\mathrm{cos}}_{X}$.
As previously noted, this correlation is however small away from feature
frequencies around the BAO scale and taking this correlation into account
would only strengthen the inferred bounds. We refer to Appendix C for an
extensive discussion and further details.

We present the resulting first constraints on linear and logarithmic features
from large-scale structure data alone, i.e. without the inclusion of any other
external datasets or information, in Fig. 7.

Figure 7: Upper limits on the feature amplitude $A_{X}$ at 95% c.l. as a
function of the feature frequency $\omega_{X}$ for linear (left) and
logarithmic oscillations (right), $X=\mathrm{lin},\mathrm{log}$, from BOSS
DR12 data alone. The dotted and dashed lines show the bounds that are
separately inferred from the low- and high-redshift bins, while the solid line
indicates the current limits from LSS by combining the two BOSS redshift bins.

Our analysis limits the amplitude of these primordial feature models,
$A_{\mathrm{lin}}$ and $A_{\mathrm{log}}$, to be less than one to two percent
of the primordial scalar amplitude $A_{\mathrm{s}}$ in the range of feature
frequencies accessible with BOSS. Moreover, we do not find any significant
detections of features as expected from Fig. 6 (see also Appendix C).

###  4.2 Comparison with Planck CMB Bounds

While we present the first limits on feature models from LSS alone,
constraints have been inferred from CMB observations for more than a decade
(cf. e.g. [84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 1]). It is
therefore interesting to compare the deduced constraints. While the frequency
coverage is wider in the CMB, our LSS-only bounds interestingly improve the
limits derived from current CMB data by up to a factor of 2.3 and 3.1 for
$\omega_{\mathrm{lin}}\gtrsim$200\text{\,}\mathrm{Mpc}$$ and
$\omega_{\mathrm{log}}\gtrsim 20$, respectively. This is illustrated in Fig.
8,

Figure 8: Comparison of the 95% upper limits on the feature amplitudes
$A_{X}$, $X=\mathrm{lin},\mathrm{log}$, from LSS and the CMB for linear (left)
and logarithmic features (right). The solid lines indicate our new BOSS-only
results and are identical to the solid lines of Fig. 7. The bounds from Planck
2015 temperature (dotted) and temperature+polarization data (dashed) are for
the first time displayed as a function of feature frequency as well. Beyond
those frequencies which show a degeneracy with the standard BAO spectrum, the
BOSS data are able to improve over the CMB.

which directly compares the constraints on the feature amplitudes from our
BOSS analysis with those deduced from current CMB temperature (TT), and
temperature and polarization (TTTEEE) data released by the Planck
collaboration in 2015 [97, 98] (see Appendix D for details on these
limits).111111We show the constraints from both TT and TTTEEE since the Planck
collaboration had labeled the results employing high-multipole polarization
data as preliminary in 2015. Having said that, the available information on
feature models released by the collaboration seems to have remained fairly
stable between their 2015 and 2018 releases [96, 1]. We will therefore use the
polarization data when deriving the joint constraints below. Since the common
focus of previous analyses was on the best-fit points or the likelihood
improvement, we note that the limits on the feature amplitudes from the CMB
have not been shown as a function of frequency before.

The improvements of our LSS bounds over those from Planck are primarily the
consequence of two effects. First, the number of signal-dominated modes over
the employed range of wavenumbers in BOSS and Planck are roughly comparable
(approximately $k_{\mathrm{max}}^{3}V$ and $\ell_{\mathrm{max}}^{2}$,
respectively). Second, the imprint of high-frequency oscillations in the CMB
power spectra is suppressed relative to that in the matter power spectrum, as
shown in Fig. 19 of Appendix D. In combination, the signal-to-noise ratio of a
high-frequency feature is somewhat larger in BOSS than in Planck which leads
to a more stringent constraint.

Finally, we can infer the best current limits on primordial linearly- and
logarithmically-oscillating feature models by combining the BOSS and Planck
data.121212Combined analyses of CMB and LSS data have previously been explored
in [99, 100, 101, 102] by employing measurements of the linear matter power
spectrum over a limited range of wavenumbers without nonlinear modeling. These
joint constraints are derived in Appendix D and shown in Fig. 9.

Figure 9: Joint BOSS and Planck upper limits at 95% c.l. on the linear (left)
and logarithmic (right) feature amplitudes $A_{X}$,
$X=\mathrm{lin},\mathrm{log}$. The best current constraints come from a
combination of BOSS DR12 and Planck 2015 TTTEEE data (solid). We also show the
BOSS+Planck TT (dashed) results and include the BOSS-only bounds (dotted) for
comparison.

As expected, we observe that these bounds are dominated by and, therefore,
closely follow our limits from galaxy clustering data of BOSS except at
smaller frequencies. Generally speaking, the bounds on features in the
discussed range of frequencies $\omega_{X}$ are now established at the one-
percent level relative to the primordial power spectrum.

###  4.3 Future LSS and CMB Constraints

With the discussed improvements in the constraints on primordial features
inferred from BOSS over those derived from Planck CMB data, it is timely to
ask how these bounds will evolve with future CMB and LSS surveys. To this end,
we performed Fisher matrix forecasts for upcoming, planned and futuristic
experiments. We extend previous LSS forecasts (cf. e.g. [35, 36, 37, 38, 40,
42, 43, 44]) in a number of ways, in particular by taking the effects of
nonlinearities, bandpowers and window functions into account, and
(conservatively) marginalizing over further uncertainties in the broadband
power spectrum. Furthermore, we compare the reach of LSS surveys to that of
future CMB missions. In this section, we focus on linear features since most
other features can be easily decomposed into a basis of linear oscillations.
Before discussing the results of these forecasts, we briefly summarize our
approach and refer to Appendices B and D for further details.

For our future LSS forecasts, we use the relative wiggle spectrum
$\mathcal{O}_{\\!g}(k)\equiv P^{\mathrm{w}}_{g}(k)/P^{\mathrm{nw}}_{g}(k)$ as
the observable up to
$k_{\mathrm{max}}=$0.5\text{\,}\mathrm{\mathit{h}}\text{\,}{\mathrm{Mpc}}^{-1}$$
based on [48] as outlined in Appendix B, including the effects of
nonlinearities, bandpowers and window function.131313In this way, we also find
that our choice of
$k_{\mathrm{max}}=$0.3\text{\,}\mathrm{\mathit{h}}\text{\,}{\mathrm{Mpc}}^{-1}$$
in the described analysis captures essentially all the information on features
available in the BOSS DR12 dataset. To estimate the sensitivity of the CMB, we
directly follow the methodology of [48], employing perfectly delensed
temperature and polarization power spectra. The fiducial point in both cases
is a featureless $\Lambda$CDM cosmology consistent with the Planck
measurements [103, 104]. After computing the Fisher matrices in the amplitude
parameterization, we obtain the forecasted 95% upper limits on
$A_{\mathrm{lin}}$ by randomly sampling from the associated Gaussian
distributions and applying the same procedure as in our BOSS analysis (see
Appendices B and C).

The resulting forecasted sensitivity of several LSS and CMB experiments is
illustrated in Fig. 10.

Figure 10: Forecasted sensitivity for the “feature spectrometer” of linear
features. The potential reach of various LSS (solid) and CMB (dashed)
experiments to constrain the feature amplitude $A_{\mathrm{lin}}$ at a
confidence level of 95% is presented as a function of their frequency
$\omega_{\mathrm{lin}}$. We refer to the main text (and Appendix B) for the
details regarding the experiments and note that the positive-semi-definite
nature of $A_{\mathrm{lin}}$ is taken into account in the displayed estimates.
Large-scale structure surveys have the potential to improve over the CMB by
more than one order of magnitude, while the CMB will always dominate the reach
in feature frequency. As we discussed in Section 2, the LSS forecasts for
$\omega_{\mathrm{lin}}\lesssim$100\text{\,}\mathrm{Mpc}$$ should be treated
with care since these frequencies might be affected by the effects of
nonlinear gravitational evolution on small scales and be generally more
sensitive to the details of signal modeling.

Apart from BOSS and Planck, we included the planned surveys DESI [105], Euclid
[106] and CMB-S3 as an umbrella for the multiple upcoming CMB missions [107,
108, 109]. In addition, we show the potential reach of CMB-S4 [110] (or,
similarly, PICO [111]) and a ‘Future’ LSS experiment which is assumed to map
about ${10}^{8}$ objects up to redshift $z=3$ over half of the sky. To get a
sense for the theoretically possible limits, we also forecast a half-sky,
cosmic-variance-limited LSS survey with $z\leq 6$ and
$k_{\mathrm{max}}=$0.75\text{\,}\mathrm{\mathit{h}}\text{\,}{\mathrm{Mpc}}^{-1}$$
(‘LSS-CVL’),141414A similar performance could in principle be achieved by a
$21\text{\,}\mathrm{c}\mathrm{m}$ intensity mapping survey [112, 113, 114].
and a CMB experiment that measures the temperature and polarization spectra to
the cosmic variance limit up to multipoles of $\ell_{\mathrm{max}}^{T}=3000$
and $\ell_{\mathrm{max}}^{P}=5000$ on 75% of the sky (‘CMB-CVL’). We note that
the peak in the forecasted limits around
$\omega_{\mathrm{lin}}=r_{s}\approx$150\text{\,}\mathrm{Mpc}$$ is due to the
(approximate) degeneracy of the feature with the BAO spectrum over the signal-
dominated range of wavenumbers as previously discussed in §3.3.

The possibly most notable aspect of these forecasts is that the coming
generation of surveys, in particular DESI and Euclid, are projected to be more
sensitive than a cosmic-variance-limited CMB experiment over a substantial
range of frequencies. In this precise sense, large-scale structure will
permanently surpass the CMB in sensitivity. Equally significant is the
potential for future LSS observations to increase their constraining power on
feature models. The new BOSS limit presented in this work and the forecasts
for LSS-CVL leave approximately two orders of magnitude that could be
achievable with a suitably designed survey. As mentioned previously, it will
however be necessary to revisit some of the aspects of the analysis that we
employed on BOSS data to credibly achieve such sensitivities.

The improvements seen in future surveys come primarily from two factors:
smaller shot noise and higher redshifts. The constraining power of a survey is
dominated by the number of signal-dominated $k$-modes. Most of these modes are
at large wavenumbers, but are limited by the shot noise of the survey. The
significant increase in the number density of objects available in upcoming
surveys substantially increases the number of modes and drives the
improvements in sensitivity. In addition, the larger redshift range of these
observations means that the nonlinear damping is reduced, increasing the size
of the signal at higher wavenumbers.151515Recent advances in density field
reconstruction (see e.g. [115, 116, 117, 118, 119, 120, 121, 122]) also have
the potential to effectively further decrease the nonlinear damping scale. An
increase of the reconstruction efficiency from about 50% to 100% would lead to
an improvement in the forecasted sensitivity to $A_{\mathrm{lin}}$ of factors
in excess of 3, 2 and 1.5 for BOSS, DESI and LSS-CVL, for instance, and the
possibility of gaining signal-to-noise for larger values of
$k_{\mathrm{max}}$. This is a considerable improvement compared to the results
shown in Fig. 10 and substantially larger than what we expect for the BAO
frequency, i.e. the search for primordial features might give new motivation
to develop more efficient reconstruction techniques. Furthermore, future
surveys will also benefit from larger survey volumes which can be seen clearly
in the larger range of feature frequencies $\omega_{\mathrm{lin}}$ that are
accessible. This is because increasing the volume allows for finer $k$-bins,
which results in a larger Nyquist frequency $\omega_{\mathrm{Ny}}$.

In summary, LSS bounds on features are currently competitive with and will
surpass those from the CMB (present and future) over an increasing range of
frequencies. Large-scale structure observations have a significantly larger
sensitivity over their available frequency range due to the large number of
modes. Furthermore, the transfer of primordial power to the matter power
spectrum is more efficient than for the CMB, which leads to a larger intrinsic
signal (see Appendix D for a more detailed discussion). On the other hand, the
CMB can cover a wider range of frequencies than will be accessible even with
futuristic LSS surveys.

##  5 Conclusions

In this paper, we explored the impact of large-scale structure data on the
search for primordial features in the power spectrum. We showed theoretically
that such analyses are promising since they are not limited by the small-scale
nonlinearities of structure formation and the exponential damping caused by
large-scale bulk flows can be reliably computed (as we explicitly did at
leading order for both linear and logarithmic features). We then applied these
results to BOSS DR12 data and found constraints comparable to (but somewhat
stronger than) the best limits from Planck. The joint bounds on these models
are therefore dominated by the galaxy clustering data. Moreover, we forecast
that near-term surveys improve on this result by up to an order of magnitude
and could out-perform a cosmic-variance-limited CMB experiment over a
substantial range of feature frequencies.

Large-scale structure surveys offer great promise for dramatically improving
our understanding of the very early universe. However, to date, these hopes
have been largely limited by the modeling uncertainties around the nonlinear
scale. In this work, we have however shown that, for the right observables,
the statistical power of current surveys is already sufficient to
significantly impact our understanding of inflation and beyond.

While our emphasis was on primordial features, in particular from an
inflationary origin, both the method and the results have significantly
broader implications. Any sufficiently sharp feature in the matter power
spectrum could be analyzed in this way and could even be decomposed in a basis
of linear oscillations. We expect that constraints from LSS will be
competitive with those derived from the CMB, provided that the signal appears
directly in the (dark) matter and is not suppressed by the baryon fraction.

Finally, the statistical power of this approach is not limited to the power
spectrum and ultimately could be extended to higher-point statistics.
Primordial features are known to have associated non-Gaussian signatures (see
e.g. [123, 124, 125, 126, 127, 128]) which should similarly be robust to the
complications presented by nonlinear evolution. This presents the unique
opportunity not only to perform joint CMB power and bispectrum analyses [129,
95, 130, 1], but to also include the respective LSS observables. Furthermore,
the three-dimensionality of galaxy surveys may allow for entirely new types of
analyses that exploit the full angular dependence of higher-point correlation
functions. The universe has given us the unprecedented power of large-scale
structure to answer the most basic questions of our cosmic origins. We
contributed a small step towards this ultimate goal.

##### Acknowledgements

We thank Daniel Baumann, Dick Bond, Raphael Flauger, Mikhail Ivanov, Marcel
Schmittfull, Zvonimir Vlah and Matias Zaldarriaga for helpful conversations,
and Anagha Vasudevan, Mikhail Ivanov, Sergey Sibiryakov and Julien Lesgourgues
for sharing a draft of their work [61] with us. F. B. is a Royal Society
University Research Fellow. M. B. acknowledges support from the Netherlands
Organisation for Scientific Research (NWO), which is funded by the Dutch
Ministry of Education, Culture and Science (OCW), under VENI grant
016.Veni.192.210. D. G. is supported by the US Department of Energy under
grant no. DE-SC0019035. B. W. gratefully acknowledges support from the Marvin
L. Goldberger Membership at the Institute for Advanced Study, NSF Grant
PHY-1820775, a Research Studentship Award of the Cambridge Philosophical
Society, the Simons Foundation and a Visiting PhD Fellowship of the Delta-ITP
consortium, a program of NWO.

This work is based on observations obtained by the Sloan Digital Sky Survey
III (SDSS-III, http://www.sdss3.org/). Funding for SDSS-III has been provided
by the Alfred P. Sloan Foundation, the Participating Institutions, the
National Science Foundation and the US Department of Energy Office of Science.
This research also partly uses observations obtained with the Planck satellite
(http://www.esa.int/Planck), an ESA science mission with instruments and
contributions directly funded by ESA Member States, NASA and Canada. Parts of
this work were performed using resources provided by the Cambridge Service for
Data Driven Discovery (CSD3) operated by the University of Cambridge Research
Computing Service (http://www.csd3.cam.ac.uk/), provided by Dell EMC and
Intel using Tier-2 funding from the Engineering and Physical Sciences Research
Council (capital grant EP/P020259/1), and DiRAC funding from the Science and
Technology Facilities Council (http://www.dirac.ac.uk/). The DiRAC component
of CSD3 was funded by BEIS capital funding via STFC capital grants
ST/P002307/1 and ST/R002452/1, and STFC operations grant ST/R00689X/1. DiRAC
is part of the National e-Infrastructure.

We acknowledge the use of CAMB [131], CLASS [132], CosmoMC/GetDist [133],
IPython [134], MultiNest [135, 136] and the Python packages emcee [78],
Matplotlib [137], nbodykit [64] and NumPy/SciPy [138].

##  Appendix A Nonlinear Damping of Logarithmic Features

The effect of gravitational nonlinearities on the BAO signal has been
considered in various ways and can easily be extended to linear features as we
discuss in Section 2. The effects of large-scale gravitational nonlinearities
on logarithmically-spaced oscillations have however not been considered
previously.161616The authors of [61] independently performed this calculation.
In this appendix, we provide additional details on the computation of the
resulting damping of these features, complementing the discussion in §2.3. As
in the main text, we first detail the perturbative treatment and then resum
the infrared contributions to all orders in perturbation theory.

###  A.1 Perturbative Treatment

We have found in (2.16) that the effect of long modes on a generic wiggle
power spectrum $P^{\mathrm{w}}(k)$ at one-loop order implies the action of the
derivative operator
$\cosh\hskip-1.0pt\left({{\vec{q}}\cdot\nabla_{\vec{k}}}\right)$ on
$P^{\mathrm{w}}(k)$. For logarithmic features, we need to consider
$P^{\mathrm{w}}_{\mathrm{log}}(k)=P^{\mathrm{nw}}(k)\,A_{\mathrm{log}}\sin\\!\left[\omega_{\mathrm{log}}\log(\kappa)+\varphi_{\mathrm{log}}\right]$,
with $\kappa=k/k_{\star}$. As in the case of linearly-spaced oscillations
(2.17), we neglect small corrections that arise from applying the derivative
operators to the smooth envelope $P^{\mathrm{nw}}(k)$. Moreover, in order to
avoid clutter, we set $\varphi_{\mathrm{log}}=0$ in the following, but note
that it is straightforward to include the phase in the calculation. The
$(2n)^{\text{th}}$ derivative of the oscillatory part,
$\sin\\!\left[\omega_{\mathrm{log}}\log(\kappa)\right]$, is given by

$\displaystyle\nabla_{k_{i_{1}}}\dots\nabla_{k_{i_{2n}}}\sin\\!\left[\omega_{\mathrm{log}}\log(\kappa)\right]$
$\displaystyle=\frac{\hat{k}_{i_{1}}\cdots\hat{k}_{i_{2n}}}{2\mathrm{i}\,k^{2n}}\left[f_{n}(\omega_{\mathrm{log}})\,\kappa^{\mathrm{i}\omega_{\mathrm{log}}}-f_{n}(-\omega_{\mathrm{log}})\kappa^{-\mathrm{i}\omega_{\mathrm{log}}}\right]$
$\displaystyle=\frac{\hat{k}_{i_{1}}\cdots\hat{k}_{i_{2n}}}{k^{2n}}\left\\{\left[f_{n}(\omega_{\mathrm{log}})+f_{n}(-\omega_{\mathrm{log}})\right]\sin\\!\left[\omega_{\mathrm{log}}\log(\kappa)\right]\right.$
(A.1)
$\displaystyle\hphantom{=\frac{\hat{k}_{i_{1}}\cdots\hat{k}_{i_{2n}}}{k^{2n}}\\{\,}\left.-\mathrm{i}\left[f_{n}(\omega_{\mathrm{log}})-f_{n}(-\omega_{\mathrm{log}})\right]\cos\\!\left[\omega_{\mathrm{log}}\log(\kappa)\right]\right\\},$

where we employed

$\sin(x\log y)=\frac{1}{2\mathrm{i}}\left(\mathrm{e}^{\mathrm{i}x\log
y}-\mathrm{e}^{-\mathrm{i}x\log
y}\right)=\frac{1}{2\mathrm{i}}\left(x^{\mathrm{i}y}-x^{-\mathrm{i}y}\right),\quad
f_{n}(\omega_{\mathrm{log}})\equiv\frac{(\mathrm{i}\omega_{\mathrm{log}})!}{(\mathrm{i}\omega_{\mathrm{log}}-2n)!}\,.$
(A.2)

We can then perform the sum over $n$ and get

$\displaystyle\Bigr{[}\cosh\\!\left({{\vec{q}}\cdot\nabla_{\vec{k}}}\right)-1\Bigr{]}P^{\mathrm{w}}_{\mathrm{log}}(k)$
$\displaystyle=\left\\{\cos\left[\omega_{\mathrm{log}}\log\\!\left(1-\frac{q\,\mu}{k}\right)\right]-1\right\\}P^{\mathrm{w}}_{\mathrm{log}}(k)$
$\displaystyle\hphantom{=\
}+\left\\{\sin\left[\omega_{\mathrm{log}}\log\\!\left(1-\frac{q\,\mu}{k}\right)\right]\right\\}\left\\{A_{\mathrm{log}}\cos\\!\left[\omega_{\mathrm{log}}\log\\!\left(\frac{k}{k_{\star}}\right)\right]\right\\}P^{\mathrm{nw}}(k)$
$\displaystyle=\left\\{\cos\left[\omega_{\mathrm{log}}\log\\!\left(1-\frac{q\,\mu}{k}\right)\right]-1\right\\}P^{\mathrm{w}}_{\mathrm{log}}(k)$
$\displaystyle\hphantom{=\
}+\left\\{\sin\left[\omega_{\mathrm{log}}\log\\!\left(1-\frac{q\,\mu}{k}\right)\right]\right\\}\frac{P^{\mathrm{nw}}(k)}{\omega_{\mathrm{log}}}\frac{\mathrm{d}\delta
P_{\zeta}^{\mathrm{log}}(k)}{\mathrm{d}\log k}\,,$ (A.3)

where we defined $\mu\equiv\hat{k}\cdot\hat{q}$. It is also useful to consider

$\displaystyle\Bigr{[}\cosh\\!\left({{\vec{q}}\cdot\nabla_{\vec{k}}}\right)-1\Bigr{]}\frac{P^{\mathrm{nw}}(k)}{\omega_{\mathrm{log}}}\frac{\mathrm{d}\delta
P_{\zeta}^{\mathrm{log}}(k)}{\mathrm{d}\log k}$
$\displaystyle=\left\\{\cos\left[\omega_{\mathrm{log}}\log\\!\left(1-\frac{q\,\mu}{k}\right)\right]-1\right\\}\frac{P^{\mathrm{nw}}(k)}{\omega_{\mathrm{log}}}\frac{\mathrm{d}\delta
P_{\zeta}^{\mathrm{log}}(k)}{\mathrm{d}\log k}$ $\displaystyle\hphantom{=\
}-\left\\{\sin\left[\omega_{\mathrm{log}}\log\\!\left(1-\frac{q\,\mu}{k}\right)\right]\right\\}P^{\mathrm{w}}_{\mathrm{log}}(k)\,,$
(A.4)

since we need this expression in the following calculation of the IR-resummed
damping.

###  A.2 Infrared Resummation

In §2.3.1, we showed that
$P^{\mathrm{w}}_{1\text{-loop}}\approx\mathcal{O}(1)P^{\mathrm{w}}_{\text{tree-
level}}$, which suggests that all higher-order terms might be equally
important corrections to the linear wiggle power spectrum of logarithmic
features. In order to resum these infrared contributions, we need to evaluate
all the higher-loop diagrams of (2.23). In contrast to linear features, there
is no straightforward way to write down the $L$-loop contribution based on the
$1$-loop result, which is why we proceed by induction. Let us start with the
two-loop contribution, $L=2$, to get some intuition:

$\displaystyle P^{\mathrm{w}}_{2\text{-loop},\mathrm{LO}}(k)$
$\displaystyle=\frac{1}{8}\int^{\Lambda}\\!\frac{\mathrm{d}^{3}q_{1}}{(2\pi)^{3}}\frac{\mathrm{d}^{3}q_{2}}{(2\pi)^{3}}\,P^{\mathrm{nw}}(q_{1})P^{\mathrm{nw}}(q_{2})\mathcal{D}_{{\vec{q}}_{1}}\mathcal{D}_{-{\vec{q}}_{1}}\,\mathcal{D}_{{\vec{q}}_{2}}\mathcal{D}_{-{\vec{q}}_{2}}P^{\mathrm{w}}(k)$
$\displaystyle=-\frac{k^{2}}{4}\int^{\Lambda}\\!\frac{\mathrm{d}^{3}q_{1}}{(2\pi)^{3}}\,P^{\mathrm{nw}}(q_{1})\mathcal{D}_{{\vec{q}}_{1}}\mathcal{D}_{-{\vec{q}}_{1}}\left[\Sigma^{2}_{\mathrm{log}}(k)P^{\mathrm{w}}_{\mathrm{log}}(k)+\hat{\Sigma}^{2}_{\mathrm{log}}(k)\,\frac{P^{\mathrm{nw}}(k)}{\omega_{\mathrm{log}}}\frac{\mathrm{d}\delta
P_{\zeta}^{\mathrm{log}}(k)}{\mathrm{d}\log k}\right]$
$\displaystyle=\frac{k^{4}}{2}\left[\Sigma^{2}_{\mathrm{log}}(k)\left(\Sigma^{2}_{\mathrm{log}}(k)P^{\mathrm{w}}_{\mathrm{log}}(k)+\hat{\Sigma}^{2}_{\mathrm{log}}(k)\,\frac{P^{\mathrm{nw}}(k)}{\omega_{\mathrm{log}}}\frac{\mathrm{d}\delta
P_{\zeta}^{\mathrm{log}}(k)}{\mathrm{d}\log k}\right)\right.$
$\displaystyle\hphantom{=\frac{k^{4}}{2}\ \
\,}+\left.\hat{\Sigma}^{2}_{\mathrm{log}}(k)\left(\Sigma^{2}_{\mathrm{log}}(k)\frac{P^{\mathrm{nw}}(k)}{\omega_{\mathrm{log}}}\frac{\mathrm{d}\delta
P_{\zeta}^{\mathrm{log}}(k)}{\mathrm{d}\log
k}-\hat{\Sigma}^{2}_{\mathrm{log}}(k)P^{\mathrm{w}}_{\mathrm{log}}(k)\right)\right]$
$\displaystyle=\frac{k^{4}}{2}\left[\left(\Sigma^{4}_{\mathrm{log}}(k)-\hat{\Sigma}^{4}_{\mathrm{log}}(k)\right)P^{\mathrm{w}}_{\mathrm{log}}(k)+2\Sigma^{2}_{\mathrm{log}}(k)\,\hat{\Sigma}^{2}_{\mathrm{log}}(k)\,\frac{P^{\mathrm{nw}}(k)}{\omega_{\mathrm{log}}}\frac{\mathrm{d}\delta
P_{\zeta}^{\mathrm{log}}(k)}{\mathrm{d}\log k}\right],$ (A.5)

where we used (A.3) and (A.1). Here, we should note that the operator
$\mathcal{D}_{{\vec{q}}}\,\mathcal{D}_{-{\vec{q}}}$ does not act on
$\Sigma^{2}_{\mathrm{log}}(k)$ or $\hat{\Sigma}^{2}_{\mathrm{log}}(k)$ since
it only acts on the wiggle power spectra [56]. The three-loop term can be
derived along the same lines to be

$\displaystyle
P^{\mathrm{w}}_{3\text{-loop},\mathrm{LO}}(k)=-\frac{k^{6}}{6}\Biggr{[}$
$\displaystyle\left(\Sigma^{6}_{\mathrm{log}}(k)-3\Sigma^{2}_{\mathrm{log}}(k)\,\hat{\Sigma}^{4}_{\mathrm{log}}(k)\right)P^{\mathrm{w}}_{\mathrm{log}}(k)$
$\displaystyle+\left(3\Sigma^{4}_{\mathrm{log}}(k)\hat{\Sigma}^{2}_{\mathrm{log}}(k)-\hat{\Sigma}^{6}_{\mathrm{log}}(k)\right)\frac{P^{\mathrm{nw}}(k)}{\omega_{\mathrm{log}}}\frac{\mathrm{d}\delta
P_{\zeta}^{\mathrm{log}}(k)}{\mathrm{d}\log k}\Biggr{]}\,.$ (A.6)

Considering the one-, two- and three-loop contributions, it becomes apparent
that the structure at $L^{\text{th}}$ order is given by

$\displaystyle
P^{\mathrm{w}}_{L\text{-loop},\mathrm{LO}}(k)=\frac{(\mathrm{i}k)^{2L}}{L!}\Biggr{\\{}$
$\displaystyle\frac{1}{2}\hskip-1.0pt\left[\left(\Sigma^{2}_{\mathrm{log}}(k)+\mathrm{i}\hat{\Sigma}^{2}_{\mathrm{log}}(k)\right)^{\\!L}\\!\\!+\left(\Sigma^{2}_{\mathrm{log}}(k)-\mathrm{i}\hat{\Sigma}^{2}_{\mathrm{log}}(k)\right)^{\\!L}\right]\\!P^{\mathrm{w}}_{\mathrm{log}}(k)$
$\displaystyle\\!\\!+\frac{1}{2\mathrm{i}}\hskip-1.0pt\left[\hskip-1.0pt\left(\Sigma^{2}_{\mathrm{log}}(k)+\mathrm{i}\hat{\Sigma}^{2}_{\mathrm{log}}(k)\right)^{\\!L}\\!\\!-\left(\Sigma^{2}_{\mathrm{log}}(k)-\mathrm{i}\hat{\Sigma}^{2}_{\mathrm{log}}(k)\right)^{\\!L}\right]\\!\frac{P^{\mathrm{nw}}(k)}{\omega_{\mathrm{log}}}\frac{\mathrm{d}\delta
P_{\zeta}^{\mathrm{log}}(k)}{\mathrm{d}\log k}\Biggr{\\}}.$ (A.7)

The IR-resummed wiggle power spectrum of (2.28) is then obtained by resumming
all the loops.

##  Appendix B Large-Scale Structure Forecasts

We employ a suite of likelihood- and Fisher-based forecasts in particular to
validate and cross-check our analysis pipeline, and investigate the potential
reach of future surveys. In this appendix, we collect further details
regarding these LSS forecasts (§B.1) and collect the utilized experimental
specifications (§B.2). Furthermore, we provide additional checks of our
feature search (§B.3) as well as supplementary information for the forecasts
of future experiments (§B.4).

###  B.1 Forecasting with the Wiggle Spectrum

As previously stated, the wiggle spectrum is the main observable in our
forecasting pipeline, which was developed in [48] for the standard BAO
spectrum. In the following, we summarize its main aspects and introduce
further advances which especially include the use of bandpowers and the
convolution with a window function. These components are not required in a
wide range of applications, such as light relics, but are important to
reliably predict the sensitivity to (highly-)oscillating features.

We use two types of forecasts in this work, which are based either on the
Fisher information matrix $F_{ij}$ or on the likelihood function $\mathcal{L}$
itself. The former are computationally efficient and are therefore very useful
in particular to cover a large space of parameters and experimental
specifications. However, they only allow to access the standard deviation
around a fixed fiducial point assuming smooth noise and have to also be taken
with care given the involved approximations. We therefore only employ these
forecasts to estimate the sensitivity of future surveys and for a limited
number of tests. The likelihood-based forecasts come with a larger
computational cost, but are much more versatile. For instance, we can not only
obtain the standard deviations, but can also extract the mean values which
allows us to estimate significances and provides more direct comparisons with
MCMC analyses. In addition, it is possible to inject random noise realizations
and/or artificial feature signals. For these reasons, the majority of
forecasts in this work are of the latter type. In the following, we first
discuss the Fisher methodology, since it is commonly employed, and especially
highlight modifications to the standard approach. We then build on this
pipeline and introduce the likelihood-based forecasts.

####  B.1.1 Fisher Matrix Forecasts

Focusing on the oscillatory part of the power spectrum, the Fisher matrix of a
galaxy survey with multiple (independent) redshift bins $z$ can generally be
approximated by [48] 171717This is based on the standard Fisher matrix for
galaxy surveys of [139] which employs the galaxy power spectrum $P_{g}(k,\mu)$
as the observable.

$F_{lm}=\sum_{z,\,k_{i}}\frac{\Delta
k\,k_{i}^{2}}{(2\pi)^{2}}\int_{-1}^{1}\\!\frac{\mathrm{d}\mu}{2}\,\frac{D_{z}(k_{i},\mu)^{2}}{(1+D_{z}(k_{i},\mu)\,O_{z}(k_{i},\mu))^{2}}\frac{\partial\mathcal{O}_{z}(k_{i},\mu)}{\partial\theta_{l}}\frac{\partial\mathcal{O}_{z}(k_{i},\mu)}{\partial\theta_{m}}\,V_{\mathrm{eff}}(k_{i},\mu;z)\,,$
(B.1)

where $\mu$ is the cosine between the wavevector ${\vec{k}}$ and the line-of-
sight,
$\mathcal{O}_{z}(k,\mu)=P^{\mathrm{w}}_{z}(k,\mu)/P^{\mathrm{nw}}_{z}(k,\mu)$
is the (linear) relative anisotropic wiggle spectrum, $D_{z}(k,\mu)$ is the
nonlinear damping function and $V_{\mathrm{eff}}(k,\mu)$ is the effective
volume. In a featureless universe, the wiggle spectrum is simply the BAO
spectrum, while it may also contain a feature signal in our case. Since we
assume isotropic clustering,181818In other words, we take the limit of a
spherically-averaged clustering measurement. This is motivated by the fact
that the primordial information that we are interested in is strictly
isotropic and most of the information in BOSS is contained in the monopole
power spectrum. this quantity is given by

$\mathcal{O}_{z}(k,\mu)=\mathcal{O}_{z}(k)\equiv
B_{z}(k)\left\\{O(k/q;z)+\big{[}1+O(k/q;z)\big{]}\delta
P_{\zeta}(k)\right\\}+A_{z}(k)\,,$ (B.2)

with the linear BAO spectrum
$O(k;z)=P^{\mathrm{w}}_{\mathrm{BAO}}(k;z)/P^{\mathrm{nw}}(k;z)$ being
evaluated at the rescaled wavenumber
$k/q=D_{V}^{\mathrm{fid}}(z)/D_{V}(z)\,k$. This rescaling with the radial BAO
dilation $D_{V}\propto(D_{A}^{2}/H)^{1/3}$ is necessary because the
wavenumbers $k$ are derived from the measured angles and redshifts in a survey
using the angular diameter distance $D_{A}^{\mathrm{fid}}(z)$ and Hubble rate
$H^{\mathrm{fid}}(z)$ in a fiducial cosmology.191919In contrast to the data
analysis, we do not additionally rescale by the fiducial ratio of the sound
horizon in our forecasts since we recompute the BAO spectrum $O(k)$ for
different cosmologies using CLASS [132]. Moreover, we introduced the free
functions $B_{z}(k)$ and $A_{z}(k)$ which are taken to be smooth polynomials
in $k$ and distinct in each redshift bin, $\sum_{m}b_{m,z}k^{2m}$ and
$\sum_{m}a_{n,z}k^{n}$, with $m=0,\ldots,3$ and $n=0,\ldots,4$. By
marginalizing over these functions with fiducial values $a_{n,z}=0$,
$b_{0,z}=1$ and $b_{m\neq 0,z}=0$ in our forecasts, we effectively discard any
information in the observable that might be affected by nonlinearities,
biasing or observational systematics so that we only use a robust signal of
the primordial features and the standard BAO imprint. Finally, the nonlinear
damping and effective volume are implemented as

$D_{z}(k,\mu)\approx\mathrm{e}^{-k^{2}\Sigma_{\mathrm{nl}}(z)^{2}/2}\,,\qquad
V_{\mathrm{eff}}(k,\mu;z)\approx\left[\frac{\bar{n}_{g}(z)\,P_{g}(k,\mu;z)}{\bar{n}_{g}(z)\,P_{g}(k,\mu;z)+1}\right]^{2}V_{z}\,,$
(B.3)

where we assumed a constant nonlinear damping scale,
$\Sigma_{\mathrm{nl}}(k,\mu;z)\approx\Sigma_{\mathrm{nl}}(z)$, and position
independence of the comoving number density of galaxies,
$n_{g}({\vec{r}}\hskip 1.0pt)\approx\bar{n}_{g}=\mathrm{const}$, in each
redshift bin. Furthermore, the survey volume in a given redshift bin with
spherical geometry is denoted by $V_{z}$ and the fiducial galaxy power
spectrum by $P_{g}(k)$ which in particular includes the linear galaxy bias. We
note that we implicitly assumed in (B.1) that the feature spectrum is
nonlinearly damped in the same way as the BAO spectrum (cf. §2.3). (We
reiterate that this is a brief summary and all details can be found in [48],
including the modeling of the galaxy power spectrum, the nonlinear damping
scale and the effects of reconstruction.)

We have already written the Fisher matrix (B.1) as a sum over discrete
wavenumbers since the finite size of a galaxy survey introduces both a minimum
accessible wavenumber 202020The minimum wavenumber, or fundamental mode, which
is available in a survey with a spherical geometry is, in principle, given by
the survey volume $V$ according to $k_{\mathrm{min}}=2\pi[3V/(4\pi)]^{-1/3}$.
$k_{\mathrm{min}}$ and a minimum binning width in Fourier space given by the
fundamental mode, $\Delta k\geq k_{\mathrm{min}}$. For many current
applications, the width $\Delta k$ has become small enough so that the power
spectrum $P(k)$ is smooth in a given band $[k_{i}-\Delta k/2,k_{i}+\Delta
k/2]$ and we can approximate it by bandpowers $P_{i}\approx P(k_{i})$.
However, highly-oscillating primordial features introduce a significant
variation within any such band so that we have to compute the finite-size
bandpowers according to

$P_{i}=\frac{1}{\Delta k}\int_{k_{i}-\Delta k/2}^{k_{i}+\Delta
k/2}\\!\mathrm{d}k\,P(k)\,.$ (B.4)

The bandpower-averaged wiggle spectrum, which contains both the BAO and the
feature spectra, is then given by
$\mathcal{O}_{i}=\left(P_{i}-P^{\mathrm{nw}}_{i}\right)/P^{\mathrm{nw}}_{i}$,
for instance. To illustrate the effect of this averaging procedure (see also
Fig. 4), the bandpass-filtered primordial power spectrum (2.2) with linear
features is given by

$P_{\zeta,i}^{\mathrm{lin}}\approx
P_{\zeta,0}(k_{i})\left[1+\mathrm{sinc}(\omega_{\mathrm{lin}}\,\Delta
k/2)\,\delta P_{\zeta}^{\mathrm{lin}}(k_{i})\right],$ (B.5)

where $\mathrm{sinc}(x)=\sin(x)/x$ and we assumed $P_{\zeta,0}(k)\approx
P_{\zeta,0}(k_{i})$ for $k\in[k_{i}-\Delta k/2,k_{i}+\Delta k/2]$. This
implies that the oscillatory features are suppressed unless
$\omega_{\mathrm{lin}}\Delta k\ll 2$, or $\omega_{\mathrm{lin}}\ll 2/\Delta
k\approx$600\text{\,}\mathrm{Mpc}$$ for $\Delta
k=$0.005\text{\,}\mathrm{\mathit{h}}\text{\,}{\mathrm{Mpc}}^{-1}$$. For
logarithmic features, we could decompose the oscillations into linear features
in a given band and arrive at an analogous conclusion.

The second effect of a finite survey volume that we have to take into account
is the convolution of the power spectrum with the window function. This is of
course directly related to the bandpass filtering in reality although we
separate them here for convenience. Whereas the former averages the power
spectrum over the wavevectors ${\vec{k}}$ in a given band, the window function
introduces a coupling between otherwise independent wavenumbers. For an all-
sky survey with redshift range $[z_{-},z_{+}]$ and effective redshift
$\bar{z}$, the spherical top-hat window function is

$W_{\bar{z}}({\vec{x}})=W_{\bar{z}}(r,\theta,\phi)=\frac{1}{V_{\bar{z}}}\left[\Theta(d_{+}-r)-\Theta(d_{-}-r)\right],$
(B.6)

where $V_{\bar{z}}=4\pi(d_{+}^{3}-d_{-}^{3})/3$ is the bin volume,
$d_{\pm}\equiv d_{c}(z_{\pm})=\int_{0}^{z_{\pm}}\mathrm{d}z\,c/H(z)$ are the
comoving distances to the edges of the survey (or, equivalently, redshift bin)
and $\Theta(x)$ is the Heaviside step function. In practice, we however do not
have access to the full sky, but only to a fraction ${f_{\mathrm{sky}}}<1$.
For the purpose of our forecasts, we therefore include an incomplete sky by
restricting the integration over the azimuthal angle $\phi$:

$W_{\bar{z}}({\vec{x}})=\frac{1}{V_{\bar{z}}}\left[\Theta(d_{+}-r)-\Theta(d_{-}-r)\right]\Theta(2\pi{f_{\mathrm{sky}}}-\phi)\,.$
(B.7)

In this case, the Fourier transform of the window function is
radial,212121Although this is an idealized form of the window function, we
explicitly checked that forecasts employing the actual NGC and SGC window
functions of BOSS lead to consistent results.

$W_{\bar{z}}({\vec{k}})=W_{\bar{z}}(k)=\frac{3}{d_{+}^{3}-d_{-}^{3}}\left[\frac{d_{+}^{3}}{kd_{+}}j_{1}(kd_{+})-\frac{d_{-}^{3}}{kd_{-}}j_{1}(kd_{-})\right],$
(B.8)

with the spherical Bessel function of the first kind $j_{n}(x)$. By
restricting the power spectrum to finite-size bandpowers $P_{i}$ and using the
fact that the window function (B.8) is radial, we can rewrite the convolved
power spectrum, which is generally given by

$P^{c}(k,\bar{z})=\int\\!\frac{\mathrm{d}^{3}k^{\prime}}{(2\pi)^{3}}\,P(k^{\prime},\bar{z})\,W_{\bar{z}}^{2}({\vec{k}}-{\vec{k}}^{\prime})\,,$
(B.9)

in terms of a matrix equation:

$P_{i}^{c}=w_{ij}P_{j}\,,\qquad w_{ij}(\bar{z})=\frac{k_{j}^{2}\,\Delta
k}{(2\pi)^{2}}\int_{-1}^{1}\\!\mathrm{d}y\,W_{\bar{z}}^{2}\\!\left(\sqrt{k_{i}^{2}+k_{j}^{2}-2k_{i}k_{j}y}\right),$
(B.10)

which we can evaluate numerically for all $k_{i},k_{j}$ of an LSS survey. As
in the case of bandpowers, we again decompose the convolved spectrum
$P_{i}^{c}$ in its smooth and oscillatory components according to (2.6). To
summarize, the main extensions to the Fisher forecasting methodology of [48]
based on the wiggle spectrum are given in (B.4) and (B.10).

####  B.1.2 Likelihood-Based Forecasts

We also implemented forecasts based on the likelihood function
$\mathcal{L}(\vec{\theta})$ itself, as previously reported in [48, 49]. While
the modeling of the observables and covariances is the same as in the Fisher
analyses, we directly evaluate the likelihood function
$\mathcal{L}(\vec{\theta})$ on a grid in the parameter space 222222We note
that it is advantageous to employ the amplitude parameterization in terms of
$(A^{\mathrm{sin}}_{X},A^{\mathrm{cos}}_{X})$ over the phase parameterization
in both the Fisher- and likelihood-based forecasts for a few reasons. For
instance, we can use a fiducial featureless power spectrum, the likelihood in
$(A^{\mathrm{sin}}_{X},A^{\mathrm{cos}}_{X})$ is close to Gaussian and we do
not have to deal with the rather flat posterior in the phase $\varphi_{X}$. of
$\vec{\theta}=(\alpha_{z},A^{\mathrm{sin}}_{X},A^{\mathrm{cos}}_{X})$ as
follows:

$-2\log\mathcal{L}(\vec{\theta})=\chi^{2}(\vec{\theta})=\sum_{z,\,k_{i}}\Delta
k\,k_{i}^{2}\left[\mathcal{O}_{z}(k_{i};\vec{\theta})-\tilde{\mathcal{O}}_{z}(k_{i})\right]^{T}C^{-1}(k_{i},z)\left[\mathcal{O}_{z}(k_{i};\vec{\theta})-\tilde{\mathcal{O}}_{z}(k_{i})\right].$
(B.11)

Here, we used the theoretical (‘model’) wiggle spectrum
$\mathcal{O}_{z}(k_{i};\vec{\theta})$, the fiducial (‘data’) spectrum
$\tilde{\mathcal{O}}_{z}(k_{i})$ and the inverse covariance
$C^{-1}_{z}(k_{i})$ of the respective experiment. The latter is computed as in
the Fisher matrix (B.1) and includes the (white) instrumental noise
contribution, cosmic variance and the exponential nonlinear damping. We note
that all spectra are generally bandpass-filtered and convolved with the window
function as discussed above,
$\mathcal{O}_{i}^{c}=(P_{i}^{c}-P_{i}^{\mathrm{nw}\hskip-1.0pt,c})/P_{i}^{\mathrm{nw}\hskip-1.0pt,c}$,
which we have however omitted in (B.11) for ease of notation.

The model spectrum $\mathcal{O}_{z}$ varies over the considered parameter
space and is defined as

$\mathcal{O}_{z}(k;\vec{\theta};a_{i},b_{i})=B_{z}(k)\left\\{O_{\mathrm{fid}}(k/\alpha_{z},z)+\left[1+O_{\mathrm{fid}}(k/\alpha_{z},z)\right]\delta
P_{\zeta}(k)\right\\}+A_{z}(k)\,,$ (B.12)

where $O_{\mathrm{fid}}$ is the linear BAO spectrum of the fiducial cosmology,
$\alpha_{z}=\alpha(z)$ is the isotropic BAO parameter, and
$A_{z}(k)=\sum_{i=0}^{n}a_{i,z}k^{i}$ and
$B_{z}(k)=\sum_{j=0}^{m}b_{j,z}k^{2j}$ are the same polynomial ‘broadband’
polynomials as above, where six terms with $m=n=2$ turn out to be sufficient.
We marginalize over these terms by minimizing $\chi^{2}$ of (B.11) for these
parameters, i.e.
$\chi^{2}(\vec{\theta})=\min_{a_{n},b_{m}}\\!\big{\\{}\chi^{2}(\vec{\theta};a_{n},b_{m})\big{\\}}$.

The data spectrum $\tilde{\mathcal{O}}_{z}$ is computed by evaluating (B.12)
for a fiducial set of parameters $\vec{\theta}_{\mathrm{fid}}$ (which can
include non-zero feature amplitudes), with $B_{z}(k)=1$ and $A_{z}(k)=0$. In
addition to the smooth data with the experimental uncertainties being simply
captured by the covariance matrix, we also perform forecasts with ‘noisy
data’. In this case, we obtain the data spectrum by randomly picking the value
of $\tilde{\mathcal{O}}_{i}$ from a one-dimensional Gaussian distribution
function with mean $\tilde{\mathcal{O}}_{i}$ and variance $C(k_{i})$. This
therefore simulates the scatter of the actual measurement due to the expected
noise of an experiment (including sample variance) as captured by the
covariance matrix. We can include this in our forecasts in order to estimate
how likely it might be that features are found in the noise instead of the
data or, in other words, that the noise mimics the presence of oscillatory
features. In the main text, this constitutes an important check of the mock
and data analyses, and provides an estimate of the actual significance of
possible feature signals.

Having computed the likelihood function $\mathcal{L}(\vec{\theta})$ over all
of parameter space in which it is non-negligible, we then infer the predicted
posterior distribution $p(\theta_{l})$ of a parameter $\theta_{l}$ by
marginalizing over all other parameters $\theta_{m\neq l}$. Since the one-
dimensional posteriors for $\alpha$, $A^{\mathrm{sin}}_{X}$ and
$A^{\mathrm{cos}}_{X}$ are very close to Gaussian, we finally obtain the mean
$\bar{\theta}_{l}$ and standard deviation $\sigma(\theta_{l})$ through a
Gaussian fit to $p(\theta_{l})$.

###  B.2 Experimental Specifications

We do not only build on the signal modeling of [48], but also its
characterization of the LSS surveys (which was derived from [140]). In
general, we can characterize a cosmological galaxy survey by the following
quantities: redshift range, sky coverage, linear galaxy bias $b$ per redshift
bin and number (density) of objects $N_{g}$ ($\bar{n}_{g}$) in each redshift
bin. Here, we neglect the redshift error in spectroscopic surveys since it is
usually small compared to the damping scales, but would need to take it into
account for photometric observations. For planned experiments, such as DESI
and Euclid, we use specific values (see Appendix B of [48]), with Tables 1

$\bar{z}$ | $z_{\mathrm{min}}$ | $z_{\mathrm{max}}$ | $b$ | $${10}^{3}$\,\bar{n}_{g}\ [${\mathrm{\mathit{h}}}^{3}\text{\,}{\mathrm{Mpc}}^{-3}$]$ | $V\ [${\mathrm{\mathit{h}}}^{-3}\text{\,}{\mathrm{Gpc}}^{3}$]$ | $\Sigma_{\mathrm{nl}}\ [${\mathrm{\mathit{h}}}^{-1}\text{\,}\mathrm{Mpc}$]$
---|---|---|---|---|---|---
$0.350$ | $0.20$ | $0.50$ | $1.63$ | $0.275$ | $2.20$ | $4.6$
$0.625$ | $0.50$ | $0.75$ | $1.88$ | $0.142$ | $4.19$ | $4.4$
(a) Forecasts for BOSS DR12 data.

$\bar{z}$ | $z_{\mathrm{min}}$ | $z_{\mathrm{max}}$ | $b$ | $${10}^{3}$\,\bar{n}_{g}\ [${\mathrm{\mathit{h}}}^{3}\text{\,}{\mathrm{Mpc}}^{-3}$]$ | $V\ [${\mathrm{\mathit{h}}}^{-3}\text{\,}{\mathrm{Gpc}}^{3}$]$ | $\Sigma_{\mathrm{nl}}\ [${\mathrm{\mathit{h}}}^{-1}\text{\,}\mathrm{Mpc}$]$
---|---|---|---|---|---|---
$0.350$ | $0.20$ | $0.50$ | $2.04$ | $0.275$ | $2.20$ | $7.0$
$0.625$ | $0.50$ | $0.75$ | $2.34$ | $0.142$ | $4.19$ | $7.0$
(b) Forecasts for BOSS mock catalogs.

Table 1: Basic specifications for BOSS (inspired by [63] as detailed in [48])
with a sky area of
$\Omega=$10\,252\text{\,}\mathrm{d}\mathrm{e}\mathrm{g}^{2}$$ resulting in
roughly $1.2\text{\times}{10}^{6}$ objects in a volume of about
$6.4\text{\,}{\mathrm{\mathit{h}}}^{-3}\text{\,}{\mathrm{Gpc}}^{3}$. We
separately list the characteristic quantities employed when comparing to (a)
the BOSS DR12 data and (b) the corresponding mock catalogs since they differ
in the linear bias $b$ and the (post-reconstruction-equivalent) nonlinear
damping scale $\Sigma_{\mathrm{nl}}$ as discussed in the main text.

and 2

$\bar{z}$ | $0.05$ | $0.15$ | $0.25$ | $0.35$ | $0.45$ | $0.65$ | $0.75$ | $0.85$ | $0.95$
---|---|---|---|---|---|---|---|---|---
$b$ | $1.40$ | $1.48$ | $1.55$ | $1.61$ | $1.67$ | $2.05$ | $1.71$ | $1.71$ | $1.53$
$${10}^{3}$\,\bar{n}_{g}\ [${\mathrm{\mathit{h}}}^{3}\text{\,}{\mathrm{Mpc}}^{-3}$]$ | $38.8$ | $15.7$ | $3.96$ | $0.883$ | $0.0992$ | $0.591$ | $1.31$ | $0.920$ | $0.779$
$V\ [${\mathrm{\mathit{h}}}^{-3}\text{\,}{\mathrm{Gpc}}^{3}$]$ | $0.0357$ | $0.229$ | $0.563$ | $0.985$ | $1.45$ | $2.41$ | $2.86$ | $3.28$ | $3.66$

$\bar{z}$ | $1.05$ | $1.15$ | $1.25$ | $1.35$ | $1.45$ | $1.55$ | $1.65$ | $1.75$ | $1.85$
---|---|---|---|---|---|---|---|---|---
$b$ | $1.45$ | $1.48$ | $1.47$ | $1.47$ | $1.69$ | $1.68$ | $2.27$ | $2.45$ | $2.47$
$${10}^{3}$\,\bar{n}_{g}\ [${\mathrm{\mathit{h}}}^{3}\text{\,}{\mathrm{Mpc}}^{-3}$]$ | $0.466$ | $0.398$ | $0.387$ | $0.180$ | $0.133$ | $0.110$ | $0.0387$ | $0.0197$ | $0.0208$
$V\ [${\mathrm{\mathit{h}}}^{-3}\text{\,}{\mathrm{Gpc}}^{3}$]$ | $4.00$ | $4.30$ | $4.56$ | $4.79$ | $4.98$ | $5.14$ | $5.28$ | $5.39$ | $5.48$

Table 2: Basic specifications for DESI (derived from [105] as explained in
[48]), covering a sky area
$\Omega=$14\,000\text{\,}\mathrm{d}\mathrm{e}\mathrm{g}^{2}$$ and resulting in
roughly $2.7\text{\times}{10}^{7}$ objects in a volume of about
$59\text{\,}{\mathrm{\mathit{h}}}^{-3}\text{\,}{\mathrm{Gpc}}^{3}$.

updating the employed parameterizations of BOSS and DESI. For more distant
surveys, we assume a constant number density $\bar{n}_{g}$ for a given total
number of objects $N_{g}$ and a linear bias of $b(z=0)=1$. Our ‘Future’ LSS
survey contains $N_{g}=${10}^{8}$$ objects distributed over half the sky up to
redshift $z_{\mathrm{max}}=3$. The experiment referred to as ‘LSS-CVL’ is
cosmic variance limited on all employed scales and is designed to survey half
of the sky for $z\leq 6$. In our forecasts for BOSS, we generally take the
maximum wavenumber to be
$k_{\mathrm{max}}=$0.3\text{\,}\mathrm{\mathit{h}}\text{\,}{\mathrm{Mpc}}^{-1}$$
to coincide with the choice in the data analysis. All other (Fisher) forecasts
use
$k_{\mathrm{max}}=$0.5\text{\,}\mathrm{\mathit{h}}\text{\,}{\mathrm{Mpc}}^{-1}$$,
except for ‘LSS-CVL’ for which we choose
$k_{\mathrm{max}}=$0.75\text{\,}\mathrm{\mathit{h}}\text{\,}{\mathrm{Mpc}}^{-1}$$
since further extending the range of wavenumbers would likely yield only minor
improvements in sensitivity due to the exponential damping.

###  B.3 Additional Tests of the Pipeline

Given the described forecasting pipeline, we can provide additional insights
into our primordial feature search and discuss some of the tests that we
performed. In the following, we study the impact of the approximations in the
theoretical damping calculation on the BOSS constraints, revisit the impact of
the finite-volume effects and, in particular, test whether injected feature
signals can be detected in the analysis.

####  B.3.1 Check of Damping Assumptions

When computing the nonlinear damping of the linear and logarithmic
oscillations from large-scale bulk flows in §2.3, we made a number of
simplifying approximations which allowed us to use a single damping scale, the
standard BAO damping scale $\Sigma_{\mathrm{BAO}}$, in our data analysis. We
can explicitly check the validity of these approximations in Fisher forecasts
that generalize (B.1) to include the full resummed expressions for the linear
and logarithmic spectra of (2.26) and (2.28), and compare with the approximate
formulas of (2.27) and (2.31), respectively.

In order to perform this test, we need to numerically evaluate the three
damping scales of (2.19), (2.21) and (2.22), while choosing an appropriate
value of the cutoff scale $\Lambda$ which separates long modes $q$ from other
wavenumbers. The crucial point of the approximations is the fact that all the
computations are strictly valid in the regime of $q/k\ll 1$, i.e. a separation
of long and short modes. The cutoff $\Lambda$ therefore needs to be smaller
than the wavenumbers $k$ of interest. At the same time, however, all long
modes within the support of the feature also experience a damping effect. This
is the reason why it is sensible to take $\Lambda=\epsilon k$ for some
$\epsilon\ll 1$ (we employ $\epsilon=0.5$).232323We note that the logarithmic
damping factors $\Sigma_{\mathrm{log}}$ and $\hat{\Sigma}_{\mathrm{log}}$ are
not well defined in the limit $\Lambda\rightarrow k$ because the argument of
the logarithms in (2.21) and (2.22) approaches zero. This is precisely the
limit in which the computation becomes invalid since it is based on the
separation of long and short modes. Interestingly, this is not the case for
the BAO damping factor $\Sigma_{\mathrm{BAO}}$, whose value asymptotes for
$\Lambda\gtrsim$0.5\text{\,}\mathrm{\mathit{h}}\text{\,}{\mathrm{Mpc}}^{-1}$$
and can be integrated to $\Lambda\to+\infty$ without significantly affecting
the value of $\Sigma_{\mathrm{BAO}}$, even though the validity of the
nonlinear damping calculation breaks down at $\Lambda\sim k$ [55]. This choice
leads to all damping scales, including $\Sigma_{\mathrm{BAO}}$, to be
effectively $k$-dependent, $\Sigma_{X}\to\Sigma_{X}(k)$. Having said that, it
is important to remark once again that any dependence of these quantities on
the specific choice of the cutoff indicates that next-to-leading-order effects
should be taken into account (see e.g. [56] for the case of the standard BAO
signal). Since we fit $\Sigma_{\mathrm{BAO}}=\mathrm{const}$ in the data
analysis (as is standard), we also compute this damping scale for a
$k$-independent cutoff. Motivated by the maximum wavenumber of
$k_{\mathrm{max}}=$0.3\text{\,}\mathrm{\mathit{h}}\text{\,}{\mathrm{Mpc}}^{-1}$$,
we take
$\Lambda=$0.15\text{\,}\mathrm{\mathit{h}}\text{\,}{\mathrm{Mpc}}^{-1}$$ in
this case.

Figure 11

Figure 11: Impact of the various approximations to the theoretical damping
scales on the BOSS constraints for linear (top) and logarithmic features
(bottom). We display the relative difference of the Fisher-forecasted standard
deviation, $\delta\sigma=\sigma/\sigma_{\mathrm{full}}$, for the BAO parameter
$\alpha$ and the feature amplitudes $A^{\mathrm{sin}}_{X}$ and
$A^{\mathrm{cos}}_{X}$, where $\sigma_{\mathrm{full}}$ is obtained using the
full theoretical result. In the considered parameter space, the constraints
are essentially unaffected by $\hat{\Sigma}_{\mathrm{log}}(k)$. Here, we used
the effective post-reconstruction damping scales inferred at $z=0$.

shows the effect of the various approximations on the estimated constraints of
BOSS. We note that we evaluate the damping scales at redshift $z=0$ for
simplicity, given that the redshift dependence is the same for all damping
terms. This however also means that we effectively exaggerate the employed
damping scales and the actual impact on the constraints is even smaller than
shown. Even with this conservative choice, we can deduce that all of our
assumptions are valid in the context of the BOSS DR12 dataset. To be more
specific, assuming $\hat{\Sigma}_{\mathrm{log}}(k)\approx 0$ has basically no
visible impact on the constraints in the displayed parameter space of interest
in this work, as expected. In addition, approximating
$\Sigma_{\mathrm{lin}}(k),\,\Sigma_{\mathrm{log}}(k)\approx\Sigma_{\mathrm{BAO}}(k)$
only results in sub-percent variations to the constraints for
$\omega_{\mathrm{lin}}$ away from the BAO scale and
$\omega_{\mathrm{log}}\gtrsim 20$, and differences at the few-percent level
for $\omega_{\mathrm{log}}\in[10,20]$. Finally, taking $\Sigma_{\mathrm{BAO}}$
to be constant instead of computing it with a $k$-dependent cutoff penalizes
the constraints by roughly 3% for all linear and logarithmic frequencies. This
implies that all of the employed approximations are justified in the context
of the BOSS DR12 dataset and the upper limits that we infer in Section 4 are
in fact conservative. Nevertheless, the constraints inferred in future surveys
will likely benefit from using the theoretically-computed forms of the damping
scales $\Sigma_{\mathrm{lin}}$ and $\Sigma_{\mathrm{log}}$.

####  B.3.2 Impact of Finite-Volume Effects

Our ability to search for highly-oscillating features is limited by the fact
that we have only access to a finite cosmic volume, as we discussed in the
main text. Apart from introducing a cutoff at the Nyquist frequency due to
aliasing, the impact of finite-size bandpowers and the window function has to
be taken into account. We illustrate the consequences of these effects on the
sensitivity of BOSS in Fig. 12.

Figure 12: Illustration of the effects of the finite survey volume of BOSS on
the sensitivity to linear (top) and logarithmic features (bottom). We compare
the likelihood-based constraints on the BAO parameter $\alpha$, and the
feature amplitudes $A^{\mathrm{sin}}_{X}$ and $A^{\mathrm{cos}}_{X}$ when
using continuous spectra, $P(k)$, bandpass-filtered spectra, $P_{i}$, and
window-function-convolved spectra, $P_{i}^{c}$.

While the constraints on the BAO parameter $\alpha$ are essentially unchanged,
as expected given the BAO scale of $150\text{\,}\mathrm{Mpc}$, we observe a
gradual decrease in sensitivity to the feature amplitudes for larger
frequencies $\omega_{X}$. In consequence, we would overestimate the
constraining power of BOSS by up to a factor of two if we neglected the finite
volume of the survey.

These results can be easily understood in the context of Fig. 4 which shows
the impact of the finite-size effects on the spectra themselves. If we could
employ continuous spectra $P(k)$, a given primordial signal would have the
same amplitude independent of the feature frequency in the analysis, resulting
in the same sensitivity on all parameters (except for the interference with
the BAO signal). Since the amplitude effectively decreases for larger
$\omega_{X}$ when bandpass filtering the power spectrum [proportional to
$\mathrm{sinc}(\omega_{\mathrm{lin}}\,\Delta k/2)$ according to (B.5) for
linear features], the constraints gradually weaken and the feature model