# Back to the features: assessing the discriminating power of future CMB
missions on inflationary models

Matteo Braglia  Xingang Chen  Dhiraj Kumar Hazra  Lucas Pinol

###### Abstract

Future Cosmic Microwave Background (CMB) experiments will deliver extremely
accurate measurements of the E-modes pattern of the CMB polarization field.
Given the sharpness of the E-modes transfer functions, such surveys make for a
powerful detector of high-frequency signals from primordial features that may
be lurking in current data sets. With a handful of toy models that increase
the fit to the latest Planck data, but are of marginal statistical
significance, we use a state-of-the-art forecast pipeline to illustrate the
promising prospects to test primordial features in the next decade. Not only
will future experiments allow us to detect such features in data, but they
will also be able to discriminate between models and narrow down the physical
mechanism originating them with high statistical significance. On the other
hand, if the anomalies in the currently measured CMB spectra are just
statistical fluctuations, all the current feature best fit candidates will be
ruled out. Either way, our results show that primordial features are a clear
target of forthcoming CMB surveys beyond the detection of tensor modes.

###### Contents

  1. 1 Introduction
  2. 2 The models
  3. 3 Ruling out features if the Universe is featureless
  4. 4 Detecting features if  they are the true model of the Universe
  5. 5 Discriminating between feature models
  6. 6 Very similar features: the case for primordial non-Gaussianities
  7. 7 An estimator for correlation between feature models
  8. 8 Conclusions
  9. A Methodology
  10. B Computation of effective parameters
    1. B.1 Bump/dip in potential
    2. B.2 Turn model
  11. C Estimate of feature signal in the CMB angular spectra
    1. C.1 Resonant features: $\sin\log k$
    2. C.2 Sharp features: $\sin k$

##  1 Introduction

Precision cosmological observations set tight constraints on the shape of the
primordial power spectrum (PPS) from inflation. The most stringent ones come
from the latest release by the Planck satellite, which produced the most
accurate maps to date of $T$ and $E$ modes of the Cosmic Microwave Background
(CMB) anisotropies from largest to small scales [1]. Planck probes scales
$k\in[10^{-3},\,10^{-1}]\,{\rm Mpc}^{-1}$, corresponding to physical
wavelength being stretched above the Hubble radius, for example, $\sim 60-55$
$e$-folds before the end of inflation models with very high energy scales. At
those scales, the analysis of Planck data favors a parameterization of the
power spectrum of the curvature perturbation $\zeta$ with only 2 parameters,
$\mathcal{P}_{0}(k)\equiv\lvert\zeta(k)\rvert^{2}\,k^{3}/2\pi^{2}=A_{s}\left(k/k_{*}\right)^{n_{s}-1}$,
and accurately constrains them to $\ln(10^{10}\,A_{s})=3.045\pm 0.016$ and
$n_{s}=0.9649\pm 0.0044$ at 68% CL, using temperature and E-mode polarization
data [2]. Here $k_{*}$ is called the pivot scale, and is conventionally set to
$0.05\,{\rm Mpc}^{-1}$.

Simplest models of inflation are said to be of the Single-Field Slow-Roll
(SFSR) kind: a scalar field slowly rolls down a sufficiently flat potential,
thereby producing a spectrum of primordial fluctuations close to scale
invariance. The exact value of the deviation from scale invariance, encoded in
the parameter $n_{s}-1$, varies from model to model. Impressively, current
data is precise enough to discriminate models of SFSR inflation based on the
favored value for the spectral index $n_{s}$. In fact, another key prediction
of inflation is the production of a background of gravitational waves [3],
often parameterized by $r$ as the ratio of the amplitude of their spectrum to
that of the scalar perturbations, for which non-observations of primordial
B-modes from BICEP-Keck and Planck set the stringent upper limit $r_{0.05{\rm
Mpc}^{-1}}<0.036$ at 95% CL [4] (see also [5, 6, 7]). Joint constraints on
$(n_{s},r)$ favor models with concave potentials featuring a low $r$, such as
the original proposal of $R^{2}$-inflation [8], its $\alpha$-attractor
generalizations [9, 10], or small-field inflation models. The hope for a
detection of primordial gravitational waves have driven efforts towards
building next-generation experiments with greatly improved sensitivity to
polarization anisotropies of the CMB. In particular, primordial $B$ modes will
be probed at an unprecedented level of precision, with forecast sensitivities
on $r$ down to the level of $r\lesssim\mathcal{O}(10^{-3})$ [11]. Such
experiments include the Simons Observatory [12], LiteBIRD [13] and CMB-S4
[14]. Their respective (tentative) timelines are displayed in Fig. 1. While
several motivated models of inflation predict $r$ to be in the ballpark of
future sensitivities, our ignorance about the physics at the very high
energies does not allow us to make definitive predictions. In fact, there are
small-field realizations of Slow-Roll (SR) that predict $r$ to be even orders
of magnitude smaller than this. It is therefore of utmost importance to
consider observables beyond $(n_{s},r)$, in order to extract as much
information as possible from the high-quality data that future experiments
will deliver in the forthcoming years.

Primordial features provide an example of such observables [15, 16, 17, 18].
They are scale-dependent corrections induced by temporary or small deviations
from the SR evolution of the inflationary background. Such deviations can be
caused by features in the inflationary potential, in the internal field space
or in the couplings between the inflaton and other fields active during
inflation. Broadly, one can distinguish two families of features: sharp and
resonant features [19], depending on whether the deviation from SR is nearly
instantaneous or periodic. The resulting modulation of the (otherwise near-
scale-invariant) power spectrum shows up as an oscillatory pattern with a
characteristic running of the frequency, as well as a model-dependent
envelope. Specifically, sharp features generate a correction of the form
$\Delta\mathcal{P}_{\rm sharp}(k)\propto\mathcal{P}_{0}\sin(k/k_{0})$ where
$k_{0}$ is the mode that crosses the horizon at the time of the sharp feature,
while resonant features take the form $\Delta\mathcal{P}_{\rm
res}(k)\propto\mathcal{P}_{0}\sin(\omega\ln k+{\rm phase})$, where the running
is linear in $\ln k$ space. Some models of inflation, such as the so-called
Classical Primordial Standard Clocks (CPSC) predict a complicated
superposition of the two kinds of features [20, 21, 22, 23, 24, 25, 26, 27,
28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38].111Recently, it has also been
pointed out that an oscillatory feature in the inflaton potential, with a
bumpy modulation may address some interesting residuals in CMB data [39, 40,
41]. In these models, in general, only one of the two runnings shows up
clearly and the other acts as a modulation of the dominant running. Measuring
the parameters describing $\Delta\mathcal{P}$ would give access to the
fundamental physics at the origin of the signal, necessarily beyond the simple
SFSR paradigm.

Figure 1:  Timeline for future CMB experiments (left panel), their E-modes
noise curves (upper right panel) and associated error bars (lower right
panel). In the lower right panel, we also show the signal (blue line)
resulting from a resonant feature model with an amplitude
$\Delta\mathcal{P}/\mathcal{P}=0.05$, typically of the order of magnitude of
the features hinted at by Planck data. Note also that the error bars have been
centered around zero for illustration purposes only, so they can be easily
compared to the amplitude of the signal.

Interestingly, current data may already hint at primordial features.
Temperature data indeed show residuals that do not agree with the standard
cosmological model, with a statistical significance up to $3\sigma$, and that
can be explained by primordial features. Such anomalies include a systematic
suppression of power in the $TT$ angular spectra compared to the $\Lambda$CDM
best-fit for large scales, $\ell\lesssim 30$, and an oscillatory-like
structure of residuals at higher $\ell$, which is mostly significant around
$\ell\sim 750$. Many feature models have been compared to Planck data in the
literature (see [17, 18] for reviews). Although interesting best-fit
candidates have been identified, none of the feature models is statistically
favored over the SFSR paradigm, due to the introduction of penalizing extra
parameters.

Given the hints for them from the data, and the interesting physical
information that would be extracted by their detection, the purpose of this
paper is to consolidate the status of features as a target for future CMB
experiments. Indeed, although those experiments were mainly designed for the
detection of primordial $B$ modes, they will also achieve an unprecedented
precision in the measurement of $E$ modes. As shown in Fig. 1, this will
result in reducing the error bars on $EE$ power spectra (and $TE$ as well) to
the level of cosmic variance, which will be much more sensitive to features
than temperature spectra, due to their narrower transfer functions in $\ell$
space. Future CMB experiments therefore offer new prospects for assessing the
statistical significance of features by jointly constraining their correlated
effects on $T$ and $E$ modes, as also noticed in Ref. [42].

Specifically, the purpose of this paper is to address the two following
questions. Will we be able to rule out the above mentioned feature models if
they are truly only statistical fluctuations, and conversely detect them if
they are not? Supposing that a feature is detected in a future experiment, how
well can we distinguish between different feature models? In particular, the
second question has never been addressed, while it would crucially tell us to
which extent features can be used to constrain the physics of inflation. To
meet our goals, we use an efficient pipeline for the analysis of primordial
feature models, introduced in Ref. [33]. We consider several feature models
that improve the fit to current Planck data and we forecast the capability of
future experiments to tell them apart. Rather than using templates for the
feature signal, using BI-spectra and Non-Gaussianity Operator (BINGO) [43] and
its extension to two-field models [44], we numerically integrate the equations
for linear perturbations and arrive at exact power spectra, in order to fully
capture model-dependent details that may be overlooked using analytical
approximations. This enables us to explore fine details of the signal, that
such experiments will be sensitive to. In our analyses, we pay particular
attention to the implications of a future detection on early Universe model
building. Our results are very promising and generically prove that we will be
able to assess the statistical significance of primordial feature models
within the next decade of CMB observations.

Our paper is structured as follows. We start in Sec. 2 with a description of
the feature models that we use in our analysis, and present their
corresponding best-fit. In Sec. 3, we assume that hints for features
constitute statistical fluctuations in the current data and that the true
model of the universe is featureless; we forecast the corresponding
constraints on feature models that we would get. In Sec. 4, we answer the
opposite question, i.e. we assume for each of these models that it provides
the correct description of the Universe, and present the constraints on its
primordial power spectrum that we would get in that case. Sec. 5 is instead
devoted to the question of whether we can distinguish amongst different
feature models if a feature is detected in the data. For this, we assume CPSC
as a fiducial model for the Universe and discuss how much we could tell apart
this model from other feature models. Taking advantage of the properties of
CPSC, we also discuss the prospects for ruling out alternative scenarios to
inflation. Sec. 6 is dedicated to the particular case of feature models that
predict very similar angular power spectra: we argue that primordial non-
Gaussianities may be used as an additional source of information to
disentangle feature models. In Sec. 7, we introduce a simple estimator that
can be used to quickly quantify the difference between two feature signals and
roughly assess the prospects for distinguishing them in light of our results
without further data analyses. We discuss the implications of our results in
Sec. 8.

While our paper aims at being self-consistent, the discussion could be hard to
follow if we described every step of our analysis in the main text. Therefore,
to facilitate the reading, we collect all technicalities in Appendices. In
App. A, we describe our forecast pipeline and methodology. In App. B we
analytically estimate the feature signal for some models to motivate the
choice of the effective parameters used in our analysis. In App. C, we give a
rough analytical estimate of the imprints of primordial features on T and E
spectra to show their sensitivity to different runnings and frequencies of the
signals. In our paper, we use $\ln$ and $\log$ interchangeably to denote
natural logarithm and we explicitly write $\log_{10}$ to denote the logarithm
in basis 10.

##  2 The models

Figure 2:  Inflationary trajectories considered in our work. [Top-Left, CPSC
model] The inflationary trajectory eventually enters a curved path of valley
that excites the oscillations of a massive field (see insert). [Top-Right,
Resonant model] The trajectory is straight, but the inflaton potential
features oscillatory modulations. [Bottom-Left, Bump/Dip model] The straight
trajectory is shortly disturbed by a bump (or dip, not shown in the figure) in
the inflaton potential. [Bottom-Right, Turn model] Two straight inflationary
trajectories are connected by a sharp turn. These plots are for illustrative
purposes only, indeed the features consistent with real data correspond to
modifications in the trajectory that are not visible by eye, we have therefore
voluntarily omitted quantitative labels.

| CPSC | Resonant | Bump | Turn
---|---|---|---|---
$\ln B$ | $-1.2\pm 0.36$ | $-2.24\pm 0.38$ | $-1.44\pm 0.36$ | $-2.31\pm 0.36$
$\Delta\chi^{2}$ | 13.4 | 10.9 | 8.9 | 8.6

Figure 3:  [Top] Bestfit candidates to Planck data. [Bottom] Bayes factors
($\ln B$) and $\Delta\chi^{2}=\chi^{2}_{\rm featureless}-\chi^{2}_{\rm model}$
obtained for the analysis of Planck data.

The feature models we investigate in this work can be described by the
following Lagrangian.

$\displaystyle\mathcal{L}=$
$\displaystyle-\frac{1}{2}\left[1+\Xi(\phi)\sigma\right]^{2}(\partial\phi)^{2}-V_{\rm
inf}\left(1-\frac{1}{2}C_{\phi}\phi^{2}\right)\left(1+\frac{\Delta
V}{V}(\phi)\right)$
$\displaystyle-\frac{1}{2}(\partial\sigma)^{2}-\frac{1}{2}m_{\sigma}^{2}\sigma^{2}~{}.$
(2.1)

The field $\phi$ in the Lagrangian above represents the inflaton, and $V_{\rm
inf}\left(1-\frac{1}{2}C_{\phi}\phi^{2}\right)$ is a small-field slow-roll
potential leading to negligible tensor modes. The field $\sigma$ is a massive
field with $m_{\sigma}\gg H$, we assume it is initially stabilized at
$\sigma=0$. Features in the inflationary trajectory are introduced by non-zero
perturbations to either the slow-roll potential, $(\Delta V/V)$, or the
kinetic term $\Xi$ of the inflaton. In the former case, we can simply neglect
the dynamics of the massive field $\sigma$, as it does not move from the
minimum at $\sigma=0$, therefore making the inflationary dynamics effectively
of the single-field kind. On the other hand, if the coupling $\Xi$ is turned
on we have to follow the full two-field dynamics of the Lagrangian (2.1).

In the following, we will denote the model without features, i.e. that with
$\Xi=\Delta V/V=0$, as the featureless model. Before describing the feature
models considered in our paper, we find it useful for the reader to remind the
relations between the parameters of the slow-roll potential in Eq. (2.1) on
one hand, and the tilt and the amplitude of the power spectrum at a given
scale $k_{0}$222Note that in this paper we compute $A_{s}$ and $n_{s}$ at the
scale $k_{0}$ corresponding to the wavenumber that exits the horizon at the
time of the feature in the inflationary trajectory, that we will denote as the
feature scale, as opposed to the usual definition at the CMB pivot scale
$k_{*}=0.05\,{\rm Mpc}^{-1}$.  on the other hand:

$n_{s}-1=-2C_{\phi}\,\,\,\,\,\,\,\,{\rm
and}\,\,\,\,\,\,\,\,\ln\left(10^{10}\,A_{s}\right)=\ln\left(10^{10}\,24\pi^{2}\epsilon_{0}\,V_{\rm
inf}\right)$ (2.2)

where the first slow-roll parameter is approximately
$\epsilon_{0}=\frac{C_{\phi}^{2}\,\phi_{0}^{2}}{2}$ and background quantities
with a subscript $0$ are evaluated at the time implictly defined by the Hubble
crossing condition $k_{0}=a_{0}H_{0}$ for the scale of the feature. The
feature signals studied in this paper are all small scale-dependent
corrections to the leading order power-law power spectrum described by the
tilt and the amplitude in Eq. (2.2).

We will be considering 4 representative examples of feature models that we now
briefly describe. All the models introduce 3 extra parameters with respect to
the SFSR paradigm. This is the minimal number of parameter needed to describe
the amplitude, frequency and phase (or envelope, as in some models the phase
is uniquely fixed) of the feature. A schematic representation of the
inflationary trajectory for these models is shown in Fig. 2 and we plot the
corresponding best fit power spectra to Planck data in Fig. 3.

  * •

CPSC model: $\Xi=\xi\,{\rm\Theta}(\phi-\phi_{0}),\,\Delta V/V=0$, where
$\Theta$ is the Heaviside step function. This is a two-field model of
inflation, where the massive field $\sigma$ is excited by a sharp feature and
starts to oscillate, before eventually settling down to the minimum of its
(effective) potential [28, 29]. In this paper, we use the construction of Ref.
[34, 35], where the sharp feature is the introduction of a bending curved path
of valley, from the time $N_{0}$ at which $\phi(N_{0})=\phi_{0}$, in the
background evolution. The correction to the power spectrum consists in a sharp
feature signal at large scales, followed by a resonant signal at smaller ones.
However, unlike in the resonant model (see below) where the oscillations
appear at all wavenumbers, the resonant feature in this model is damped and
shows up only for a few wavenumbers This is easily explained as the background
oscillations of the damped field $\sigma$ only last for a finite range of
$e$-folds. The feature signal is described by 3 parameters: the scale of the
feature parameterized by $k_{0}$ which is the mode crossing the horizon at
time $N_{0}$, the frequency of the resonant signal $m_{\sigma}/H$ and the
maximum amplitude of the signal that can be written in terms of the model
parameters as [29, 35]:

$\left.\frac{\Delta\mathcal{P}_{\rm CPSC}}{\mathcal{P}_{0}}\right\rvert_{\rm
maximum}=\sqrt{4\pi}\,\xi^{2}\epsilon_{0}\left(\frac{m_{\sigma}}{H}\right)^{3/2}.$
(2.3)

The bestfit feature signal to Planck data has an amplitude
$\Delta\mathcal{P}/\mathcal{P}_{0}=0.023$, a frequency $m_{\sigma}/H=18.2$ and
the sharp feature exciting the massive field oscillations starts at
$k_{0}=2\times 10^{-3}\,{\rm Mpc}^{-1}$.

  * •

Resonant model: $\Xi=0,\,\Delta V/V=A\,\sin(\phi/\Lambda)$. We introduce a
periodic modulation to the slow-roll potential, introducing a typical physical
frequency $\omega>H$ related to $\Lambda$ and the slow-roll potential [19].
All physical wavenumbers $k/a$ eventually become comparable to the scale
$\omega$, therefore each mode experiences a brief sub-horizon resonance that
leads to a characteristic resonant feature correction to the PPS. The discrete
shift symmetry in the periodic oscillation of the potential determines the
running of the resonant signal to be of the form
$\sin\left(\frac{\omega}{H}\ln k\,{\rm Mpc}+{\rm phase}\right)$, with an
approximate constant amplitude. So, only for this simple model we simply use
the following template for the PPS in our analysis, rather than resorting to a
full numerical resolution of the inflationary dynamics from the potential:

$\frac{\Delta\mathcal{P}_{\rm
res}(k)}{\mathcal{P}_{0}}=\frac{\Delta\mathcal{P}}{\mathcal{P}_{0}}\sin\left(\frac{\omega}{H}\ln
k\,{\rm Mpc}+\varphi\right).$ (2.4)

The bestfit feature signal to Planck data has an amplitude
$\Delta\mathcal{P}/\mathcal{P}_{0}=0.015$, a frequency $\omega/H=18.1$ and a
phase $\varphi=5.9$. It is worth noting that the bestfit frequency of the
signal is the same as in the CPSC bestfit. In fact, the two spectra overlap in
the region $0.03\,{\rm Mpc}^{-1}\,\lesssim k\lesssim 0.07\,{\rm Mpc}^{-1}\,$
and fit similar residuals at the multipoles $400\lesssim\ell\lesssim 1000$,
providing a better improvement to the fit compared to the sharp feature models
that we are going to introduce. Taken at face value, this is a possible hint
at a mechanism active during inflation that produces a resonant signal with a
characteristic frequency of $\omega/H\sim 18$ at the scales of
$400\lesssim\ell\lesssim 1000$.

  * •

Bump or dip model: $\Xi=0,\,\Delta
V/V=A\,\exp\left[-\left(\phi-\phi_{0}\right)^{2}/\Delta^{2}\right]$. We
introduce a sharp feature with a Gaussian profile in the inflaton potential,
centered around the inflaton value $\phi_{0}$ and with a typical field-space
width of $\Delta$ [45]. Note that the profile of the feature is crucial as the
signal is very sensitive to it. For example, popular choices also include step
[46] or kink-like features [47]. These models have been studied extensively in
the literature, but they result in features that can fit well low-$\ell$
anomalies, but not high-$\ell$ ones, which are the focus of this paper. Models
that can address both large and small-scale anomalies do exist, but require
more than 3 extra parameters [48, 34, 35, 49]. As other sharp features, the
bump model results in oscillations linear in $k$-space. It can be described by
3 parameters: the mode $k_{0}$ crossing the horizon at the time of the sharp
feature $N_{0}$, the duration inflaton spends crossing the feature $\Delta
N\simeq 6\Delta/\sqrt{2\epsilon_{0}}$ and the maximum feature amplitude that
is found as (see App. B):

$\left.\frac{\Delta\mathcal{P}_{\rm bump}}{\mathcal{P}_{0}}\right\rvert_{\rm
maximum}=\frac{3A}{\epsilon_{0}}\sqrt{\frac{2\pi}{e}}.$ (2.5)

Note that the sign of the parameter $A$ in the potential is free (the feature
can be either a bump or a dip), so the effective parameter
$\Delta\mathcal{P}_{\rm bump}/\mathcal{P}_{0}$ can take both positive and
negative values, effectively amounting to a phase shift of $\pi$ in the
oscillations.

The bestfit feature signal to Planck data has an amplitude
$\Delta\mathcal{P}/P_{0}=0.04$, a frequency $k_{0}=4.1\times 10^{-3}\,{\rm
Mpc}^{-1}$ and a very short duration $\Delta N=0.09$ $e$-folds (note the
positive sign of the amplitude which means that the feature is due to a bump,
not a dip). Interestingly, this is the candidate with the largest amplitude.
However, due to the relatively large $k_{0}$, the amplitude of the feature
signal in $\mathcal{P}(k)$ is growing in the observable range and peaks at
very small scales to which Planck is not sensitive. Note, however, that since
the sensitivity of CMB angular spectra to sharp features decreases with $\ell$
(see Appendix C), this compensates with the growing amplitude and the CMB
residuals do not grow towards high multipoles. We note that, while the
multimodal nature of the (unconstrained) feature parameters is expected, the
reader might be surprised to see an apparent bimodality of the overall
amplitude $\ln 10^{10}\,A_{s}$, that presents two peaks at $\ln
10^{10}\,A_{s}^{A}\sim 3.135$ and $\ln 10^{10}\,A_{s}^{B}\sim 3.057$, the
latter being close to the bestfit value. In fact, this does not mean at all
that the amplitude is not constrained. Rather, it is a consequence of our
choice of computing $\ln 10^{10}\,A_{s}$ at the feature scale $k_{0}$, instead
of a fixed reference scale $k_{*}$, as done more often. While this is
convenient for all the other models, it leads to a slightly bimodal posterior
for $\ln 10^{10}\,A_{s}$ in the bump model, in proximity of two different
values of $N_{0}-N_{*}$. The bump/dip model attempts to fit the low-$\ell$ dip
around $\ell\sim 20-30$ in the Planck temperature data, which is also clear
from the bottom-left panel of Fig. 6 in Section 4. Therefore $k_{0}^{A}\ll
k_{0}^{B}$ and $N_{0}^{A}<N_{0}^{B}$, with the mode $B$ being the bestfit
candidate that fits high-$\ell$ residuals. To maintain the overall amplitude
of $\mathcal{P}(k)$ across a wide range of scales, we must have $\ln
10^{10}\,A_{s}^{A}>\ln 10^{10}\,A_{s}^{B}$, hence the apparent bi-modality.

  * •

Turn model:
$\Xi=\xi\exp\left[-\left(\phi-\phi_{0}\right)^{2}/\Delta^{2}\right],\,\Delta
V/V=0$. We introduce a localized kinetic coupling between the inflaton and the
massive field, which is stabilized at $\sigma=0$. The coupling has a Gaussian
profile and is switched on around the time $N_{0}$, where we used the same
notations as before. This model is phenomenologically similar, although not
exactly the same, as the turn model in Ref. [50]. We also note that the
difference between the turn model and CPSC model is that, in the former case,
the coupling is turned on only for a very short time leading to a purely sharp
feature, with no resonant signal imprinted, because the $\sigma$ field is
either not oscillating, or oscillating along a straight trajectory without
direct coupling to the inflaton so that the resonant signal is too small to
show up [51, 24]. In terms of field-space trajectory, while the CPSC is
described by a straight line followed by a bent trajectory, here the
trajectory simply consists in two straight lines joined by a sharp turn, as
can be seen from Fig. 2. As shown in App. B, the feature signal is very
similar to the dip one, except for the first oscillation at the largest scales
which is a bit more pronounced in this model. For this reason, it is
instructive to consider this model to investigate whether future experiments
will be sufficiently accurate such as to distinguish it from the dip one. The
turn model can be described by $k_{0}$, the duration of the feature $\Delta
N=6\Delta/\sqrt{2\epsilon_{0}}$ and the maximum amplitude (see App. B):

$\left.\frac{\Delta\mathcal{P}_{\rm turn}}{\mathcal{P}_{0}}\right\rvert_{\rm
maximum}=2\epsilon_{0}\left(\frac{\xi}{m_{\sigma}/H}\right)^{2}\sqrt{\frac{2\pi}{e}}.$
(2.6)

Note that, as opposed to the bump feature, here $\Delta\mathcal{P}_{\rm
turn}/\mathcal{P}_{0}$ can only take positive values, since it is
quadratically dependent on the coupling constant $\xi$. For this reason, the
model can only be similar to half of the feature signals produced by the
bump/dip model, i.e. those from a dip in the potential.

The bestfit feature signal to Planck data has an amplitude
$\Delta\mathcal{P}/P_{0}=0.02$, a frequency $k_{0}=2.5\times 10^{-3}\,{\rm
Mpc}^{-1}$ and duration $\Delta N=0.58$ $e$-folds, so the feature is not as
sharp as for the bump case.

We note that feature parameters are not constrained333For constraints on the
feature parameters from Planck data see triangle plots in the next Sections.,
consistently with the Bayes factors reported in the Table in Fig. 3, that show
no evidence for features. Nevertheless, Planck data do set stringent bounds on
the amplitude of such features. While depending quantitatively on the model at
hand, Planck bounds on the feature amplitude are generically
$\Delta\mathcal{P}/\mathcal{P}_{0}<\mathcal{O}(0.05)$ (see next Sections and
Ref. [2]). As can be seen from Fig. 3, the best-fit candidates for each model
address the oscillatory residuals in the high-$\ell$ temperature Planck data.
While the $\chi^{2}$ are different for each best fit, we see from the Bayes
factors that none of the models is statistically favored over the featureless
case. This can be explained by the fact that they require extra parameters,
which penalizes their Bayesian evidence. It is also interesting to see that
the correlated features in the $TE$ and $EE$ spectra are very different
depending on the model, but that we can only partially tap into such
information because the error bars are much larger than those in the $TT$
spectrum. This motivates the analysis of the next Sections, where we show that
shrinking such error bars will indeed open up a window on distinguishing
feature models one from the others, and assessing their statistical
significance.

As mentioned in the Introduction, our forecast analysis is based on picking
the best-fit candidate of these models, one model after the other, and
assuming it as a fiducial model of the universe. Then, other models are tested
against this fiducial universe by performing a nested sampling. In particular,
we will be using three combinations of present and future datasets, namely,
Planck + Simons Observatory (PL + SO), Planck + LiteBIRD + Simons Observatory
(PL + LB + SO) and Planck + LiteBIRD + CMB Stage-4 (PL + LB + S4)). We provide
more details on our forecast analysis in Appendix A. For an in-depth
discussion of all the technical details, as well as the noise specifications
of the experiments we use, we refer the reader to [35], which we follow
closely. If we wanted to exhaust all possible combinations of the 4 feature
models, plus the featureless one, we would have to perform 25 separate
analysis, which would unnecessarily complicate the understanding of the
physical results presented in this work. Therefore, we decided to perform only
some representative comparisons that highlight a specific information that we
will get access to with future CMB experiments. We will provide more details
in each of the next Sections.

##  3 Ruling out features if the Universe is featureless

Figure 4:  [Left] Constraints on $A_{s}$ and $n_{s}$ when the featureless
bestfit model is assumed as fiducial. We mark the Planck bestfit with a small
circle. In the legend we provide the values for the gain in the ${\rm FoM}$ in
going from one survey to another. [Right] Reconstructed signal from Planck
data (gray colormap) and from our forecast analysis for PL + LB + S4 (red
colormap). We plot different shades in going from 0 (inner regions, darker
shades) to 3$\sigma$ (outer regions, lighter shades). On top of the
reconstructed signal, we plot the current feature best fits to show that they
are consistent with Planck data, but also that they will be all excluded by
future measurements if the featureless fiducial is the true model of the
Universe.

Figure 5:  Constraints on the CPSC (top-left), resonant (top-right) and bump
(bottom-left) model. We assume the featureless bestfit model as fiducial. We
mark the Planck bestfit of each model with small circles. Projected Bayes
factors ($\ln B$) obtained for different experiments w.r.t. the Featureless
bestfit model, which is assumed as fiducial cosmology (bottom-right).

We start our analysis by addressing the most simple question, which is also
the most conservative one given the success of SFSR models. This question is
the following: would be able to rule out features in case they are just
statistical fluctuations? We therefore assume the featureless best fit as our
fiducial and compare our features model to it. We stress, however, that
presently feature models are only slightly disfavored with respect to the
featureless one, so, in the Bayesian interpretation, this hypothesis has no
special status. In order to see the power of future experiments to rule out
each of the feature best fit candidates, we perform the analysis against the
fiducial for all of the models but the Turn one, since it is very similar to
the Dip one and therefore we do not expect significant differences.

Let us first discuss the constraints on the featureless power spectrum that we
will get from future experiments. These are shown in Fig. 4. The improved
error bars on the E-mode measurements will significantly shrink the posterior
with an improvement of the FoM of a factor of $\mathcal{O}(10^{2})$ to
$\mathcal{O}(10^{3})$, with the better improvement which naturally comes by
adding both large scale and small scale measurements from LiteBird and S4. In
particular, the constraint on the spectral index $n_{s}=0.9661\pm 0.0019$ at
68 % CL from PL + LB + S4 will improve current constraints by roughly a factor
of two, setting under more pressure some of the popular single-field slow-roll
models444Note that the fiducial value of $n_{s}$ that we adopt in this Section
($n_{s}=0.9661$) is slightly different from the one quoted in the Introduction
($n_{s}=0.9649$). The reason is that the former is the fiducial from the
analysis with the EG20 likelihood [52] (adopted in this work), while the
latter comes from the analysis in the official Planck paper [2]. The two
values are consistent for current errors on $n_{s}$ and this does not affect
our analysis. .

In the right panel of Fig. 4, we also show the so-called predictive posterior
distribution on the power spectrum on which we overplot the feature best fit
candidates to Planck data555We use the code fgivenx [53], which was also used
for plotting the functional posterior distribution of feature models in [54]
and for free form reconstructions in [55]. . This clearly illustrates why
those candidates, although being perfectly consistent with current data, will
be completely excluded if the featureless model is really the one describing
our Universe. Indeed, all the candidates fall outside the red band
corresponding to the 1 to 3$\sigma$ constraint on the featureless
$\mathcal{P}(k)$.

We show in Fig. 5 the constraints on each feature model for each of the 3
experimental setups considered. All the contour plots share the same
characteristics. First of all, the parameter controlling the amplitude of the
feature is bounded more with respect to its current constraint. Second, all
the modes in the frequency or location parameters are washed out and the
corresponding forecast posteriors are almost flat, signaling that we will not
get any hints for features. From the Bayes factor analysis (see bottom-right
panel), the evidence against features will be generally strong, but not
decisive. This seemingly curious fact is instead rather easy to understand.
While feature models will be put under strong pressure, and we will rule out
all current best fits to Planck data, features with small amplitude will still
be consistent with future data. So although the Bayesian evidence will be
penalized by the number of extra parameters, the fit will always be slightly
better or equal to the one from the featureless model, never leading to a
decisive evidence against features.

##  4 Detecting features if  they are the true model of the Universe

Figure 6:  Reconstructed spectra for our feature models assuming in each case
that their best fits to P18 represent the true model of the Universe. Gray and
red-shaded regions correspond to 0 to 3$\sigma$ reconstructed CL from Planck
data and PL + LB + S4 forecasts respectively.

We now take a further step towards addressing the main question of this paper,
i.e. whether we will be able to distinguish different feature models using
future observations of the CMB polarization. In this short Section, we show
the evolution of the constraints on the power spectrum of our models going
from current data to forecasts for future experiments. This step will be very
helpful to interpret the results of the next Sections.

Technically, we run our nested sampling forecast for each feature model
assuming the corresponding Planck best fit, presented in Section 2, as
fiducial. Then, we plot the resulting predictive posterior distribution for
$\mathcal{P}(k)$. The constraints on the model parameters obtained from our
runs are shown in the next Sections. The results for each of our 4 feature
models are shown in Fig. 6, where, for simplicity, we only plot the predictive
posterior from our PL + LB + S4 analysis (blue, orange, green and red contours
for CPSC, resonant, bump and turn models respectively) on top of that from
real Planck data (gray contours)666Note that the reconstructed contours from
Planck in the case of the bump model show hints for the dip anomaly around
$\ell\sim 20-30$, which are not present in the case of the turn, which has a
different phase of the feature signal..

At this stage of our paper, the idea that we would like to convey is that the
constraints on the primordial features will become extremely tight in the
future, to the point that candidates from other feature models fall outside
the 3$\sigma$ posterior and thus will likely be ruled out. The shrinking of
the posterior is accompanied by an increase in the FoM which ranges from
$\mathcal{O}(10^{5})$ to $\mathcal{O}(10^{9})$ depending on the model (see
next Sections for the precise values), which reflects the compression of the
parameter space with respect to the analysis with real Planck data. We also
note that for each analysis, the featureless spectrum is ruled out with very
strong or decisive evidence depending on the model considered. We now move to
a more rigorous analysis, confirming all the statements we made in this
Section.

##  5 Discriminating between feature models

We now assume the CPSC best fit and compare to it the other feature models.
First of all, let us show in the top-left panel of Fig. 7 the projected
constraints on the CPSC model itself. Future experiments will be able to
identify, out of the various peaks of the multimodal posterior distribution
from Planck data, the main peak; the secondary ones will be washed out. The
constraints around the fiducial point will be extremely tight, with a
significant increase of the FoM up to $3.4\times 10^{5}$ for PL + LB + S4. In
particular, we will have a clear detection of the amplitude
$\Delta\mathcal{P}_{\rm CPSC}/\mathcal{P}_{0}$. We see that measurements at
high-$\ell$ from the SO will already provide a clear-cut hint to a non-zero
feature amplitude within the next decade. The addition of large-scale data
from LB will push the detection level beyond 4$\sigma$, leading to a
constraint $\Delta\mathcal{P}_{\rm CPSC}/\mathcal{P}_{0}=0.0241\pm 0.0057$,
i.e. a clear detection at 4.2$\sigma$, which will be further tightened by S4
to $\Delta\mathcal{P}_{\rm CPSC}/\mathcal{P}_{0}=0.0241\pm 0.0053$, reaching
the 4.5$\sigma$ detection [35]. Furthermore, once the amplitude
$\Delta\mathcal{P}_{\rm CPSC}/\mathcal{P}_{0}$ is detected, we will be able to
identify the frequency of the clock signal with extreme accuracy for
$m_{\sigma}/H=17.8^{+1.9}_{-1.5}$ for PL + SO, $m_{\sigma}/H=18.0\pm 1.4$ PL +
LB + SO and $m_{\sigma}/H=18.0\pm 1.1$ PL + LB + S4 respectively777Note that
these results are slightly different from those obtained in Ref. [35]. This is
because the models considered are not exactly the same. The one in Ref. [35]
was the particular limit of a more sophisticated model able to fit also large
scale anomalies. While the best fit of the model match to a very high
accuracy, the mapping of the effective amplitude parameter to its exact
numerical amplitude only worked as an order of magnitude estimate in Ref.
[35], while it is almost perfect here, leading to a small difference in the
constraints derived in the two papers. Furthermore, the location of the sharp
feature is parameterized in a different way in the two papers..

If CPSC can be distinguished from other feature models, these results will
amount to a clear detection of a massive particle with its mass precisely
measured relative to the Hubble rate and provide a direct evidence for the
phase of inflation responsible for the production of the anisotropies in the
CMB. However, while it is clear that the featureless form of the power
spectrum will be excluded if the CPSC best fit is really the model of the
Universe, it is less clear whether we could mistake other primordial feature
signals for the CPSC best fit. In fact, this is a general question that would
have to be asked in the event of a detection of any kind of features in the
density perturbations – to what extent are we able to distinguish different
feature models?

Figure 7:  Constraints on the CPSC(top-left), resonant (top-right), bump/dip
(bottom-left) and turn (bottom-right) model. We assume the CPSC bestfit model
as fiducial. We mark the Planck bestfit of each model with small circles. In
the bottom-most panel we plot the projected Bayes factors ($\ln B$) obtained
for different experiments w.r.t. the fiducial CPSC bestfit model.

To answer this question in our context, we thus compare our feature models to
the fiducial one (the CPSC current bestfit) and check whether the resulting
posterior distributions hint at any other model candidate and if the Bayesian
evidence shows preference for these models. We show the results of our
analysis in Fig. 7. Clearly, both the Bump and Turn models (bottom panels,
respectively left and right) show no constraints on the feature parameters,
except for a tighter constraint on the feature amplitude w.r.t. the one from
Planck data. On the other hand, it is interesting to see that the Resonant
model (upper right panel) does show a peak in the posterior distribution
corresponding to the frequency of the clock signal in the fiducial model. The
peak in the posterior for the amplitude however is consistent with $0$ at
2$\sigma$, if we approximate the slightly skewed posterior on
$\Delta\mathcal{P}_{\rm res}/\mathcal{P}$ to a Gaussian. This can be explained
as, although the model decently fits the first 2 oscillations of the clock
signal, its signal has a constant amplitude which does not fit well the
resonant part of the CPSC fiducial, whose oscillations are instead damped.
Furthermore, the Resonant model is not able to fit the sharp feature part of
the CPSC best fit. The plot in the bottom-most panel of Fig. 7 shows that
there would be mounting evidence in favor of the CPSC best fit, consistent
with what we have just discussed. We thus conclude that we will be able to
rule out sharp features as the origin of the feature signal and put under very
strong pressure a purely resonant signal, if the true model of the Universe is
the CPSC best fit to the current data.

In Sec. 7 we will show that, globally over all scales, models used in the
above analyses are all $40\%$ to $80\%$ correlated with the fiducial CPSC
model. Despite these correlations, we see that we will still be able to
distinguish them using CMB alone due to the rich information feature signals
carry in their phases and envelops across different scales. So our analysis
shows the considerable discriminating power of upcoming CMB experiments for
pinpointing the physical mechanism during inflation responsible for the
feature signal. We expect such a prospect to become even more promising when
information of other tracers of density perturbations, such as galaxies,
cosmic hydrogen atoms and stochastic gravitational waves are added.

Figure 8:  Constraints on the general clock template from a forecast for
PL+LB+S4 assuming the CPSC best fit.

Ruling out alternative scenarios to inflation. A very interesting property of
CPSC models is that the clock part of the signal encodes information about the
time-dependence of the scale factor itself, $a(t)\sim t^{p}$, during the early
Universe evolution [20]. The signal can be described in a model independent
way by the following template:

$\frac{\Delta\mathcal{P}_{\rm
clock}}{\mathcal{P}_{0}}=A\left(\frac{2k}{k_{r}}\right)^{-\frac{3}{2}+\frac{1}{2p}}\sin\left[p\,\frac{m_{\sigma}}{H}\left(\frac{2k}{k_{r}}\right)^{\frac{1}{p}}+\varphi\right],$
(5.1)

where $k_{r}$ is the first resonant mode and we see that the running of the
clock signal, i.e. the functional form of the argument of the oscillation, is
essentially the inverse function of the scale factor: $t(a)\sim a^{1/p}$ which
is $\sim k^{1/p}$ at horizon crossing for the mode $k$. The above template is
valid for $k<2\,k_{r}$ for $0<p<1$ and for $k>2\,k_{r}$ otherwise. We note
that inflation corresponds to the large $\lvert p|$ limit, where
$p(m_{\sigma}/H)(2k/k_{r})^{1/p}\rightarrow\delta\varphi+(m_{\sigma}/H)\ln(2k/k_{r})$,
where $\delta\varphi$ simply shifts the phase of the signal to
$\varphi\rightarrow\varphi+\delta\varphi$ and the oscillations become linear
in $\ln k$-space with a clock frequency $(m_{\sigma}/H)$. Ekpyrotic scenarios
correspond to $0<p<1/3$ and matter contraction to $p=3/2$.

Current Planck data does not show preference for any values of $p$ as no clock
signal is detected [29, 36]. It is therefore crucial to investigate whether we
will be able to do that and to rule out alternative scenarios to inflation
based on measurements of the clock signal. The procedure we adopt is the
following. We take our CPSC best fit. From our forecast analysis above, we
know that we will be able to find very strong to decisive evidence in favor of
this model if this is the true model. By looking at the reconstructed power
spectrum and at the CMB residuals, we can identify the beginning of the clock
as $\ell_{r}\sim 700$. Therefore, we fit the template in Eq. (5.1) to the
fiducial considering only $\ell>700$. The goal is to check whether we can get
a constraint on the parameter $p$. To make the analysis faster, and to avoid
degeneracies with the cosmological parameters that will be induced by cutting
$\ell<700$, only in this analysis we fix
$\omega_{c},\,\omega_{b},\,100\theta_{s},\,\tau$ to the CPSC best fit values,
but we do vary the amplitude and tilt of the power spectrum. For simplicity we
only explore the most optimistic experimental configuration, i.e. PL + LB +
S4. We use the following priors on the clock parameters (the priors on the
amplitude, tilt and frequency $m_{\sigma}/H$ are left unchanged):
$\log_{10}\,k_{r}\in[-1.25,\,0],\,$ $\varphi\in[0,\,2\pi],\,$
$A\in[0,\,0.3],\,$ $\log_{10}(\lvert p\rvert)\in[-3,\,2.7]$. We perform two
runs, one for positive and one for negative values of $p$. Since inflation
corresponds to large values of $\lvert p\rvert\gg 1$ and we want to exclude
values $\lvert p\rvert<1$, it is very useful to adopt a logarithmic prior on
$p$. However, it is complicated to do so for a parameter that switches sign.

The results are shown in Fig. 8, where we only plot the posteriors for the
relevant parameters, i.e. the frequency and amplitude of the signal together
with $p$. As can be seen, we will get a lower bound on the parameter $p$, i.e.
$\log_{10}\,p>0.328$ and $\log_{10}\,(-p)>0.323$ at 3$\sigma$ for positive and
negative values of $p$ respectively. As also shown in Fig. 8, this will rule
out alternative scenarios such as Ekpyrosis or matter contraction in a model
independent way which is complementary to the detection of B-modes,
traditionally considered the smoking gun of Inflation.

We also comment that the above method is suggestive but incomplete, because a
complete method of comparing CPSC signals of different scenarios should use
full predictions that include both sharp feature and clock signals. However,
for alternative scenarios to inflation, these results are unavailable to date.

##  6 Very similar features: the case for primordial non-Gaussianities

Figure 9:  Constraints on the Bump and Turn models assuming the bump fiducial
(top) and for the Turn and Dip models assuming the turn fiducial (bottom). We
also report the Bayes factors from our analysis in the bottom panel.

In the previous Section, using the fiducial CPSC model as a representative
example, we saw that the prospects for distinguishing features of different
nature are promising. As a final example, we address the following question:
can we distinguish, based only on the power spectra, between two feature
models that can produce very similar features at high $\ell$?

To answer that, we focus on the Bump/Dip and the Turn models. We remind the
reader that a dip in the single-field potential produces a sharp feature
signal which is very similar to that produced by a sharp turn in the field
space of a multifield model of inflation, except at the largest scales where
nevertheless its amplitude is small. It is therefore interesting to
investigate if the precision of future experiments is enough to distinguish
between those two models based only on such a small difference at large
scales. We also remind that the Bump model produces features with opposite
phase compared to the Dip and the Turn models.

There are therefore two interesting situations. First, if the Bump best fit is
the true model of our Universe, we should be able to rule out the Turn
model.888 Here we take Bump instead of Dip as the fiducial because, if we fit
Bump/Dip model to data, currently the data prefers Bump.  Second, if the Turn
best fit is the true model of our Universe, we expect the Dip model to fit it
quite well (again, except for the largest scales). In the following, we will
refer to these two cases as TB (Turn vs Bump) and DT (Dip vs Turn)
respectively. The results for TB and DT are shown in the top and bottom panels
of Fig. 9 respectively. In the same figure, we also provide the Bayes factors
from our analyses.

#### Turn v.s. Bump (TB).

Let us first discuss the TB case. The Bump model can be constrained very well,
and we see a detection of all the parameters. Interestingly, because a feature
is detected with very high significance, the parameters of the Turn model are
very well constrained, even though it cannot fit the data as well as the Bump
model and is eventually decisively excluded. The feature parameters of the
Turn model are constrained to (we take PL + LB + S4 constraints for
illustration) $N_{0}-N*=-1.712\pm 0.012,\,$ $\Delta N_{\rm
turn}=0.228^{+0.031}_{-0.11}\simeq\,\Delta N_{\rm
bump}\,\sqrt{2}^{+0.031}_{-0.11},\,$ $\Delta\mathcal{P}_{\rm
turn}/\mathcal{P}_{0}=0.0197^{+0.0048}_{-0.0064}$. We thus have a detection of
non-zero feature amplitude, although the mean value is smaller than the bump
one: $\Delta\mathcal{P}_{\rm bump}/\mathcal{P}_{0}=0.0325^{+0.0066}_{-0.016}$.
Note that, apart for the difference in the phase of the oscillations that we
already mentioned, another fact explains why the sharp feature from a turn is
unable to fit the one from the Bump best fit. The Bump best fit is described
by a very sharp feature, with $\Delta N=0.09$. For such a small value of
$\Delta N$, needed to match the envelope of the sinusoidal oscillations, the
Turn model develops a quite large bump at large scales, spoiling even more the
fit to the fiducial. Despite the hint to a best fit, the $\chi^{2}$ to the
data is much worse for the Turn model than for the Bump one. This, together
with a large fraction of the prior volume which ends up wasted, leads to a
very bad Bayes factor for the Turn model, and therefore the Bump model is
decisively favored. We thus provided one with an example where two models,
which are degenerate in a fraction of parameter space, can be told apart with
a very high statistical significance. Physically, our results imply that we
would be able to reconstruct a tiny bump in a single-field inflationary
potential, without the risk to confuse it with a sharp turn in a multi-
dimensional field space.

#### Dip v.s. Turn (DT).

We now turn to the second case, i.e DT. The situation is now very different
since the Dip model is able to fit virtually all the feature signal produced
by the Turn one in the primordial power spectrum, see the left panel of Fig.
10. More in detail, as can be seen in the lower panel of Fig. 9, both
amplitudes are constrained to a very similar value $\Delta\mathcal{P}_{\rm
turn}/\mathcal{P}_{0}=0.0197^{+0.0048}_{-0.0064}$ and $\Delta\mathcal{P}_{\rm
bump}/\mathcal{P}_{0}=-0.0167\pm 0.0042$, where again we quote results from PL
+ LB + S4. The mean values of $N_{0}$ and $\Delta N$ also match, keeping in
mind that to get the same feature envelope we must have $\Delta N_{\rm
turn}=\Delta N_{\rm bump}\,\sqrt{2}$. The two models will thus be
indistinguishable at the level of the power spectrum, as can also be seen from
the Bayes factors in Fig. 9. Physically, our results imply that, using only
the power spectrum, we would not be able to claim the detection of a sharp
turn in a multi-dimensional field space, as the same feature could be produced
by a tiny dip in a single-field inflationary potential.

Figure 10:  [Left] Power spectra for the best fit Turn and for a Dip model
matching the oscillations after the first one. They are the best fit from the
TD analysis. [Right] Associated $f_{\rm NL}^{\mathrm{eq}}(k)$ in the
equilateral limit.

#### Primordial non-Gaussianities.

Of course, it would be desirable to distinguish the Turn model from the Dip
one in order to draw conclusions about the micro-physical details of
inflation. Since we have just shown that this was impossible using only the
power spectrum, we turn to the comparison of the oscillatory features these
models leave in the primordial three-point function, the bispectrum (see,
e.g., Refs. [56, 19]). Using the transport approach implemented in the
PyTransport code [57, 58], we numerically compute the value of the bispectrum
for the two best fits. More precisely, the bispectrum
$B_{\zeta}(k_{1},k_{2},k_{3})$ is defined as
$\braket{\hat{\zeta}_{\vec{k}_{1}}\hat{\zeta}_{\vec{k}_{2}}\hat{\zeta}_{\vec{k}_{3}}}=(2\pi)^{3}\delta^{(3)}\left(\vec{k}_{1}+\vec{k}_{2}+\vec{k}_{3}\right)B_{\zeta}(k_{1},k_{2},k_{3})$,
and we vary the overall scale $k$ of the equilateral configuration
corresponding to $k_{1}=k_{2}=k_{3}=k$, quoting the dimensionless number
$f_{\mathrm{NL}}^{\mathrm{eq}}(k)=(10/9)\times
k^{6}B_{\zeta}(k,k,k)/\left[(2\pi)^{4}\mathcal{P}_{0}(k)^{2}\right]$. Note
that in this definition we have used the power-law power spectrum
$\mathcal{P}_{0}(k)$ containing the $n_{s}$ scale-dependence but not the
feature one. The scale-dependent quantities $f_{\mathrm{NL}}^{\mathrm{eq}}(k)$
for the Turn and the Dip best fits are represented in the right panel of Fig.
10. Two remarks are in order. First, the running of the signal in the
equilateral configuration of the bispectrum is very similar in both cases,
just like in the power spectrum. Second, the Dip best fit results in a
bispectrum with an amplitude roughly 3 times larger than the Turn one.
Although the $f_{\mathrm{NL}}^{\mathrm{eq}}$ of order unity found for the Turn
and Dip best fits will not be constrained with upcoming CMB experiments, these
results show that primordial non-Gaussianities may play a crucial role in
disentangling models of inflation that produce very similar features in the
power spectrum, if they are detectable by observations with higher
sensitivities. We leave for future works a more thorough analysis of
primordial non-Gaussianities in the context of the feature models and the
comparison methodology presented in this work, in order to correctly conclude
about the discriminating power of oscillatory features in higher-order
correlation functions in the CMB [59, 60, 61, 62]. In particular, the
information contained in the full shape of the bispectrum should be relevant,
beyond the exact equilateral configuration. In this respect, upcoming Large-
Scale-Structure experiments like DESI [63], Euclid [64] and SPHEREx [65], as
well as more futuristic hydrogen atom 21cm map experiments may bring a crucial
complementary constraining power.

##  7 An estimator for correlation between feature models

A = CPSC $(4.5\,\sigma)$
---
|  Forecast (as in Fig. 7) |  $\gamma_{A,B}\times$ 4.5$\,\sigma$
Resonant |  Between 2 and 3$\sigma$ |  3.2$\sigma$
Bump/Dip |  Between 1 and 2$\sigma$ |  1.8$\sigma$
Turn |  Between 1 and 2$\sigma$ |  2.7$\sigma$

A = Turn $(3.9\,\sigma)$
---
|  Forecast (as in Fig. 9) |  $\gamma_{A,B}\times$3.9$\,\sigma$
Dip |  3.9$\sigma$ |  3.9$\,\sigma$

Figure 11:  We plot the normalized $\Delta C_{\ell}^{XY}$. [Top-Left] we take
A to be the CPSC bestfit and B as the various best fits to the CPSC fiducial
from our forecast. [Top-Right] we take A to be the turn bestfit and B as the
bestfit of the bump model to the turn fiducial from our forecast. The colors
are the same as in Fig. 3, i.e., blue–CPSC, yellow–resonant, green–bump,
red–turn. In the lower tables, using the method of this section, we estimate
the detection levels of these models B based on the detection level of the
fiducial model A, and compare them to the exact detection level from our
forecast for PL+LB+S4. Since the $\gamma$s in the top panels are very similar
to each others, we simply estimate the detection levels by using the lowest
$\gamma$ for each model.

As we have seen from previous examples, some feature models share more
similarities than others. It may be useful to construct a simple estimator
quantifying how closely two models are correlated. If a new feature model is
constructed that is not analyzed above, one can use this estimator to give a
quick estimate of a lower bound of statistical significance based on its
correlation with the models analyzed in this work, without having to perform
the full data analysis.

We define a correlator $\gamma$ between two different feature profiles,
$\Delta\mathcal{P}_{A}$ and $\Delta\mathcal{P}_{B}$, as

$\gamma^{XY}_{A,B}\equiv\frac{(\Delta\mathcal{P}_{A},\Delta\mathcal{P}_{B})}{\sqrt{(\Delta\mathcal{P}_{A},\Delta\mathcal{P}_{A})(\Delta\mathcal{P}_{B},\Delta\mathcal{P}_{B})}}~{},$
(7.1)

where the scalar product $(\cdot\,,\,\cdot)$ is defined as

$(\Delta\mathcal{P}_{A},\Delta\mathcal{P}_{B})\equiv\int_{\ell_{\rm
min}}^{\ell_{\rm max}}\mathrm{d}\ell~{}w(\ell)~{}\frac{\Delta
C^{XY}_{\ell,A}}{C_{\ell}^{XY}}\frac{\Delta
C^{XY}_{\ell,B}}{C_{\ell}^{XY}}~{},$ (7.2)

in which $XY$ stands for $TT$, $TE$ or $EE$. $\Delta C_{\ell}^{XY}$ is
computed through the transfer function $T^{XY}_{\ell}(k)$,

$\Delta C_{\ell}^{XY}=\int\mathrm{d}\log
k~{}T^{XY}_{\ell}(k)\Delta\mathcal{P}~{}.$ (7.3)

Also note that, in the case $X,\,Y=T,\,E$, instead of considering
$C_{\ell}^{TE}$ at the denominator which vanishes at some multipoles, we use
the quantity
$\sqrt{C_{\ell}^{TE}\,C_{\ell}^{TE}+C_{\ell}^{TT}\,C_{\ell}^{EE}}$. The weight
function $w(\ell)$ is determined by the sensitivity of the experiment. Since
both the temperature and polarization will be eventually measured up to cosmic
variance limited precision, we will assume the cosmic variance limit when
defining the weight function. As the signal-to-noise ratio of the fractional
correction to the CMB power spectrum introduced by each feature is given by
$\left(\Delta
C_{\ell}^{XY}/C_{\ell}^{XY}\right)/\left(1/\sqrt{2\ell+1}\right)$, the weight
function can be chosen as

$w(\ell)=2\ell+1~{}.$ (7.4)

Several examples of this correlation estimator are shown in the top panels of
Fig. 11, where we also plot the quantities $\Delta
C^{XY}_{\ell,A,B}/C^{XY}_{\ell,A,B}$ for a better understanding of the
results. Combined with Fig. 7 and 9, from these examples we see that, in the
event of a feature detection, we will be able to distinguish two models using
CMB even if they are up to $\sim 80\%$ correlated, although we caution that
there may be model-dependent exceptions.

Another use of this estimator is the following. If a model A can be detected
in the CMB XY data with a statistical significance of $x^{XY}$-sigma and a new
model B, which has not been compared with data, is $\gamma_{A,B}^{XY}$
correlated with A, then one can estimate a rough lower bound of
$|\gamma^{XY}_{A,B}|~{}x^{XY}$-sigma detection for model B with the XY data.
This estimate is for the individual TT, TE or EE data. In the analyses of
previous sections, the constraints we got combined TT, TE and EE, so, without
knowing the individual constraints on model A, generally speaking, we do not
know how to get the combined estimate for model B. On the other hand, in a
simple situation where all three $\gamma_{A,B}^{XY}$’s are similar
$\sim\gamma_{A,B}$ (as in the several examples shown in Fig. 11), we can
approximate the lower bound as $|\gamma_{A,B}|~{}x$-sigma if the combined
significance for model A is $x$-sigma. Since we already have done the model
comparison in the CPSC and Turn model fiducial data, in the tables in Fig. 11
we test the values of $|\gamma_{A,B}|~{}x$-sigma for several models.

We have chosen to construct the estimator in terms of the CMB multipole space
$\ell$ instead of the momentum space $k$ of the original primordial feature
profile, because the signals in the CMB multipole space is what is relevant to
measurements. This choice is experiment specific; if we are concerned with
galaxy survey experiments instead of CMB, the scalar product can be defined in
the 3D momentum space,

$(\Delta\mathcal{P}_{A},\Delta\mathcal{P}_{B})=\int_{k_{\rm min}}^{k_{\rm
max}}\mathrm{d}k~{}w_{\rm LSS}(k)\frac{\Delta P_{A}}{P}(k)\,\frac{\Delta
P_{B}}{P}(k),$ (7.5)

with a different weight function $w_{\rm LSS}$.

##  8 Conclusions

Primordial features carry rich information on the primordial universe and can
potentially deliver answers to some of the most important questions about
inflation and alternative theories. These properties have been pointed out,
studied and classified in many recent works [15, 16, 17, 18]. Unfortunately,
theoretically, the amplitudes of such signals have no generic lower bound and
are highly model-dependent, and experimentally no statistically significant
signals have been detected. Thus, the future of this field critically depends
on the prospects of observational sensitivities for these signals in the near
future. In this paper, we addressed this question in the context of CMB E-mode
polarization data, which will soon be updated by experiments such as Simons
Observatory, LiteBIRD, and CMB S4. Due to their sharper transfer functions,
E-modes carry more information about features than T-modes and are therefore
expected to improve our knowledge about features.

In this work, we first summarized the theoretical predictions of several
representative feature models on the power spectrum using a unified formalism.
All these models currently give good fits to some small-scale glitches in the
Planck residual data. Then, using a state-of-the-art pipeline optimized to
compare feature models with data and using Bayesian inference, we showed that,
if the featureless model is the true model of the universe, these feature
candidates from the Planck data will be ruled out with strong evidences.
Conversely, if one of these feature models is the true model of the Universe,
the featureless one will be ruled out with very strong evidence.

Next, we addressed the main question of this paper. If one of the candidate
feature models is confirmed by upcoming CMB polarization data, we will have to
consider the possibility that this feature can be explained by several models
that have very different underlying physics. Will these data be strong enough
to tell the differences? We found that, due to the rich information that
feature-model predictions carry with their phases and envelops across all
scales, future observations will be able to distinguish them with very strong
to decisive Bayesian evidences. This result is even more impressive when we
realize that the feature models we considered are often highly correlated to
each others. However, there do exist exceptional cases in which
microscopically very different models predict almost identical features in the
power spectrum. Still, we showed for an example, that two models with
identical power spectra lead to primordial non-Gaussianities of different
amplitudes. A more thorough analysis of the bispectrum shape functions for
these very similar feature models, and the study of the detectability of these
non-Gaussian signals, are left for future study.

These findings clearly demonstrate that the study of primordial features
should be considered as one of the main tools to understand the primordial
universe in the next decades. With the addition of the large-scale structure
data from various galaxy surveys [66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76,
77, 78] in a similar time frame, as well as 21cm hydrogen line surveys [79,
80] and stochastic gravitational wave background mapping [81, 82, 37] in a
more distant future, we can hope for discoveries that will shake our
understanding of the origin of the Universe.

## Acknowledgments

MB is supported by the Spanish Atracción de Talento contract no.
2019-T1/TIC-13177 granted by Comunidad de Madrid, the I+D grant
PID2020-118159GA-C42 of the Spanish Ministry of Science and Innovation and the
i-LINK 2021 grant LINKA20416 of CSIC. MB and DH thanks the Indian Institute of
Technology Madras, where part of this work was carried out, for hospitality
and acknowledge travel support through the India-Italy mobility program
(INT/Italy/P-39/2022 (ER)). L.P. would like to acknowledge support from the
“Atracción de Talento” grant 2019-T1/TIC15784, his work is partially supported
by the Spanish Research Agency (Agencia Estatal de Investigación) through the
Grant IFT Centro de Excelencia Severo Ochoa No CEX2020-001007-S, funded by
MCIN/AEI/10.13039/501100011033.

##  Appendix A Methodology

To solve for the inflationary dynamics, we use BINGO and its two-field
extension. The details of the numerical methods are outlined for single-field
models in [43] and for two-field models in [44]. Note that while for single-
field models we integrate the perturbation equations from sub-Hubble to super-
Hubble scales, for two-field models we integrate till the $e$-fold when the
curvature perturbation freezes, to avoid any inaccuracies due to unaccounted
super-Hubble evolution. CAMB [83, 84] is used to calculate the angular power
spectrum. The primordial power spectrum, for single-field models, is evaluated
for all the modes required by CAMB without interpolation, while for two-field
models, we use an adaptive interpolation as outlined in [33]. We would like to
note that we evaluate angular power spectra at all multipoles as our priors in
the potential parameters allow for high-frequency oscillations that can be
missed if an interpolation is used in the multipole space.

For comparison with data and forecast we use CosmoChord (PolyChord add-on of
CosmoMC999https://cosmologist.info/cosmomc/ [85]) [86, 87]. We used the best
fits obtained against the 2020 CamSpec release
v12.5101010https://people.ast.cam.ac.uk/ stg20/camspec (hereafter EG20, from
the names of its authors [52]) datasets as this unbinned dataset covers larger
sky fraction in temperature and polarization with a different cleaning. In all
analyses, the nuisance parameters are allowed to vary and priors on the
nuisance parameters are used accordingly. The priors on the model parameters
are shown in Table 1.

The forecast pipeline follows the methodology discussed in our earlier
analysis in [35]. For completeness, we briefly discuss the pipeline here as
well. We work with the sensitivities of the upcoming and proposed CMB
observations. The multipole range for the fiducial spectra for the forecasts
are fixed to scales probed by Planck, i.e. between $\ell=2-2500$. Since the
best fits to the present data are obtained by comparing the models in this
multipole range, our choice ensures conservative forecasts and a comparison of
models within the same scales. Depending on the sensitivities, we combine the
present and upcoming observations spanning the multipole range. We consider
three combinations PL + SO, PL + LB +SO and PL + LB + S4. For the PL + SO
combination, we assume 80% sky coverage from Planck at large scales
($\ell=2-39$) and 20% at scales corresponding to $\ell=40-2500$. While Planck
already has data with significantly larger sky coverage at small scales, for
forecast, in order to avoid the overlap with SO we use a smaller fraction. We
assume SO to cover $\ell=40-2500$ with 40% sky fraction. In PL + LB + SO,
keeping the multipole range and sky coverage same for SO, we redistribute the
scales between Planck and LiteBIRD. PL noise is used here between
$\ell=800-2500$ (20% sky coverage) and LB noise is used for multipoles
$\ell=2-39$ (80% sky coverage) and $40-1350$ (20% sky coverage). For PL + LB +
S4 we use the same distributions in the multipole range with the SO noise and
sky coverage replaced by that of CMB-S4. Where the signal-to-noise ratios of
these observations are comparable, this distribution ensures insignificant
overlap between them.

For all our forecast analyses, an inverted Wishart likelihood has been used
following [88]. The noise power spectrum is calculated using the r.m.s. noise
per pixel, full width at half maxima of the assumed Gaussian beam and the sky
fraction. For SO, the baseline noise configuration of the large aperture
telescope has been used111111The temperature and polarization noise power
spectra for SO that we use in our analysis are available at
https://github.com/simonsobs.. Note that we have not used the B-mode noise as
we are not considering tensor modes.

Cosmological Parameters | Priors
---|---
$\Omega_{\mathrm{b}}h^{2}$ | $[0.02,\,0.0265]$
$\Omega_{\mathrm{CDM}}h^{2}$ | $[0.1,\,0.135]$
$100*\theta_{s}$ | $[1.03,\,1.05]$
$\tau$ | $[0.03,\,0.08]$

Feature parameters | Priors
---|---
| Resonant | CPSC | Turn | Bump/Dip
$\ln\,10^{10}A_{s}$ | $[2.763,\,4.375]$
$n_{s}$ | $[0.920,\,0.996]$
$\Delta\mathcal{P}/\mathcal{P}_{0}$ | $[0,\,0.35]$ | $[-0.35,\,0.35]$
$N_{0}-N_{*}$ | $-$ | $[-4.3,\,-1.0]$
$\log_{10}\,\omega/H$ | $[0,\,1.9]$ | $-$
$\Delta N$ | $-$ | $[0.05,\,1.5]$
$\varphi$ | $[0,\,2\pi]$ | $-$

Table 1:  Priors on the cosmological and feature parameters used in our
analysis in each models. The priors on foreground and calibration parameters
are kept unchanged from their default values in the Polychord analysis. Note
that, to unify the nomenclature, we define $\omega=m_{\sigma}$ for the CPSC
model.

To describe the statistical power of future CMB surveys to constrain the
parameter space of the models, we use the so called Figure of Merit (FoM),
defined as:

${\rm FoM}\equiv\lvert{\rm Fi}\rvert^{-1/2},$ (A.1)

where ${\rm Fi}$ is the covariance matrix obtained from our data analysis. An
increase in the FoM corresponds to a compression of the parameter space going
from a survey to another one, so a better constraining power of the latter. We
quote

$F\equiv\frac{{\rm FoM}_{\rm experiment}}{{\rm FoM}_{\rm Planck}}.$ (A.2)

##  Appendix B Computation of effective parameters

In this Appendix, we provide a detailed computation of the feature signals in
the power spectrum for the Bump/Dip and the Turn model using the in-in
formalism (see e.g. [15] for a review). We also compare the analytic
prediction with the numerical spectra to show to what extent they agree. The
purpose is to derive the effective parameters used in the analysis in the main
text. Since the effective parameters of the resonant and CPSC model have been
already derived elsewhere, we restrict our discussion to the sharp feature
models mentioned above.

We remind readers that the expressions of the effective parameters are
approximate relations that serve as an efficient bridge between observables
and raw model parameters. The purpose is to facilitate the search of the
parameter space, especially for models with a large number of parameters and
with complicated relations between observables and model parameters. All
theoretical predictions used in the data analyses and forecasts are numerical
results and not approximations. The full methodology is described in Ref.
[33].

###  B.1 Bump/dip in potential

This model is defined by $\Xi=0,\,\Delta
V/V=A\,\exp\left[-\left(\phi-\phi_{0}\right)^{2}/\Delta^{2}\right]$ (see Sec.
2). The signal will therefore depend on the width $\Delta$, location
$\phi_{0}$ and amplitude $A$ of the perturbation to the potential. The width
$\Delta$ can be easily related to the time duration of the sharp feature.
Assuming that the feature lasts $\Delta\phi=6\Delta$, which is approximately
$100\%$ of the area below the Gaussian perturbation to the potential, we
obtain $\Delta N=6\Delta/\sqrt{2\epsilon_{0}}$. On the other hand, the
relation between $\phi_{0}$ and $N_{0}$ is easily computed numerically with
the help of a shooting algorithm, see e.g. Refs. [33, 35] for details. In this
Section, we derive an approximate expression for the maximum amplitude of the
feature, which is the only remaining effective parameter that we need for our
data analysis. To do that, we use the in-in formalism. Here, we just sketch
the calculation referring the reader to reviews, see for example Ref. [15] for
the technicalities.

To compute the correction to the SFSR power spectrum, we need to separate the
second order Hamiltonian for the curvature perturbation, i.e.

$H_{2}(\epsilon)=\int\,\mathrm{d}^{3}\vec{x}\,a\,\epsilon\,M_{\mathrm{Pl}}^{2}\left[\zeta^{\prime
2}+\left(\partial\zeta\right)^{2}\right].$ (B.1)

into a free and interacting parts defined as $H_{2,\,0}\equiv
H_{2}(\epsilon_{0})$ and $H_{\text{int},2}\equiv H_{2}(\Delta\epsilon)$
respectively, where $\epsilon\equiv\epsilon_{0}+\Delta\epsilon$ and
$\Delta\epsilon$ is a small perturbation. During the bump/dip, we can ignore
the Hubble friction term in equation of motion for the inflaton and we have
the simpler energy conservation relation $\Delta(\dot{\phi}/2+V)=0$. With the
definition $\epsilon=\dot{\phi}^{2}/(2H^{2})$, we can easily compute
$\Delta\epsilon$ as:

$\Delta\epsilon=-3A\exp\left[-\frac{(\phi-\phi_{0})^{2}}{\Delta^{2}}\right]=-3A\exp\left[-\alpha(x-1)^{2}\right],$
(B.2)

where $\alpha=\frac{2\epsilon_{0}}{\Delta^{2}}=\frac{36}{\Delta N^{2}}$ and
$x\equiv\tau/\tau_{0}$, where $\tau_{0}$ is the conformal time of the sharp
feature. Using the in-in formalism, we compute the correction to the
primordial power spectrum as:

$\frac{\Delta\mathcal{P}(k)}{\mathcal{P}_{0}}=\frac{k^{3}}{2\pi^{2}}\frac{1}{\mathcal{P}_{0}}(-i)\int_{-\infty}^{0_{-}}{\rm
d}\tau a\,\langle 0\lvert\left[\zeta(k,\,0)\zeta(k,\,0),\,H_{{\rm
int},\,2}\right]\rvert 0\rangle,$ (B.3)

where the state $\lvert 0\rangle$ is the state annihilated by the annihilation
operators in the expansion of the curvature perturbation $\zeta$. Note that
$\zeta$ is in the interaction picture where operators are evolved using the
free Hamiltonian $H_{2,\,0}$.

Using the expression for $\Delta\epsilon$ found above, the integral becomes:

$\frac{\Delta\mathcal{P}(k)}{\mathcal{P}_{0}}=-\frac{3A}{\epsilon_{0}}{\rm
Im}\int^{\infty}_{0_{+}}{\rm
d}x\,\frac{1}{x^{2}}e^{-\alpha(x-1)^{2}}e^{2ixx_{0}}\left(2x_{0}x^{2}+2ix-\frac{1}{xx_{0}}\right)$
(B.4)

where $x_{0}\equiv\lvert k\tau_{0}\rvert$. The integral has a saddle point at
$x=1$. So we can simply evaluate the slow-varying terms inside the round
brackets at $x=1$ and take them out the integral. The remaining integral is
simply the integral of a Gaussian once we complete the square in the argument
of the exponential. The result is just:

$\frac{\Delta\mathcal{P}(x_{0})}{\mathcal{P}_{0}}=-\frac{3A}{\epsilon_{0}}\Delta\sqrt{\frac{\pi}{2\epsilon_{0}}}\left[2\cos
2x_{0}+(2x_{0}-\frac{1}{x_{0}})\sin 2x_{0}\right]e^{-x_{0}^{2}/\alpha}.$ (B.5)

Figure 12:  Comparison of the analytical formula Eq. (B.5).

Fig. 12 shows the comparison between the analytical formula and the numerical
results. Apart from some a slight discrepancy for $x_{0}\sim 1$, probably due
to the saddle point approximation we use, the two match very well.

Note that the maximum amplitude of the feature depends on $A$ through a simple
multiplicative factor. However the dependence on $\alpha$ is more complicated.
In fact, the spectrum peaks at different values depending on the value of
$\alpha$ through the relation $x_{0}^{\rm peak}=\sqrt{\alpha/2}$. Evaluating
the expression above at $x_{0}^{\rm peak}$, we get

$\frac{\Delta P_{\rm
max}}{P_{0}}=\frac{3A}{\epsilon_{0}}\sqrt{\frac{2\pi}{e}}.$ (B.6)

This is the effective parameter that we use in our analysis.

###  B.2 Turn model

This model is defined by
$\Xi=\xi\exp\left[-\left(\phi-\phi_{0}\right)^{2}/\Delta^{2}\right],\,\Delta
V/V=0$ (see Sec. 2). Also in this case, the signal will depend on the width
$\Delta$, location $\phi_{0}$ and amplitude $\xi$ of the perturbation to the
kinetic term of the inflaton. As before, the width $\Delta$ can be easily
related to the time duration of the sharp feature. Assuming that the feature
lasts $\Delta\phi=6\Delta$, we obtain again $\Delta
N=6\Delta/\sqrt{2\epsilon_{0}}$. To derive an approximate expression for the
maximum amplitude of the feature, using the in-in formalism, we integrate out
the massive field $\sigma$ and use an effectively single field description,
where the second order interacting Hamiltonian now contains a new term, in
addition to the one computed previously due to $\Delta\epsilon$:

$H_{{\rm
int},\,2}\subset\int\,\mathrm{d}^{3}\vec{x}\,a\,\epsilon_{0}\,M_{\mathrm{Pl}}^{2}\,\frac{\left(\partial\zeta\right)^{2}}{c_{s}^{2}},$
(B.7)

where $c_{s}^{2}$ is the squared speed of sound of the inflaton, which we
compute below as in [89, 50].

We start by deriving analytical predictions for the background. We set the
massive field initially at rest at the bottom of its potential $\sigma_{i}=0$.
Eventually, when $\phi\sim\phi_{0}$, it gets momentarily displaced from the
minimum. Neglecting the acceleration and friction terms in the Klein-Gordon
equation for $\sigma$, we get121212We define
$f\equiv\left[1+\sigma\Xi(\phi)\right]$.:

$V_{\sigma}=m_{\sigma}^{2}\sigma\simeq\Xi(\phi)f\dot{\phi}^{2}=2H^{2}\epsilon_{0}\Xi(\phi)/f$
(B.8)

where we used that $\epsilon_{0}=f^{2}\dot{\phi}^{2}/2H^{2}$. A crucial
quantity is the bending rate $\eta_{\perp}$, wich measures the acceleration
perpendicular to the background trajectory. Its expression is given by:

$H^{2}\eta_{\perp}^{2}=\frac{V_{s}^{2}}{\dot{\sigma}_{A}^{2}}$ (B.9)

where $\dot{\sigma}_{A}=\sqrt{f^{2}\dot{\phi}^{2}+\dot{\sigma}^{2}}$ and
$V_{s}\equiv e_{s}^{I}V_{I}$ where $I=\sigma,\,\phi$ and $e_{s}^{I}$ is the
entropic vielbein, i.e the unit vector orthogonal to the adiabatic vielbein.
In our model, the bending is given by:

$H^{2}\eta_{\perp}^{2}=\frac{1}{\dot{\sigma}_{A}^{4}}\left[V_{\sigma}^{2}f^{2}\dot{\phi}^{2}-2V_{\phi}V_{\sigma}\dot{\phi}\dot{\sigma}+V_{\sigma}^{2}\dot{\sigma}^{2}f^{-2}\right]\,.$
(B.10)

Using the slow-roll approximation for $\phi$, this becomes

$H^{2}\eta_{\perp}^{2}=\frac{f^{2}}{\dot{\phi}^{2}}\left[V_{\sigma}-3H\dot{\sigma}\right]^{2}.$
(B.11)

We can further assume that $3H\dot{\sigma}\ll V_{\sigma}$, to arrive to a very
simple expression131313Note that the assumption $3H\dot{\sigma}\ll V_{\sigma}$
is not guaranteed to be always satisfied, but it is so when the turn is not so
sharp and it lasts shorter than the typical oscillation period of the massive
field.:

$H^{2}\eta_{\perp}^{2}\simeq\frac{f^{2}}{\dot{\phi}^{2}}V_{\sigma}^{2}\simeq
2H^{2}\Xi(\phi)^{2}\epsilon_{0}f^{2}\simeq
2H^{2}\epsilon_{0}\xi^{2}\exp\left[-\frac{2\left(\phi-\phi_{0}\right)^{2}}{\Delta^{2}}\right].$
(B.12)

We investigate the agreement of this analytical expression with numerical
results in the top 3 panels of Fig. 13. As can be seen, Eq. (B.12) matches
quite well the numerics unless $\Delta N$ is very small. The analytical
approximation for the feature in the power spectrum that we are going to
present in the following does not apply to that case, for which the
approximations that we will assume break down. Nevertheless, we do consider
spectra that are not captured by our analytical approximation in our data
analysis, which is based on numerical integration.

In the single-field effective field theory (EFT), the turn in the trajectory
induces a sudden drop in the speed of sound of the inflaton. Integrating out
the entropic perturbation, well approximated by the $\sigma$ direction in
field space, $\zeta$ acquires a speed of sound given by:

$\frac{1}{c_{s}^{2}}-1=4\frac{H^{2}\eta_{\perp}^{2}}{m_{s}^{2}},$ (B.13)

where

$m_{s}^{2}=V_{ss}-H^{2}\eta_{\perp}^{2}+\epsilon H^{2}M_{\rm Pl}^{2}R_{\rm
fs}$ (B.14)

is the squared entropic mass and $R_{\rm fs}$ is the curvature of the field-
space metric, which is zero in our case $R_{\rm fs}=0$. Note that the last
term, which represents the curvature of the field space identically vanishes
in our model. The single-field effective theory is valid as long as
$m_{s}^{2}/H^{2}\gg 1/c_{s}^{2}$ and $\dot{c_{s}}/c_{s}\ll m_{s}$ (see e.g.
Refs. [89, 90, 91, 92]). These assumptions are very well satisfied in the
center and right panels of Fig. 13, but we have not checked their validity for
the left panel. In fact, in that case, the analytical expression for
$H\eta_{\perp}$ is drastically different from the numerical one, so, even if
the EFT held, we would still expect the analytical result obtained with (B.12)
not to describe the numerical one. The explicit expression for the speed of
sound in our model is:

$\displaystyle\frac{1}{c_{s}^{2}}-1=4\frac{H^{2}\eta_{\perp}^{2}}{m_{s}^{2}}$
$\displaystyle\simeq
8\left(\frac{\xi}{m_{\sigma}/H}\right)^{2}\epsilon_{0}\exp\left[-\frac{2\left(\phi-\phi_{0}\right)^{2}}{\Delta^{2}}\right]$
$\displaystyle\simeq
8\left(\frac{\xi}{m_{\sigma}/H}\right)^{2}\epsilon_{0}\exp\left[-\frac{4\epsilon_{0}\left(x-1\right)^{2}}{\Delta^{2}}\right].$
(B.15)

The correction to $\epsilon$ can be computed in a similar way and we get

$\left(\frac{1}{c_{s}^{2}}-1\right)=-4\frac{\Delta\epsilon}{\epsilon_{0}}$
(B.16)

The correction to the power spectrum is

$\frac{\Delta\mathcal{P}(k)}{\mathcal{P}_{0}}=\frac{k^{3}}{2\pi^{2}}\frac{1}{\mathcal{P}_{0}}(-i)\int_{-\infty}^{0_{-}}{\rm
d}\tau a\,\langle 0\lvert\left[\zeta(k,\,0)\zeta(k,\,0),\,H_{{\rm
2},\,\Delta\epsilon}+H_{{\rm 2},\,\Delta c_{s}}\right]\rvert 0\rangle,$ (B.17)

which can be easily computed as before, leading to:

$\frac{\Delta\mathcal{P}(x_{0})}{\mathcal{P}_{0}}=4\epsilon_{0}\sqrt{\frac{\pi}{\alpha}}\left(\frac{\xi}{m_{\sigma}/H}\right)^{2}\,\left[-\cos
2x_{0}+\left(x_{0}-\frac{1}{2x_{0}}\right)\sin
2x_{0}\right]e^{-x_{0}^{2}/\alpha},$ (B.18)

where now $\alpha=72/\Delta N^{2}=4\epsilon_{0}/\Delta^{2}$. As before, we
have that the maximum amplitude is at $x_{0}=\sqrt{\alpha/2}$. Approximating
trigonometric terms in the expression evaluated at $x_{0}=\sqrt{\alpha/2}$ to
1, the maximum amplitude thus becomes:

$\frac{\Delta\mathcal{P}}{\mathcal{P}_{0}}=4\epsilon_{0}\sqrt{\frac{\pi}{2e}}\left(\frac{\xi}{m_{\sigma}/H}\right)^{2}.$
(B.19)

As mentioned above, unless $\Delta N\ll 1$, the analytical solution matches
the numerics quite well, as can be seen by the bottom-center and right panels
of Fig. 13. Compared to the amplitude in the case of the bump/dip, there are
two main differences. First, $\alpha_{\rm turn}=2\alpha_{\rm bump/dip}$, so,
for the same $\Delta N$, the exponential suppression $\exp(-x_{0}^{2}/\alpha)$
is sharper for the bump/dip model. Second, while for the bump/dip model the
amplitude switches sign depending on whether we have a bump or a dip in the
potential, in the turn model the amplitude is proportional to $\xi^{2}$, so it
always has the same sign.

Figure 13:  Comparison of the analytical formulae Eq. (B.12) and (B.18).

##  Appendix C Estimate of feature signal in the CMB angular spectra

In this Appendix, we present a rough estimate of the effects of primordial
features on the CMB angular power spectra. Although the derivation below makes
use of several quite crude approximations, the results are very useful to
intuitively understand how primordial features are processed by the
convolution with the CMB transfer functions and how the $C_{\ell}$’s are
sensitive to different running of primordial features.

The CMB angular power spectra are obtained by convolving the primordial power
spectrum with transfer functions describing the late time evolution of the
Universe:

$C_{\ell}^{X\,Y}=\int\mathrm{d}\,\log k\,T^{XY}_{\ell}(k)\mathcal{P}(k)$ (C.1)

where $X,\,Y=T,\,E$ and $T^{XY}_{\ell}(k)$ are the transfer functions, which
are integrals along the line of sight of CMB photons of the source terms of
$T$ and $E$ modes. The transfer functions are function of $k$ and are labelled
by the multipole $\ell$. Some illustrative examples are shown in Fig. 14. The
transfer functions are strongly peaked at $\ell\simeq kD_{*}$, where $D_{*}$
is the comoving angular diameter distance from the observer to the redshift of
recombination. Furthermore, it is important to realize that the E-modes
transfer functions are more peaked than the T-modes ones. This is why E-modes
are so effective in constraining features, as we will see shortly.

To estimate the feature signal in the CMB spectra, we will make a series of
rough assumptions that are nevertheless reasonable for understanding the main
physics. All of them stem from the peaked shape of the transfer functions.

  * •

We will assume the one-to-one relation $\ell(k)=kD_{*}$ to be always valid.
Equivalently $k(\ell)=\ell/D_{*}$.

  * •

We model the transfer functions as top-hat functions in $\log k$ space as

$T^{XY}_{\ell}(k)=T_{\ell,\,{\rm max}}^{XY}\,\Theta\left[\log k-\left(\log
k(\ell)\,-\,\Delta_{\ell}^{XY}/2\right)\right]\,\Theta\left[-\log k+\left(\log
k(\ell)\,+\,\Delta_{\ell}^{XY}/2\right)\right],$

with an $\ell$-dependent width $\Delta_{\ell}^{XY}$, which is also different
for each pair $XY$. This assumption is quite reasonable as, again, the
transfer functions are strongly peaked around $\ell(k)$. We show the transfer
functions computed with the Einstein-Boltzmann solver for for $X,\,Y=T,\,T$
and their top-hat approximation for $\ell=2$ in the top panel of Fig. 14.

Roughly speaking, at a given $\ell$ (or $k$), if the variation of the
primordial spectrum is slow in the log-interval $\Delta_{\ell}^{XY}$, its
effect will show up in the angular spectrum $C_{\ell}^{XY}$. On the other
hand, if, within $\Delta_{\ell}^{XY}$, the PPS oscillates fast, the
oscillations will be washed out in the integral over $\ln k$. We now proceed
to give a more detailed proof of this naive argument.

Figure 14:  CMB transfer functions for T-modes (top) and E-modes (bottom).

For simplicity, we focus here on features with a constant envelope.
Nevertheless, as clear from the argument above, as long as the ($k$-dependent)
envelope is smooth enough over $\Delta_{\ell}$, the calculation below can be
straightforwardly generalized by simply multiplying our results by the
envelope. We parameterize the PPS as

$\mathcal{P}(k)=\mathcal{P}_{0}(k)(1+\delta_{\mathcal{P}})$ (C.2)

where $\mathcal{P}_{0}(k)=A_{s}\left(k/k_{*}\right)^{n_{s}-1}$ and

$\delta_{\mathcal{P}}(k)=C\begin{cases}\,\,\sin\omega\log\left(k/k_{p}\right)&\text{resonant}\\\
\,\,\sin k/k_{0}&\text{sharp}\end{cases},$ (C.3)

For simplicity, we neglect the phases and any $k$-dependent envelops of the
features in our calculation.

Our goal is to estimate the correction to the CMB spectra due to
$\delta_{\mathcal{P}}$. To do that, we write

$C_{\ell}^{X\,Y}\equiv
C_{\ell}^{X\,Y,\,0}(1+\delta_{\ell})=\int\mathrm{d}\,\log
k\,T^{XY}_{\ell}(k)\mathcal{P}_{0}(k)(1+\delta_{\mathcal{P}})$ (C.4)

The calculation of the leading order contribution is very simple. Since
$\mathcal{P}_{0}(k)$ is almost constant, it can be factored out the integral
and we are left with:

$\displaystyle C_{\ell}^{X\,Y,\,0}$ $\displaystyle=\int\mathrm{d}\,\log
k\,T^{XY}_{\ell}(k)\mathcal{P}_{0}(k)$ $\displaystyle=T_{\ell,\,{\rm
max}}^{XY}\mathcal{P}_{0}(\ell)\int\mathrm{d}x\,\Theta\left[x-\left(\log
k(\ell)\,-\,\Delta_{\ell}^{XY}/2\right)\right]\,\Theta\left[-x+\left(\log
k(\ell)\,+\,\Delta_{\ell}^{XY}/2\right)\right]$ $\displaystyle=T_{\ell,\,{\rm
max}}^{XY}\,A_{s}\left(\frac{\ell}{\ell_{*}}\right)^{n_{s}-1}\Delta_{\ell}^{XY}.$
(C.5)

For example, at large scales, we have the very well known Sachs-Wolfe plateau
for which $T_{\ell,\,{\rm max}}^{TT}\,\Delta_{\ell}^{TT}\propto
1/\ell(\ell+1)$. Anyway, the dominant contribution $C^{0}$ is not relevant for
our discussion.

We are interested in $\delta_{\ell}$. We have

$\displaystyle C_{\ell}^{X\,Y,\,0}\,\delta_{\ell}$
$\displaystyle=T_{\ell,\,{\rm
max}}^{XY}\mathcal{P}_{0}(\ell)\int\mathrm{d}\,\log k\,\cdots$
$\displaystyle=T_{\ell,\,{\rm
max}}^{XY}\mathcal{P}_{0}(\ell)\Delta_{\ell}^{XY}\frac{1}{\Delta_{\ell}^{XY}}\int\mathrm{d}\,\log
k\,\cdots$
$\displaystyle=C_{\ell}^{X\,Y,\,0}\frac{1}{\Delta_{\ell}^{XY}}\int\mathrm{d}\,\log
k\,\cdots,$ (C.6)

so

$\delta_{\ell}=\frac{C}{\Delta_{\ell}^{XY}}\int^{\log k_{\ell}+\Delta/2}_{\log
k_{\ell}-\Delta/2}\mathrm{d}x\,\,\begin{cases}\sin(\omega
x)&\text{resonant},\\\ \sin(e^{x}/k_{0})&\text{sharp}.\end{cases},$ (C.7)

where, for simplicity, we forgot the $X,\,Y,\,\ell$ labels on $\Delta$ and
denoted $k_{\ell}\equiv k(\ell)$. The effects of sharp and resonant features
with constant envelope $C=0.1$ are shown in Fig. 15. In the following
subsections we discuss each of them separately. Let us stress that, in the
case where the features have a $k$-dependent envelope, the result above still
applies by substituting $C$ with $C(k)$.

Figure 15:  Residuals due to the $\sin\log\,k$ (left) and $\sin\,k$ (right)
features in the PPS. Note that we keep the cosmological parameters fixed, so
the correction is only due the features in the PPS.

###  C.1 Resonant features: $\sin\log k$

Solving the integral (C.7), we obtain:

$\delta_{\ell}^{\rm res}\simeq
C\sin\left[\omega\log\,k(\ell)/k_{p}\right]\,\frac{\sin\left(\omega\Delta_{\ell}\,/2\right)}{\omega\Delta_{\ell}\,/2}.$
(C.8)

In the limit where the transfer functions are very narrow, i.e.
$\omega\Delta_{\ell}\to 0$ at all multipoles, the correction is simply the
same as the correction to the PPS $\delta_{\ell}^{\rm res}\simeq
C\sin\left(\omega\log\,k(\ell)/k_{p}\right)$. This can also be obtained by
simply setting the transfer function in the original integral as a Dirac
delta. In Fig. 15, we show an example that satisfies the condition
$\omega\,\Delta_{\ell}/2\ll 1$ with $\omega=4$ and $X,\,Y=E,\,E$. See the
orange line of the bottom-left panel of Fig. 15. Note that the plots show the
numerical result for $\delta_{\ell}$ obtained with the Einstein-Boltzmann
solver.

If the limit $\omega\,\Delta_{\ell}/2\ll 1$ is not satisfied, as for all the
other cases in the left panel of the figure, Eq. C.8 can still be used to
understand the impact of resonant features on the CMB spectra. Remember that,
at a given $\ell$, $\Delta^{TT}_{\ell}>\Delta^{EE}_{\ell}$ and
$\Delta^{XY}_{\ell}$ decreases with $\ell$. Therefore, the amplitude of the
feature is less washed out in $EE$; and, at a given feature frequency
$\omega$, $\delta_{\ell}$ increases with $\ell$. If we instead fix $\ell$, we
see that $\delta_{\ell}$ is larger for smaller frequencies $\omega$.

To summarize, the sensitivity to resonant features:

  * •

increases with $\ell$

  * •

decreases with $\omega$.

As we will see in the next subsection, the first of this conclusions does not
apply to sharp features, for which, in fact, the opposite holds.

###  C.2 Sharp features: $\sin k$

Solving the integral (C.7), we obtain:

$\delta_{\ell}^{\rm sharp}\simeq C\frac{{\rm
Si}\left(k\,\frac{e^{\Delta_{\ell}/2}}{k_{0}}\right)-{\rm
Si}\left(k\,\frac{e^{-\Delta_{\ell}/2}}{k_{0}}\right)}{\omega\Delta_{\ell}\,/2}.$
(C.9)

where ${\rm Si}$ is the sine integral function. Also here, in the limit where
the transfer functions are very narrow, i.e. $\Delta_{\ell}\to 0$ at all
multipoles, the correction is simply the same as the correction to the PPS
$\delta_{\ell}^{\rm res}\simeq C\sin\left(k(\ell)/k_{0}\right)$. Some examples
are shown in the right panel of Fig. 15.

In the limit of large frequencies $\omega_{\rm lin}\equiv k_{0}^{-1}\ll 1$ we
can expand the formula above as

$\delta_{\ell}^{\rm
sharp}\simeq\frac{C}{\Delta_{\ell}\,k/k_{0}}\left[e^{\Delta_{\ell}/2}\cos\left(e^{-\Delta_{\ell}/2}k/k_{0}\right)-e^{-\Delta_{\ell}/2}\cos\left(e^{\Delta_{\ell}/2}k/k_{0}\right)\right].$
(C.10)

The envelope of $\delta_{\ell}$ is very well approximated by $\sim
2e^{\Delta_{\ell}/2}k_{0}/k\Delta_{\ell}=2e^{\Delta_{\ell}/2}/\Delta_{\ell}\omega_{\rm
lin}\,k$. The consequence is that $\delta_{\ell}$ decreases at smaller scales
because of the suppression factor $k^{-1}$. Note that $\delta_{\ell}$
increases too because $\Delta_{\ell}^{-1}$ grows for large wavenumbers $k$,
but it does not grow enough to compensate for the $k^{-1}$ factor. Furthermore
$\delta_{\ell}$ is suppressed for large values of $\omega_{\rm lin}$, as in
the case for resonant features (see above).

To summarize, the sensitivity to sharp features:

  * •

decreases with $\ell$

  * •

decreases with $\omega$.

## References

  * [1] Planck collaboration, N. Aghanim et al., _Planck 2018 results. I. Overview and the cosmological legacy of Planck_ , _Astron. Astrophys._ 641 (2020) A1 [1807.06205].
  * [2] Planck collaboration, Y. Akrami et al., _Planck 2018 results. X. Constraints on inflation_ , _Astron. Astrophys._ 641 (2020) A10 [1807.06211].
  * [3] A. A. Starobinsky, _Spectrum of relict gravitational radiation and the early state of the universe_ , _JETP Lett._ 30 (1979) 682.
  * [4] BICEP, Keck collaboration, P. A. R. Ade et al., _Improved Constraints on Primordial Gravitational Waves using Planck, WMAP, and BICEP/Keck Observations through the 2018 Observing Season_ , _Phys. Rev. Lett._ 127 (2021) 151301 [2110.00483].
  * [5] P. Campeti and E. Komatsu, _New Constraint on the Tensor-to-scalar Ratio from the Planck and BICEP/Keck Array Data Using the Profile Likelihood_ , _Astrophys. J._ 941 (2022) 110 [2205.05617].
  * [6] M. Tristram et al., _Improved limits on the tensor-to-scalar ratio using BICEP and Planck data_ , _Phys. Rev. D_ 105 (2022) 083524 [2112.07961].
  * [7] D. Paoletti, F. Finelli, J. Valiviita and M. Hazumi, _Planck and BICEP/Keck Array 2018 constraints on primordial gravitational waves and perspectives for future B-mode polarization measurements_ , 2208.10482.
  * [8] A. A. Starobinsky, _A New Type of Isotropic Cosmological Models Without Singularity_ , _Phys. Lett. B_ 91 (1980) 99.
  * [9] R. Kallosh, A. Linde and D. Roest, _Superconformal Inflationary $\alpha$-Attractors_, _JHEP_ 11 (2013) 198 [1311.0472].
  * [10] J. J. M. Carrasco, R. Kallosh, A. Linde and D. Roest, _Hyperbolic geometry of cosmological attractors_ , _Phys. Rev. D_ 92 (2015) 041301 [1504.05557].
  * [11] R. Kallosh and A. Linde, _B-mode Targets_ , _Phys. Lett. B_ 798 (2019) 134970 [1906.04729].
  * [12] Simons Observatory collaboration, P. Ade et al., _The Simons Observatory: Science goals and forecasts_ , _JCAP_ 02 (2019) 056 [1808.07445].
  * [13] LiteBIRD collaboration, M. Hazumi et al., _LiteBIRD: JAXA’s new strategic L-class mission for all-sky surveys of cosmic microwave background polarization_ , _Proc. SPIE Int. Soc. Opt. Eng._ 11443 (2020) 114432F [2101.12449].
  * [14] CMB-S4 collaboration, K. N. Abazajian et al., _CMB-S4 Science Book, First Edition_ , 1610.02743.
  * [15] X. Chen, _Primordial Non-Gaussianities from Inflation Models_ , _Adv. Astron._ 2010 (2010) 638979 [1002.1416].
  * [16] J. Chluba, J. Hamann and S. P. Patil, _Features and New Physical Scales in Primordial Observables: Theory and Observation_ , _Int. J. Mod. Phys. D_ 24 (2015) 1530023 [1505.01834].
  * [17] A. Slosar et al., _Scratches from the Past: Inflationary Archaeology through Features in the Power Spectrum of Primordial Fluctuations_ , _Bull. Am. Astron. Soc._ 51 (2019) 98 [1903.09883].
  * [18] A. Achúcarro et al., _Inflation: Theory and Observations_ , 2203.08128.
  * [19] X. Chen, R. Easther and E. A. Lim, _Generation and Characterization of Large Non-Gaussianities in Single Field Inflation_ , _JCAP_ 04 (2008) 010 [0801.3295].
  * [20] X. Chen, _Primordial Features as Evidence for Inflation_ , _JCAP_ 01 (2012) 038 [1104.1323].
  * [21] X. Chen, _Fingerprints of Primordial Universe Paradigms as Features in Density Perturbations_ , _Phys. Lett. B_ 706 (2011) 111 [1106.1635].
  * [22] X. Chen and C. Ringeval, _Searching for Standard Clocks in the Primordial Universe_ , _JCAP_ 08 (2012) 014 [1205.6085].
  * [23] T. Battefeld, J. C. Niemeyer and D. Vlaykov, _Probing Two-Field Open Inflation by Resonant Signals in Correlation Functions_ , _JCAP_ 05 (2013) 006 [1302.3877].
  * [24] X. Gao, D. Langlois and S. Mizuno, _Oscillatory features in the curvature power spectrum after a sudden turn of the inflationary trajectory_ , _JCAP_ 10 (2013) 023 [1306.5680].
  * [25] T. Noumi and M. Yamaguchi, _Primordial spectra from sudden turning trajectory_ , _JCAP_ 12 (2013) 038 [1307.7110].
  * [26] R. Saito, M. Nakashima, Y.-i. Takamizu and J. Yokoyama, _Resonant Signatures of Heavy Scalar Fields in the Cosmic Microwave Background_ , _JCAP_ 11 (2012) 036 [1206.2164].
  * [27] R. Saito and Y.-i. Takamizu, _Localized Features in Non-Gaussianity from Heavy Physics_ , _JCAP_ 06 (2013) 031 [1303.3839].
  * [28] X. Chen and M. H. Namjoo, _Standard Clock in Primordial Density Perturbations and Cosmic Microwave Background_ , _Phys. Lett. B_ 739 (2014) 285 [1404.1536].
  * [29] X. Chen, M. H. Namjoo and Y. Wang, _Models of the Primordial Standard Clock_ , _JCAP_ 1502 (2015) 027 [1411.2349].
  * [30] Q.-G. Huang and S. Pi, _Power-law modulation of the scalar power spectrum from a heavy field with a monomial potential_ , _JCAP_ 04 (2018) 001 [1610.00115].
  * [31] G. Domènech, J. Rubio and J. Wons, _Mimicking features in alternatives to inflation with interacting spectator fields_ , _Phys. Lett. B_ 790 (2019) 263 [1811.08224].
  * [32] Y. Wang, Z. Wang and Y. Zhu, _Non-standard primordial clocks from induced mass in alternative to inflation scenarios_ , _JCAP_ 11 (2020) 026 [2007.09677].
  * [33] M. Braglia, X. Chen and D. K. Hazra, _Comparing multi-field primordial feature models with the Planck data_ , _JCAP_ 06 (2021) 005 [2103.03025].
  * [34] M. Braglia, X. Chen and D. K. Hazra, _Uncovering the history of cosmic inflation from anomalies in cosmic microwave background spectra_ , _Eur. Phys. J. C_ 82 (2022) 498 [2106.07546].
  * [35] M. Braglia, X. Chen and D. K. Hazra, _Primordial standard clock models and CMB residual anomalies_ , _Phys. Rev. D_ 105 (2022) 103523 [2108.10110].
  * [36] J. Hamann and J. Wons, _Optimising inflationary features the Bayesian way_ , _JCAP_ 03 (2022) 036 [2112.08571].
  * [37] A. Bodas and R. Sundrum, _Primordial Clocks within Stochastic Gravitational Wave Anisotropies_ , 2205.04482.
  * [38] X. Chen, R. Ebadi and S. Kumar, _Classical cosmological collider physics and primordial features_ , _JCAP_ 08 (2022) 083 [2205.01107].
  * [39] A. Antony and S. Jain, _Effect of damped oscillations in the inflationary potential_ , _Eur. Phys. J. C_ 82 (2022) 687 [2110.06837].
  * [40] D. K. Hazra, A. Antony and A. Shafieloo, _One spectrum to cure them all: signature from early Universe solves major anomalies and tensions in cosmology_ , _JCAP_ 08 (2022) 063 [2201.12000].
  * [41] A. Antony, F. Finelli, D. K. Hazra and A. Shafieloo, _Discordances in cosmology and the violation of slow-roll inflationary dynamics_ , 2202.14028.
  * [42] CORE collaboration, F. Finelli et al., _Exploring cosmic origins with CORE: Inflation_ , _JCAP_ 04 (2018) 016 [1612.08270].
  * [43] D. K. Hazra, L. Sriramkumar and J. Martin, _BINGO: A code for the efficient computation of the scalar bi-spectrum_ , _JCAP_ 05 (2013) 026 [1201.0926].
  * [44] M. Braglia, D. K. Hazra, L. Sriramkumar and F. Finelli, _Generating primordial features at large scales in two field models of inflation_ , _JCAP_ 08 (2020) 025 [2004.00672].
  * [45] WMAP collaboration, H. V. Peiris et al., _First year Wilkinson Microwave Anisotropy Probe (WMAP) observations: Implications for inflation_ , _Astrophys. J. Suppl._ 148 (2003) 213 [astro-ph/0302225].
  * [46] J. A. Adams, B. Cresswell and R. Easther, _Inflationary perturbations from a potential with a step_ , _Phys. Rev. D_ 64 (2001) 123514 [astro-ph/0102236].
  * [47] A. A. Starobinsky, _Spectrum of adiabatic perturbations in the universe when there are singularities in the inflation potential_ , _JETP Lett._ 55 (1992) 489.
  * [48] D. K. Hazra, A. Shafieloo, G. F. Smoot and A. A. Starobinsky, _Wiggly Whipped Inflation_ , _JCAP_ 08 (2014) 048 [1405.2012].
  * [49] D. K. Hazra, D. Paoletti, I. Debono, A. Shafieloo, G. F. Smoot and A. A. Starobinsky, _Inflation story: slow-roll and beyond_ , _JCAP_ 12 (2021) 038 [2107.09460].
  * [50] A. Achucarro, J.-O. Gong, S. Hardeman, G. A. Palma and S. P. Patil, _Features of heavy physics in the CMB power spectrum_ , _JCAP_ 01 (2011) 030 [1010.3693].
  * [51] X. Gao, D. Langlois and S. Mizuno, _Influence of heavy modes on perturbations in multiple field inflation_ , _JCAP_ 10 (2012) 040 [1205.5275].
  * [52] G. Efstathiou and S. Gratton, _A Detailed Description of the CamSpec Likelihood Pipeline and a Reanalysis of the Planck High Frequency Maps_ , 1910.00483.
  * [53] W. Handley, _fgivenx: Functional posterior plotter_ , _The Journal of Open Source Software_ 3 (2018) .
  * [54] L. T. Hergt, W. J. Handley, M. P. Hobson and A. N. Lasenby, _Constraining the kinetically dominated Universe_ , _Phys. Rev. D_ 100 (2019) 023501 [1809.07737].
  * [55] W. J. Handley, A. N. Lasenby, H. V. Peiris and M. P. Hobson, _Bayesian inflationary reconstructions from Planck 2018 data_ , _Phys. Rev. D_ 100 (2019) 103511 [1908.00906].
  * [56] X. Chen, R. Easther and E. A. Lim, _Large Non-Gaussianities in Single Field Inflation_ , _JCAP_ 06 (2007) 023 [astro-ph/0611645].
  * [57] M. Dias, J. Frazer, D. J. Mulryne and D. Seery, _Numerical evaluation of the bispectrum in multiple field inflation—the transport approach with code_ , _JCAP_ 12 (2016) 033 [1609.00379].
  * [58] D. J. Mulryne and J. W. Ronayne, _PyTransport: A Python package for the calculation of inflationary correlation functions_ , _J. Open Source Softw._ 3 (2018) 494 [1609.00381].
  * [59] J. R. Fergusson and E. P. S. Shellard, _The shape of primordial non-Gaussianity and the CMB bispectrum_ , _Phys. Rev. D_ 80 (2009) 043510 [0812.3413].
  * [60] J. R. Fergusson, H. F. Gruetjen, E. P. S. Shellard and M. Liguori, _Combining power spectrum and bispectrum measurements to detect oscillatory features_ , _Phys. Rev. D_ 91 (2015) 023502 [1410.5114].
  * [61] J. R. Fergusson, H. F. Gruetjen, E. P. S. Shellard and B. Wallisch, _Polyspectra searches for sharp oscillatory features in cosmic microwave sky data_ , _Phys. Rev. D_ 91 (2015) 123506 [1412.6152].
  * [62] Planck collaboration, Y. Akrami et al., _Planck 2018 results. IX. Constraints on primordial non-Gaussianity_ , _Astron. Astrophys._ 641 (2020) A9 [1905.05697].
  * [63] DESI collaboration, A. Aghamousa et al., _The DESI Experiment Part I: Science,Targeting, and Survey Design_ , 1611.00036.
  * [64] L. Amendola et al., _Cosmology and fundamental physics with the Euclid satellite_ , _Living Rev. Rel._ 21 (2018) 2 [1606.00180].
  * [65] O. Doré et al., _Cosmology with the SPHEREX All-Sky Spectral Survey_ , 1412.4872.
  * [66] Z. Huang, L. Verde and F. Vernizzi, _Constraining inflation with future galaxy redshift surveys_ , _JCAP_ 04 (2012) 005 [1201.5955].
  * [67] D. K. Hazra, _Changes in the halo formation rates due to features in the primordial spectrum_ , _JCAP_ 03 (2013) 003 [1210.7170].
  * [68] X. Chen, C. Dvorkin, Z. Huang, M. H. Namjoo and L. Verde, _The Future of Primordial Features with Large-Scale Structure Surveys_ , _JCAP_ 11 (2016) 014 [1605.09365].
  * [69] M. Ballardini, F. Finelli, C. Fedeli and L. Moscardini, _Probing primordial features with future galaxy surveys_ , _JCAP_ 10 (2016) 041 [1606.03747].
  * [70] G. A. Palma, D. Sapone and S. Sypsas, _Constraints on inflation with LSS surveys: features in the primordial power spectrum_ , _JCAP_ 06 (2018) 004 [1710.02570].
  * [71] B. L’Huillier, A. Shafieloo, D. K. Hazra, G. F. Smoot and A. A. Starobinsky, _Probing features in the primordial perturbation spectrum with large-scale structure data_ , _Mon. Not. Roy. Astron. Soc._ 477 (2018) 2503 [1710.10987].
  * [72] M. Ballardini, F. Finelli, R. Maartens and L. Moscardini, _Probing primordial features with next-generation photometric and radio surveys_ , _JCAP_ 04 (2018) 044 [1712.07425].
  * [73] F. Beutler, M. Biagetti, D. Green, A. Slosar and B. Wallisch, _Primordial Features from Linear to Nonlinear Scales_ , _Phys. Rev. Res._ 1 (2019) 033209 [1906.08758].
  * [74] M. Ballardini, R. Murgia, M. Baldi, F. Finelli and M. Viel, _Non-linear damping of superimposed primordial oscillations on the matter power spectrum in galaxy surveys_ , _JCAP_ 04 (2020) 030 [1912.12499].
  * [75] I. Debono, D. K. Hazra, A. Shafieloo, G. F. Smoot and A. A. Starobinsky, _Constraints on features in the inflationary potential from future Euclid data_ , _Mon. Not. Roy. Astron. Soc._ 496 (2020) 3448 [2003.05262].
  * [76] Y. Li, H.-M. Zhu and B. Li, _Non-linear reconstruction of features in the primordial power spectrum from large-scale structure_ , _Mon. Not. Roy. Astron. Soc._ 514 (2022) 4363 [2102.09007].
  * [77] M. Ballardini, F. Finelli, F. Marulli, L. Moscardini and A. Veropalumbo, _New constraints on primordial features from the galaxy two-point correlation function_ , 2202.08819.
  * [78] D. Chandra and S. Pal, _Investigating the constraints on primordial features with future cosmic microwave background and galaxy surveys_ , _JCAP_ 09 (2022) 024 [2205.01164].
  * [79] X. Chen, P. D. Meerburg and M. Münchmeyer, _The Future of Primordial Features with 21 cm Tomography_ , _JCAP_ 09 (2016) 023 [1605.09364].
  * [80] Y. Xu, J. Hamann and X. Chen, _Precise measurements of inflationary features with 21 cm observations_ , _Phys. Rev. D_ 94 (2016) 123518 [1607.00817].
  * [81] J. Fumagalli, S. Renaux-Petel and L. T. Witkowski, _Oscillations in the stochastic gravitational wave background from sharp features and particle production during inflation_ , _JCAP_ 08 (2021) 030 [2012.02761].
  * [82] M. Braglia, X. Chen and D. K. Hazra, _Probing Primordial Features with the Stochastic Gravitational Wave Background_ , _JCAP_ 03 (2021) 005 [2012.05821].
  * [83] A. Lewis, A. Challinor and A. Lasenby, _Efficient computation of CMB anisotropies in closed FRW models_ , _Astrophys. J._ 538 (2000) 473 [astro-ph/9911177].
  * [84] C. Howlett, A. Lewis, A. Hall and A. Challinor, _CMB power spectrum parameter degeneracies in the era of precision cosmology_ , _JCAP_ 04 (2012) 027 [1201.3654].
  * [85] A. Lewis and S. Bridle, _Cosmological parameters from CMB and other data: A Monte Carlo approach_ , _Phys. Rev. D_ 66 (2002) 103511 [astro-ph/0205436].
  * [86] W. J. Handley, M. P. Hobson and A. N. Lasenby, _PolyChord: nested sampling for cosmology_ , _Mon. Not. Roy. Astron. Soc._ 450 (2015) L61 [1502.01856].
  * [87] W. J. Handley, M. P. Hobson and A. N. Lasenby, _polychord: next-generation nested sampling_ , _Mon. Not. Roy. Astron. Soc._ 453 (2015) 4385 [1506.00171].
  * [88] S. Hamimeche and A. Lewis, _Likelihood Analysis of CMB Temperature and Polarization Power Spectra_ , _Phys. Rev. D_ 77 (2008) 103013 [0801.0554].
  * [89] A. J. Tolley and M. Wyman, _The Gelaton Scenario: Equilateral non-Gaussianity from multi-field dynamics_ , _Phys. Rev. D_ 81 (2010) 043502 [0910.1853].
  * [90] S. Cespedes, V. Atal and G. A. Palma, _On the importance of heavy fields during inflation_ , _JCAP_ 05 (2012) 008 [1201.4848].
  * [91] S. Garcia-Saenz, L. Pinol and S. Renaux-Petel, _Revisiting non-Gaussianity in multifield inflation with curved field space_ , _JHEP_ 01 (2020) 073 [1907.10403].
  * [92] L. Pinol, _Multifield inflation beyond $N_{\mathrm{field}}=2$: non-Gaussianities and single-field effective theory_, _JCAP_ 04 (2021) 002 [2011.05930].