Euclid: The search for primordial features

Primordial features, in particular oscillatory signals, imprinted in the primordial power spectrum of density perturbations represent a clear window of opportunity for detecting new physics at high-energy scales. Future spectroscopic and photometric measurements from the $Euclid$ space mission will provide unique constraints on the primordial power spectrum, thanks to the redshift coverage and high-accuracy measurement of nonlinear scales, thus allowing us to investigate deviations from the standard power-law primordial power spectrum. We consider two models with primordial undamped oscillations superimposed on the matter power spectrum, one linearly spaced in $k$-space the other logarithmically spaced in $k$-space. We forecast uncertainties applying a Fisher matrix method to spectroscopic galaxy clustering, weak lensing, photometric galaxy clustering, cross correlation between photometric probes, spectroscopic galaxy clustering bispectrum, CMB temperature and $E$-mode polarization, temperature-polarization cross correlation, and CMB weak lensing. We also study a nonlinear density reconstruction method to retrieve the oscillatory signals in the primordial power spectrum. We find the following percentage relative errors in the feature amplitude with $Euclid$ primary probes for the linear (logarithmic) feature model: 21% (22%) in the pessimistic settings and 18% (18%) in the optimistic settings at 68.3% confidence level (CL) using GC$_{\rm sp}$+WL+GC$_{\rm ph}$+XC. Combining all the sources of information explored expected from $Euclid$ in combination with future SO-like CMB experiment, we forecast ${\cal A}_{\rm lin} \simeq 0.010 \pm 0.001$ at 68.3% CL and ${\cal A}_{\rm log} \simeq 0.010 \pm 0.001$ for GC$_{\rm sp}$(PS rec + BS)+WL+GC$_{\rm ph}$+XC+SO-like both for the optimistic and pessimistic settings over the frequency range $(1,\,10^{2.1})$.


Introduction
Future galaxy surveys are the new frontier in the study of initial conditions of the Universe.They are expected to improve signif-icantly the constraints on many of the parameters characterising the physics of the early Universe.The Euclid satellite, simultaneously performing a spectroscopic survey of galaxies and an imaging survey (targeting weak lensing and galaxy clustering using photometric redshifts), will have the unique opportunity of measuring ultra-large scales thanks to its large observed volume, as well as small scales of matter distribution in the full nonlinear regime.This opportunity will drastically improve our understanding of the early-Universe physics and of cosmic inflation (Starobinsky 1980;Guth 1981;Sato 1981;Linde 1982;Albrecht & Steinhardt 1982;Hawking et al. 1982;Linde 1983) through the study of the statistics hidden in the scalar density perturbations (Mukhanov & Chibisov 1981).Euclid measurements are indeed expected to reduce uncertainties on the amplitude of primordial scalar perturbations A s , the scalar spectral index n s and its derivatives (or scalar runnings), the amplitude of primordial non-Gaussianity f NL (in particular of a local type), and the spatial curvature parameter Ω K ; and to enable us to perform significantly more stringent tests of extended models beyond singlefield slow-roll inflation (Laureijs et al. 2011;Amendola et al. 2018;Euclid Collaboration: Blanchard et al. 2020).
In the context of primordial features, large-volume photometric and spectroscopic surveys not only provide a complementary look at the structure on very large scales, but they can also provide a precise measurement of the power spectrum on small scales, complementing the CMB measurements and improving the sensitivity for high-frequency signals.This has been studied extensively with forecast analyses in Wang et al. (1999), Zhan et al. (2006), Huang et al. (2012), Chen et al. (2016a), Chen et al. (2016b), Ballardini et al. (2016), Xu et al. (2016), Ansari Fard & Baghram (2018), Palma et al. (2018), Ballardini et al. (2018), Ballardini (2019), Beutler et al. (2019), Ballardini et al. (2020), Debono et al. (2020), and Li et al. (2022) and demonstrated on real data in Beutler et al. (2019), Ballardini et al. (2023), and Mergulhão et al. (2023) showing that current galaxy clustering data of the Baryon Oscillation Spectroscopic Survey (BOSS) alone can already provide constraints that are competitive with those derived from the Planck CMB data.Given the potential of the Euclid mission, it is imperative to have a good description of the LSS observables on all observed scales, and to do so, nonlinear corrections need to be correctly accounted for (Vlah et al. 2016;Vasudevan et al. 2019;Beutler et al. 2019;Ballardini et al. 2020;Chen et al. 2020;Li et al. 2022;Ballardini & Finelli 2022).
In this paper, we forecast how well the surveys from the Euclid mission can constrain two templates of primordial undamped linear and logarithmic oscillations.In addition to the Euclid's primary probes, i.e. the spectroscopic galaxy clustering and the combination of photometric surveys, we study the further constraints that will be added by Euclid's measurements of the galaxy bispectrum and the information provided by future CMB experiments.The bispectrum is shown to be able to provide essential information on the parameters of the feature models, highlighting the importance of a high-order statistics analysis on the spectroscopic sample.
We structure the paper as follows.In Sect.2, we introduce the physics of the adiabatic mode and we briefly review the diverse theoretical mechanisms that naturally generate features in the PPS, in particular realisations of inflation.In Sect.3, we review Euclid specifications and calculation of the theoretical observables used in the Fisher forecast analysis.In Sect. 4 we start introducing the parametrised templates for primordial features with linear and logarithmic undamped oscillations.We then study the behaviour of these models on nonlinear scales of the matter power spectrum comparing the results from perturbation theory in Sect.4.2 to N-body simulations in Sect.4.3, and we study the numerical reconstruction of primordial features in Sect.4.4.We present our results from primary Euclid observables in Sect. 5. We complement those by adding the information from the galaxy bispectrum in Sect.6 and combining Euclid results with future CMB measurements in Sect.7. We conclude in Sect.8.

Primordial features and the early Universe
As mentioned in Sect. 1, the search for features in the primordial power spectrum has so far been statistically unsuccessful.On the largest scales, CMB observations agree with a red-tilted powerlaw PPS, P R, 0 (k), with the following parametrisation where A s and n s are the amplitude and the spectral index of the comoving curvature perturbations R on superhorizon scales at the pivot scale k * = 0.05 Mpc −1 .
The simplest class of models that can produce such a spectrum is known as canonical single-field slow-roll inflation.In these models a single scalar field -the inflaton -is responsible for both the primordial perturbations and the background evolution during inflation. 1 Unless the inflaton is the Higgs field (Bezrukov & Shaposhnikov 2008), inflation always takes us beyond the Standard Model (SM) of particle physics.While these scenarios can support cosmic inflation, it is rarely of the simplest kind.Beyond-SM models often involve multiple fields that interact with each other, and with the inflaton, and they can leave detectable imprints on the primordial perturbations including, for example, isocurvature modes and localised features.It is essential to have a bottom-up theoretical framework enabling us to analyse and classify all possible departures from the simplest scenario of Eq. ( 1) that are compatible with observations.Here we focus on the fluctuations as those are what we observe in the CMB and LSS.An important distinction is whether the curvature perturbations on superhorizon scales are generated by a single, possibly effective, degree of freedom.These are called effectively single field or 'single clock' models.
It can be shown that in the effective single-field slow-roll models the effective action for the comoving curvature perturbations is, to quadratic order, where a(t) is the scale factor, c s (t) is the speed of sound for curvature perturbations, and ϵ 1 (t) ≡ − Ḣ/H 2 is the first slow-roll parameter; H ≡ ȧ/a is the Hubble expansion rate, where the overdot denotes derivative with respect to the cosmic time t.Eq. ( 2) is well known in the context of single-field slow-roll inflation, but it is much more general than that.It is the first term in a perturbative expansion -the effective field theory (EFT) of inflationary perturbations (Cheung et al. 2008) -where all the information about the background is systematically encoded in a set of functions (ϵ 1 and c s being the first two) that describe what is seen by the perturbations.
This action describes almost all models of slow-roll inflation in which the primordial perturbations are generated by a single quantum field.It includes canonical single-field slow-roll models as a particular case (c s = 1) but also any multi-field model in which the primordial perturbations are generated by a single quantum field (typically an effective low energy degree of freedom involving multiple high energy fields).The action in Eq. ( 2) is then obtained by integrating out the fast, or heavy, high energy degrees of freedom; see Achúcarro et al. (2012).
The idea behind the EFT approach is that what the perturbations see is an almost time-independent background, and this symmetry under time translations completely dictates the form of the action in Eq. ( 2), to all orders.The red tilt is a small deviation from perfect scale invariance, measured by the smallness of the slow-roll parameters.Any deviations from timeindependence in the background functions result in features in the spectrum that can be calculated and cross-correlated among different observables; see Bartolo et al. (2013); Cannone et al. (2014) for applications of the EFT in the context of primordial oscillatory features.
In general, and depending on their origin, oscillatory features fall into two main classes which are the focus of this paper (see Achúcarro et al. 2022, for details and references).A small, abrupt change in the background functions ϵ 1 (t) and c s (t) at a particular moment during inflation results in a transient oscillatory feature, linearly spaced in k, superimposed on the power spectrum (1).Two well studied examples are a step in the inflationary potential (Starobinsky 1992;Adams et al. 2001) and coupling to gravity and the Bunch-Davies quantum vacuum for the perturbations.
a localised turn in the inflationary trajectory (Achúcarro et al. 2011;Chen 2012).The range of k where the oscillation persists increases with the sharpness of the change.But the sharpness cannot be arbitrarily large, since it should not excite the high energy modes that have been integrated out, nor enhance higher order terms that have been neglected in the perturbative expansion.On the other hand, oscillations logarithmically spaced in k are usually associated with periodicity in the background functions, for instance if there is a periodic modulation in the inflationary potential (Freese et al. 1990;Chen et al. 2008).More complicated combinations are also possible, that require a more tailored analysis.The upshot is that, if a feature is detected in the power spectrum, establishing its high energy origin will require stringent self-consistency checks, as well as correlated detections in other observables.
It is worth emphasising that the power of the EFT approach is that, since it is completely agnostic about the origin of the background expansion, it tests entire classes of models as opposed to individual ones.For example in multi-field models, ϵ 1 is related to the flatness of the potential along the inflationary trajectory, and c s is related to the rate of turning of this trajectory in multi-field space.But the high energy origin of ϵ 1 and c s may be completely different in other models.

Euclid probes
Euclid is one of the next generation deep-and large-field galaxy surveys.It will be able to measure up to 30 million spectroscopic redshifts, which can be used for galaxy clustering measurements, and 2 billion photometric galaxy images, which can be used for weak lensing observations.
Euclid will use near-IR (NIR) slitless spectroscopy to collect large samples of emission-line galaxies and to perform spectroscopic galaxy clustering measurements.The great advantage of this spectroscopic method is the high precision on measuring the redshift of the sources.However, one of the difficulties in inferring the cosmological parameters is given to the knowledge of the number density of the Hα targets.The total number of galaxies that can be used in the analysis mostly depends on two factors: completeness and purity.Completeness is the fraction of objects correctly identified relative to the true number, whereas purity is defined as the fraction of objects correctly identified relative to the total number of detected objects.In practice, completeness represents the size of a sample, whereas purity is its quality.A lower value of the number density n(z) of observed galaxies leads to an increase of the shot noise, resulting as a degradation on the constraints of the cosmological parameters.
One of the expected lensing observations is the cosmic shear, which is the distortion in the observed shapes of distant galaxies due to weak gravitational lensing by the LSS while the other probe comes from galaxy clustering using the positions of objects detected by the photometric measurements of redshifts.Given that both depend on the LSS density fluctuations as well as the geometry of the Universe, this will allow us to constrain the cosmological parameters.One of the main difficulties in using the weak lensing observable is the ability to accurately model the intrinsic alignments (IA) of galaxies that mimic the cosmological lensing signal.While for the galaxy clustering from photometric measurements, the main effects that need to be taken into account are the galaxy bias, accounting for the relation between the galaxy distribution and the underlying total matter distribution, and the photometric-redshift uncertainties, since redshift in this case is estimated from observing through multi-band filters instead of the full spectral energy distribution.Finally, an im-portant issue for both the weak lensing and galaxy clustering probes is the modelling of the small-scale nonlinear clustering on the weak lensing two-point statistics that the high density of detected galaxies will allow us to reach.
In the following we detail on the prescription used in this work for our model for these two probes.

Spectroscopic galaxy clustering
Following Euclid Collaboration: Blanchard et al. (2020), we define the observed galaxy power spectrum as F z (k, µ; z) where µ is the cosine of the angle of the wave mode with respect to the line of sight pointing into the direction r.Here b(z) is the linear clustering bias, f (z) is the growth rate, and σ 8 (z) is the root mean square linearly evolved density fluctuations in spheres of 8 h −1 Mpc. 2he Alcock-Paczynski (AP) effect (Alcock & Paczynski 1979), used to account for deviations of the cosmological models from the fiducial one, is parametrized through rescaling of the angular diameter distance D A (z) and the Hubble parameter H(z), and enters as multiplicative factor through where the subscript ref refers to the fiducial cosmology.This leads also to a rescaling of the wavevector components as Using Eqs. ( 4) and ( 5) we can convert the known reference cosmology (k ref and µ ref ) to the true unknown cosmology (k and µ), where The relations above can be used to define the effect of the choice of reference cosmology on the observed power spectrum via Galaxy bias, connecting the underlying dark matter (DM) power spectrum to the Hα-line emitter galaxies detected by Euclid, is modelled by a simple linear clustering bias redshiftdependent coefficient b(z).The anisotropic distortion to the density field due to the line-of-sight effects of the peculiar velocity of the observed galaxy redshifts, known as redshift-space distortions (RSD), is modelled in the linear regime by the contribution in the square brackets introduced by Kaiser (1987).Linear RSD are corrected for the nonlinear finger-of-God (FoG) effect under the assumption of an exponential galaxy velocity distribution function as a Lorentzian (Hamilton 1998).The nonlinearities due to gravitational instabilities have the effect of smearing the baryon acoustic oscillations (BAO) signal, this is modelled through time-sliced perturbation theory (TSPT); we discuss the modelling of P IR res, LO in Sect.4.2.The pairwise velocity dispersion, σ p (z), is evaluated from the linear matter power spectrum as The total galaxy power spectrum in Eq. ( 3) includes errors on the measurement of the redshift through the exponential factor where σ 2 r (z) = c(1 + z)σ 0,z /H(z) with σ 0,z the error on the measured redshifts.Finally, we introduce a shot-noise term P s (z) due to imperfect removal of the Poisson sampling and the imperfect modelling of small scales.
The final Fisher matrix for the spectroscopic galaxy clustering (GC sp ) observable for one redshift bin z i is where the derivatives are evaluated at the parameter values of the fiducial model and V eff is the effective volume of the survey, given by where V s is the volume of the survey and n(z) is the number of galaxies in a redshift bin.Assuming that the observed power spectrum follows a Gaussian distribution, we write in Eq. ( 12) its covariance matrix as where δ D (k − k ′ ) is the Dirac delta function.We model the observed matter power spectrum of Eq. ( 3) in bandpowers averaged over a bandwidth of ∆k with a top-hat window function as in Huang et al. (2012); Ballardini et al. (2016), The size of the k-bin ∆k should be equal or larger than the effective fundamental frequency defined as k eff ≡ 2π/V 1/3 for an ideal cubic volume in order to guarantee uncorrelated measurements of the observed power spectrum.The values of k eff span in the range 0.0025-0.0031h Mpc −1 for a cubic volume for the four redshift bins.We assume a bin width of ∆k = 0.004 h Mpc −1 with k min = 0.002 h Mpc −1 .Finally, we can rewrite the Fisher matrix of Eq. 12 for the discrete number of averaged bandpower for the redshift bin z i as The total spectroscopic Fisher matrix is then calculated by summing over the redshift bins as F sp αβ = i F αβ (z i ), assuming that the redshift bins are independent.Since we are interested in the standard cosmological parameters and the feature parameters, we marginalise the GC sp Fisher matrix over two redshiftdependent parameters bσ 8 (z i ) and P s (z i ) for each of the four redshift bins.

Photometric galaxy clustering and weak lensing
Both the forecasting method and the tools used for the photometric analyses are the same as the ones in Euclid Collaboration: Blanchard et al. (2020) apart from the changes in the power spectrum due to the primordial features in the predicted Euclid observables.Here we only remind the reader of the main steps.
The observables we consider are the angular power spectra C XY i j,ℓ between probe X in the i-th redshift bin and probe Y in the jth redshift bin, where the probes X and Y are L for weak lensing or G for photometric galaxy clustering; C XY i j,ℓ therefore refers to both auto-and cross-correlations of these probes.Relying on the Limber approximation and within the flat sky limit, the spectra are given by with k ℓ = (ℓ + 1/2)/r(z), r(z) the comoving distance to redshift z, and P NL (k ℓ , z) the nonlinear power spectrum of matter density fluctuations at wave number k ℓ and redshift z.The GC ph and WL window functions read where b i (z) is the galaxy bias in the i-th redshift bin, and W IA i (z) encodes the contribution of IA to the WL power spectrum.The normalised number density distribution n i (z) of observed galaxies in the i-th redshift bin is given by where (z − i , z + i ) are the edges of the i-th redshift bin and n(z) is the underlying true redshift distribution.The true number density of galaxies is convolved with the probability distribution function p ph (z ′ |z) is given following the parameterization in Euclid Collaboration: Blanchard et al. (2020).
The IA contribution is computed following the extended nonlinear alignment (eNLA) model adopted in Euclid Collaboration: Blanchard et al. (2020) so that the corresponding window function is where with ⟨L⟩(z) the redshift-dependent mean, and L ⋆ (z) the characteristic luminosity of source galaxies as computed from the luminosity function.A IA , β IA and η IA are the nuisance parameters of the model, and C IA is a constant accounting for dimensional units.
We consider a Gaussian-only covariance whose elements are given by where the upper-and lower-case Latin indices run over L and G (all tomographic bins), δ K ℓℓ ′ is the Kronecker delta coming from the lack of correlation between different multipoles (ℓ, ℓ ′ ), f sky is the survey's sky fraction, and ∆ℓ denotes the width of the logarithmic equi-spaced multipole bins.We consider a white noise, where σ ϵ is the variance of observed ellipticities.
For evaluating the Fisher matrix F ph αβ for the observed galaxy power spectrum, we use

Survey specifications
Following Euclid Collaboration: Blanchard et al. (2020), we consider for the spectroscopic sample four redshift bins centred at {1.0, 1.2, 1.4, 1.65}, whose widths are ∆z = 0.2 for the first three bins and ∆z = 0.3 for the last bin.The Hα bias, evaluated at the central redshift of the bins, corresponds to {1.46, 1.61, 1.75, 1.90}.The number density of galaxies n(z) corresponds to {6.86, 5.58, 4.21, 2.61} × 10 −4 .For the photometric probes, the sources are split into 10 equipopulated redshift bins whose limits are obtained from the redshift distribution with z 0 = 0.9/ √ 2 and the normalisation set by the requirement that the surface density of galaxies is ng = 30 arcmin −2 .This is then convolved with the sum of two Gaussians to account for the effect of photometric redshift (see Euclid Collaboration: Blanchard et al. 2020, for details).The galaxy bias is assumed to be constant within each redshift bin, with fiducial values b i = √ 1 + zi , where zi is the bin centre.We consider 100 logarithmic equi-spaced multipole bins with σ ϵ = 0.3 the variance of observed ellipticities.We set (z min , z max ) = (0.001, 4), which spans the full range where the source redshift distributions n i (z) are non-vanishing.
For all the Euclid probes the survey's sky fraction is f sky ≃ 0.36 (Euclid Collaboration: Scaramella et al. 2022).Power law LIN LOG Fig. 1.Primordial power spectrum of curvature perturbations for linear (LIN) and logarithmic (LOG) feature models with fiducial parameter values given in ( 29) and (30).

Theoretical modelling of primordial features in galaxy surveys
As representative models for features in primordial fluctuations in this paper we consider oscillations linearly and logarithmically spaced in Fourier space with a constant amplitude superimposed to a power spectrum described by a power-law.

Models
We consider a superimposed pattern of oscillations as where P R, 0 (k) is the standard power-law PPS of the comoving curvature perturbations R on superhorizon scales, given in Eq. ( 1).We study the following templates with superimposed oscillations on the PPS where X = {lin, log} and Ξ X = {k/k * , ln(k/k * )}.We choose the fiducial cosmological parameters according to Planck DR3 mean values of the marginalised posterior distributions (Planck Collaboration: Aghanim et al. 2020b): ω b = 0.02237, ω c = 0.1200, H 0 = 67.36km s −1 Mpc −1 , ln(10 10 A s ) = 3.044, n s = 0.9649, σ 8 = 0.8107, and m ν = 0.06 eV (with one massive neutrino and two massless ones).In what follows we specify the fiducial values of the model parameters.
1. Linear oscillations: 2. Logarithmic oscillations: We show the effect of the feature parameters on the PPS in Fig. 1.

Modelling nonlinear scales in perturbation theory in the presence of primordial oscillatory features
The linearly propagated matter power spectrum is given by where ), with T (k) the matter transfer function normalised to unity at large scales (i.e.k → 0) and c the speed of light.The gravitational potential power spectrum will then be derived from the two oscillatory models considered: P Φ = 9/25(2π 2 )P R (k)/k 3 .
Our analytic model for the matter power spectrum is based on TSPT and it takes into account the damping of oscillations by infrared (IR) resummation of the large-scale bulk flows (Crocce & Scoccimarro 2008;Creminelli et al. 2014;Baldauf et al. 2015;Blas et al. 2016b;Senatore & Trevisan 2018).The implementation of the IR resummation can be done following the approach to BAO in the context of TSPT (Blas et al. 2016a,b).We start by decomposing the linear matter power spectrum into a smooth (or no-wiggle; nw) and an oscillating (w) contribution, where P nw is the no-wiggle power spectrum, and the oscillatory part P w describes both the BAO feature and the primordial oscillating feature as Here we have factored out the time-dependence given by the growth factor D(z).We have neglected the cross term in Eq. ( 33) as it is proportional to A BAO × A lin and therefore subdominant.We filter the BAO feature from the linear matter power spectrum using a Savitzky-Golay filter; see Boyle & Komatsu (2018).
At next-to-leading order (NLO), the IR resummed power spectrum for linear oscillations can be written as (Blas et al. 2016b;Beutler et al. 2019) and for logarithmic oscillations as (Vasudevan et al. 2019;Beutler et al. 2019) P 1−loop is the standard one-loop result, but computed with the leading-order (LO) IR resummed power spectrum.We take the usual expression P 1−loop = P 22 + 2P 13 with while evaluating the loop integrals P 22 and P 13 with the input spectrum P IR res, LO instead of the linear spectrum (Blas et al. 2016a).F n are the usual perturbation theory (PT) kernels (Bernardeau et al. 2002).We calculate these quantities, i.e.Eqs. ( 36) and ( 37), with the publicly available code FAST-PT (McEwen et al. 2016;Fang et al. 2017). 3The IR resummed power spectrum at leading order is given by for the linear oscillations (Blas et al. 2016b;Beutler et al. 2019) and by for the logarithmic oscillations (Vasudevan et al. 2019;Beutler et al. 2019).The result of the IR resummation at LO is given by a first contribution corresponding to the smooth part of the linear power spectrum and a second contribution corrected by the exponential damping of the oscillatory part due to the effect of IR enhanced loop contributions.The result for the logarithmic oscillations has additional contributions, on top of scale-dependent damping factors, due to the non-trivial oscillatory behaviour in real space.The damping factors above correspond to where j n are spherical Bessel functions, r s ≃ 147 Mpc is the scale setting the period of the BAO (Planck Collaboration: Aghanim et al. 2020b), and k S is the separation scale controlling the modes which are to be resummed.The dependence on k S can be connected with an estimate of the theoretical perturbative uncertainties.For this reason and since IR expansions are valid for q ≪ k, we assume k S = ϵk with ϵ ∈ [0.3, 0.7] (Baldauf et al. 2015;Vasudevan et al. 2019).
In redshift space, the damping factor becomes also dependent on the cosine µ since peculiar velocities additionally wash out wiggle signals along the line of sight.In this case, we have at leading order (Eisenstein et al. 2007) in Eqs. ( 38) and ( 39) 3 https://github.com/JoeMcEwen/FAST-PTWe note that while in the next subsection, Sect.4.3, we present both PT results at LO and NLO, that is Eqs.( 38)-( 39)-( 34)-( 35), and we compare them to N-body simulations.Furthermore, we only used the LO to derive the Fisher-matrix uncertainties in Sect. 5.
For the nonlinear modelling for photometric probes, entering Eq. ( 17), we use a modified version of HMCODE (Mead et al. 2016) where the ΛCDM nonlinear matter power spectrum is dressed with the primordial oscillations according to

Comparing predictions of perturbation theory with N-body simulations
We produce cosmological simulations based on the COmoving Lagrangian Approximation (COLA) method (Tassev et al. 2013(Tassev et al. , 2015;;Winther et al. 2017;Wright et al. 2017) in order to assess the accuracy of the predictions of PT at LO and NLO for the template (28) using a modified version of the publicly available code L-PICOLA (Howlett et al. 2015). 4We fix the standard cosmological parameters and the amplitude and phase of the feature, and we study the effect of varying the frequency of the primordial oscillations over the set of values ω X ∈ {0.2, 0.4, 0.6, 0.8, 1, 1.2, 1.4, 1.6, 1.8, 2}.Following the analysis done by Ballardini & Finelli (2022), we run each simulation with 1024 3 dark matter particles in a comoving box with side length of 1024 h −1 Mpc evolved with 30 time steps.We set at redshift z = 9 the initial conditions that we generate using second-order Lagrangian perturbation theory, with the 2LPTic code (Crocce et al. 2006).Finally, we use spectra averaged over pairs of simulations with the same initial seeds and inverted initial conditions, and with amplitude fixing in order to minimise the cosmic variance (Viel et al. 2010;Villaescusa-Navarro et al. 2018).In Fig. 2, we show a comparison between the nonlinear matter power spectrum at z = 0 simulated with COLA and the one obtained by Ballardini et al. (2020) at the same resolution with the N-body code GADGET-3, a modified version of the publicly available code GADGET-2 (Springel et al. 2001;Springel 2005), 5 for A lin = 0.03 and ω lin = 10.).We also show the results of the linear theory (grey) and the ones obtained from N-body simulations (green).
We show in Fig. 3 the ratio between the matter power spectrum with primordial oscillations and the one with featureless PPS calculated at redshift z = 0 for the results both at LO and NLO for the two templates with linear and logarithmic oscillations.We collect in Appendix A the results for different frequencies; see Figs.A.1 and A.2.
For the template with undamped linear oscillations, the agreement between N-body simulations and PT results is remarkable for the entire range of frequencies considered, i.e. log 10 ω lin ∈ (0.2, 2).At LO, we find differences for the matter power spectrum in Fourier space at redshift z = 0 less than 1% for log 10 ω lin > 0.6 and between 1% and 2% for log 10 ω lin ≤ 0.6.At NLO, we find an agreement better than 1% for all the frequencies and differences < 0.5% for log 10 ω lin > 0.4.For the template with undamped logarithmic oscillations, we find a similar agreement for the frequencies log 10 ω log ≥ 1, less than 1% for the LO and less than 0.5% for the NLO.In order to improve the accuracy of the results for low frequencies (see Beutler et al. 2019, for a discussion on the validity of perturbation approach for small frequencies) we can fit a Gaussian damping to the Nbody simulation as done by Ballardini et al. (2020).
While we restrict the validation of results to the matter power spectrum in real space, a validation for the halo power spectrum in redshift space has been performed by Chen et al. (2020), showing that current results from PT are able to describe redshift-space clustering observables also in the presence of primordial features.

Nonlinear reconstruction of primordial features
An alternative to the modelling approach described above is through reconstruction.This takes the opposite logic, by starting from an evolved, nonlinear, matter field where the primor-dial oscillations have been damped, it tries to recover their undamped states.Reconstruction has long been adopted as a useful technique to undo the effect of structure evolution and retrieve sharper BAO features, hence improving their use as a standard ruler to constrain cosmological parameters (Eisenstein et al. 2007;Wang et al. 2017;Sarpa et al. 2019), with an even longer history of applications in other subfields of cosmology and astrophysics (Peebles 1989;Croft & Gaztanaga 1997;Monaco & Efstathiou 1999;Nusser & Branchini 2000;Brenier et al. 2003;Mohayaee et al. 2003;Schmittfull et al. 2017;Zhu et al. 2017;Shi et al. 2018;Birkin et al. 2019;Mao et al. 2021).
From the point of view of reconstruction, there is no fundamental difference between BAO features and features that exist in the primordial density fluctuations on similar length scales: both are present at the time of last scattering, the preservation of both is negatively impacted by late-time cosmic structure formation, and both can be at least partially recovered if such an impact is undone through processes such as reconstruction (Beutler et al. 2019;Li et al. 2022).Therefore, we expect that reconstruction can find uses in the study of primordial features, at least when the oscillations have similar frequencies as the BAO peaks, 0.05 ≲ k/(h Mpc −1 ) ≲ 0.5 -at scales much larger than this range any significant damping of the features due to cosmological evolution is yet to happen, because structure formation is hierarchical and affects smaller scales first, while on much smaller scales reconstruction becomes unreliable.
Another useful property of reconstruction is that it is possible to include galaxy bias (Birkin et al. 2019) and RSD (Monaco & Efstathiou 1999;Burden et al. 2014;Hada & Eisenstein 2018, 2019;Wang et al. 2020) in the pipeline, so that two of the most important theoretical systematic effects in analysing the galaxy redshift power spectrum are taken into account.Of course, caution is due here to ensure that the range of length scales where reconstruction from a catalogue of redshift-space biased tracers can be done reliably covers the range where we hope reconstruction to also benefit the recovery of the primordial features of interest.Checking this naturally requires a more detailed analysis, which may also depend on the reconstruction algorithm in actual use: while a quantitative analysis of this issue is not the topic here, our experience is that the scale range k ≲ 0.2 h Mpc −1 or s ≳ 15-20 h −1 Mpc satisfies these requirements.
In this work, we exemplify using the nonlinear reconstruction method described in Shi et al. (2018); Birkin et al. (2019); Wang et al. (2020).Other algorithms may give quantitatively different results, though we do not expect a huge variation among the ones that are at an advanced stage of development.The problem of reconstruction can be reduced to identifying a mapping between the initial, Lagrangian coordinate of some particle, q, and its Eulerian coordinate after evolution, x(t), at a later time t.If trajectory crossings of particles have not happened, this mapping is unique and can be obtained by solving the mass conservation equation ρ(x) where ρ(q) and ρ(x) are, respectively, the initial density field in some infinitesimal volume element dq, and the density field in the same volume element at time t, which now has been moved, and possibly deformed, and described by dx.As the density field in the early Universe is homogeneous to a very good approximation, one can take ρ(q) to be a constant, i.e. ρ(q) ≃ ρ.We define the displacement field Ψ between the final and initial positions of a particle as Ψ(x) = x − q, which can be rewritten as with Θ(x) being the displacement potential.Here we have assumed the displacement field to be curl-free, i.e. ∇ × Ψ = 0, which is valid only on large scales. 6Combining Eq. ( 47) and Eq. ( 46) gives where δ(x) = ρ(x)/ ρ − 1 is the density contrast at time t, 'det' is the determinant of a matrix, here the Hessian of Θ(x), and i, j = 1, 2, 3. Shi et al. (2018) proposed to solve this equation as a nonlinear partial differential equation with cubic power of second-order derivatives of Θ.This was later extended by Birkin et al. (2019) to cases where δ(x) in the above can be the density contrast of some biased tracers of the dark matter field, such as galaxies and dark matter haloes.
Once Θ(x) is solved, one can obtain Ψ(x) and the reconstructed density field is given by where we express Ψ in terms of the Lagrangian coordinate q, since the divergence ∇ q therein is with respect to q.This calculation means that we need to have Ψ(q) on a regular q-grid, which can be done using, for example the Delaunay  2011), by interpolating Ψ(x) to the target q-grid.
The quantity δ r obtained above is an approximation to the initial linear matter density contrast, linearly extrapolated to time t (bearing in mind that all the discussion in this subsection is for ΛCDM, for which the linear growth factor is scale-independent).As a result, part of the structure-formation-induced damping of the BAO or primordial features imprinted in the density field can be undone this way.
The damping of the BAO or primordial wiggles can be described by a Gaussian fitting function (cf.Sect. 4.3;Vasudevan et al. 2019;Beutler et al. 2019;Ballardini et al. 2020), effectively multiplying the undamped power spectrum wiggles O undamped w by a Gaussian function where ζ(z) is the parameter quantifying the level of damping.
In the linearly evolved density field with no damping, one has ζ = 0.For the evolved nonlinear matter or galaxy field, ζ > 0 with its values depending on redshift, tracer type, tracer number density etc., and the reconstructed density field generally features a smaller yet positive value of ζ.
We generate a suite of N-body simulations using the initial conditions described in Sect.4.3, which includes a no-wiggle model and 20 wiggled models.The simulations are performed using the parallel N-body code ramses (Teyssier 2002).We output 4 snapshots of the DM respectively at z = 1.5, 1.0, 0.5 and 0, and measure the matter power spectrum of each of the snapshots.We follow the same reconstruction pipeline employed in Li et al. (2022) to reconstruct the density fields from each of the snapshots for all the models (shown in Fig. B.1).We have (k, z).We apply the Hilbert transform to measure the wiggle envelopes for each wiggled model and each redshift.Nevertheless, we expect that the damping parameter ζ(z), for both the unreconstructed and the reconstructed cases, should mostly depend on the nonlinear structure formation and be insensitive to the wiggle model, given that the wiggles are weak.We have explicitly checked this by comparing the envelopes of O damped w and O undamped w of the different models, and finding them to agree very well.As a result, to get ζ(z), the O w envelopes of all wiggle models at redshift z could be combined to do a single least-squares fitting.In practice, because the measurement of the envelopes is not very reliable for the lowfrequency featured models due to the very few wiggles in the k range of fitting, k = (0.05-1) h Mpc −1 , we select only the 10 high-frequency feature models to fit ζ(z), and the result is given in Table 1.We nevertheless have checked that the fitted envelope also well describes the wiggles of low-frequency models.

Expected constraints from Euclid primary probes
In this section, we show the results of the Fisher analysis for the cosmological parameters of interest corresponding to a flat ΛCDM, after marginalisation over spectroscopic and photometric nuisance parameters.We compute the spectroscopic galaxy clustering (GC sp ), the photometric galaxy clustering (GC ph ), the weak lensing (WL), and the cross-correlation between the two photometric probes (XC) in two configurations: a pessimistic setting and an optimistic one.For the pessimistic setting, for GC sp we have k max = 0.25 h Mpc −1 , for GC ph and XC we have ℓ max = 750, and for WL we have ℓ max = 1500.In addition, we impose a cut in redshift of z < 0.9 to GC ph to limit any possible cross-correlation between spectroscopic and photometric data.For the optimistic setting, we extend GC sp to k max = 0.30 h Mpc −1 , GC ph and XC to ℓ max = 3000, and WL to ℓ max = 5000, without imposing any redshift cut to GC ph .
We start by presenting results for a fiducial scenario with A X = 0.01, ω X = 10, and ϕ X = 0.The marginalised uncertainties on the cosmological parameters and primordial feature parameters, percentages relative to the corresponding fiducial values, are collected in Table 2 for the pessimistic (top panel) and optimistic (bottom panel) settings.The marginalised 68.3% and 95.5% confidence level (CL) contours for the primordial feature parameters are shown in Fig. 4. For the fiducial value parameters, uncertainties on the amplitude of the primordial feature oscillations are dominated by GC sp measurements resulting in A lin = 0.0100 ± 0.0023 (±0.0023) and A log = 0.0100 ± 0.0025 (±0.0025) at a 68.3% CL for the pessimistic (optimistic) setting.Combining GC sp with the combination of photometric information, i.e.WL+GC ph +XC, we find A lin = 0.0100 ± 0.0021 (±0.0018) and A log = 0.0100 ± 0.0022 (±0.0018) at a 68.3% CL for the pessimistic (optimistic) setting.
For GC sp , bounds on the primordial amplitude parameter A X strongly depend on the primordial frequency value ω X ; see Huang et al. (2012); Ballardini et al. (2016); Slosar et al. (2019); Beutler et al. (2019); Ballardini et al. (2020).Low frequencies, lower than the BAO one, leave a broad modification on the matter power spectrum and are expected to be better constrained by CMB measurements.On the other hand, the imprint from highfrequency primordial oscillations is smoothed by projection effects on the observed matter power spectrum.Moreover, for the linear model the uncertainties for frequencies around the BAO frequency, i.e. log 10 ω lin ∼ 0.87, are degraded.For the combination of photometric probes, bounds on the primordial amplitude parameter A X are smoothed by projection effects in angular space and high-frequency primordial oscillations are severely washed out on the observed matter power spectrum.On the other hand, uncertainties for frequencies around the BAO scale and lower frequencies are tighter compared to the results from GC sp .
This represents an important result considering the lower accuracy of the PT predictions for low frequencies.
In Fig. 5, we show the uncertainties at a 68.3% CL for different values of the primordial frequency log 10 ω X within the range (0.1, 2.1).The combination of the spectroscopic and photometric probes allows us to reach uncertainties of 0.002-0.003on the primordial feature amplitude for the entire range of frequencies.

Power spectrum and bispectrum combination
In this section, we present the forecast results on primordial features by combining the power spectrum and bispectrum signals.
The presence of a sharp feature in the inflaton potential violates the slow-roll evolution, generating a linearly spaced oscillating primordial bispectrum.In this case the primordial bispectrum can be characterised by the standard amplitude parameter f lin NL , the phase ϕ B lin , and the frequency ω lin , and it is given by (Chen et al. 2007(Chen et al. , 2008) where A s is the normalisation parameter of P Φ = A/k (4−n s ) .A periodically modulated inflaton potential leads to resonances in the inflationary fluctuations with logarithmicallyspaced oscillations (Chen et al. 2008).This generates oscillatory features in the primordial power spectrum (see Sect. 4) and bispectrum (Flauger et al. 2010;Hannestad et al. 2010;Barnaby et al. 2012).In this case, the primordial bispectrum is given by (Chen 2010) where ϕ B log and f log NL are the phase and the amplitude, respectively, of the generated logarithmic features on the primordial bispectrum.Note that here we assume the frequency of the primordial feature bispectrum to be the same as in the case of the power spectrum (see Sect. 4.1) The redshift space model, for the galaxy bispectrum, is given after considering non-Gaussian initial conditions, as dictated by the oscillatory features, up to second order terms in RSD, bias, and matter expansions.The tree-level modelling can be written as where the linearly propagated primordial oscillatory bispectrum is Φ given by Eqs. ( 52) and ( 53) for the linear and logarithmic oscillations, respectively.The redshift space kernels are given by where µ i = k i • ẑ/k i is the cosine of the angle between the wave vector k i and the line of sight, µ = µ 1 + . . .µ n , and The kernels F 2 (k i , k j ) and G 2 (k i , k j ) are the second-order symmetric PT kernels (Bernardeau et al. 2002) − 1/3 is the tidal kernel (McDonald & Roy 2009;Baldauf et al. 2012).The second-order tidal field bias term, following the convention of (Baldauf et al. 2012), is given by b s 2 (z) = −4[b 1 (z) − 1]/7.The redshift-space bispectrum is characterised by five variables: three to define the triangle shape (the sides k 1 , k 2 , k 3 ) and two that characterise the orientation of the triangle relative to the line of sight, i.e.
The forecasts presented here use the information of the bispectrum monopole obtained after taking the average over all angles.
Note that the presence of an oscillatory primordial bispectrum generates a scale-dependent correction to the linear bias, in a similar manner to the local primordial non-Gaussian case (see Karagiannis et al. 2018, for discussions).Cabass et al. (2018) show that the scale-dependent correction oscillates with scale with an envelope similar to that of equilateral primordial non-Gaussianity (PNG) (see Schmidt & Kamionkowski 2010;Scoccimarro et al. 2012;Schmid & Hui 2013;Assassi et al. 2015, for a discussion), while their findings indicate that the amplitude of such a term is very small to be detected by upcoming surveys.Therefore, we do not consider the scale-dependent terms in our analysis and exclude them from the expressions of the redshift kernels presented above.This means that the degeneracy between the scale-dependent bias coefficient and the primordial bispectrum amplitude ( f X NL ), shown in Barreira (2020), does not affect our forecasts.
The multiplicative factors D P (k, z) and D B (k 1 , k 2 , k 3 , z), which incorporate errors on redshift measurements and FOG effect, are D P (k, z) = exp[−kµ(σ 2 p (z) + σ 2 r (z))] and 3 )(σ 2 p (z)+σ 2 r (z))/2] for the power spectrum and bispectrum, respectively (see Sect. 3.1 for details).The fiducial values of the stochastic terms in Eqs. ( 54) and ( 55) are taken to be those predicted by Poisson statistics (Schmidt 2016;Desjacques et al. 2018), i.e.P ε = 1/n g , P εε δ = b 1 /(2n g ), and B ε = 1/n 2 g .The AP effect is taken into account also for the bispectrum similar to the power spectrum case (see Sect. 3.1).The observed bispectrum is thus given by (see Karagiannis et al. 2022) Here we again use the Fisher matrix formalism to produce forecasts on the model parameters that control the primordial oscillations.The Fisher matrix of the redshift galaxy power spectrum is given by Eq. ( 12), while for the bispectrum the Fisher matrix for one redshift bin z i is given by where the sum over triangles has and k 1 , k 2 and k 3 satisfy the triangle inequality.The bin size ∆k is taken to be the fundamental frequency of the survey, k f = 2π/L, where for simplicity we approximate the survey volume as a cube, L = V 1/3 survey .The k mode is binned with a bin size of ∆k, and here we consider ∆k = k f , between its minimum and maximum values k min = k f and k max .Gaussian approximation is used for the covariance of bispectrum, i.e. the off-diagonal terms are considered to be zero. 7The variance for the bispectrum is then (Sefusatti et al. 2007) where s 123 = 6, 2, 1 for equilateral, isosceles, and non-isosceles triangles, respectively.In addition, for degenerate configurations, i.e. k i = k j + k m , the bispectrum variance should be multiplied by a factor of 2 (Chan & Blot 2017;Desjacques et al. 2018).We use the Fisher matrix of Eq. ( 60) to generate bispectrum forecasts on the parameters of interest The initial parameter vector considers all the parameters of the bispectrum tree-level modelling (i.e.cosmological parameters, biases, AP parameters, growth rate, and FOG amplitudes) to be free, and after marginalisation we derive the constraints for Θ B final .The forecasted marginalised 1σ errors on the feature model parameters, coming from the spectroscopic galaxy clustering, are presented in Table 3.
The galaxy bispectrum is capable of providing constraints not only on the feature parameters related to the primordial bispectrum but also on those controlling the primordial power spectrum feature model.The latter is due to the P lin (k i )P lin (k j ) terms that are present at tree level in the galaxy bispectrum (55) and can provide constraints equivalent to or even better than the galaxy power spectrum, especially in the case of a high k max .This can be seen in Table 3, where adding the power spectrum to the bispectrum improves the forecasts by a few percent.This highlights the importance of a bispectrum analysis for the type of feature models discussed in this work.On the other hand, the sole contribution to the constraints on the f X NL and ϕ B X feature parameters is the primordial component B I of the galaxy bispectrum, justifying the less rigorous forecasts presented here.

Combination with CMB data
Future CMB polarisation data from ground-based experiments, such as Simons Observatory (Ade et al. 2019) and CMB-S4 (Abazajian et al. 2019), and satellites, such as LiteBIRD (Lite-BIRD Collaboration: Allys et al. 2023), will be able to reduce the uncertainties on the PPS, in particular in the presence of oscillatory signals.Indeed, E-mode polarisation is sourced by velocity gradient (scattering only) leading to sharper transfer functions compared to the temperature ones helping to characterise primordial features in the PPS; see Mortonson et al. (2009)  Notes.We show results for LIN and LOG models in the pessimistic and optimistic settings in the case, using the Euclid GC sp bispectrum and the joint power spectrum+bispectrum signal.
We consider the information available in the CMB data by using power spectra from temperature, E-mode polarisation, and CMB lensing in combination with Euclid observations.We do not consider B-mode polarisation since it does not add information for the specific models examined here.We also neglect the cross-spectra between CMB fields and Euclid primary probes; see Euclid Collaboration: Ilić et al. ( 2022) for a complete characterisation and quantification of the importance of crosscorrelation between CMB and Euclid measurements.Following Euclid Collaboration: Ilić et al. (2022), we consider three CMB experiments: Planck-like, Simons Observatory (SO), and CMB-S4.Noise curves correspond to isotropic noise deconvolved with instrumental beam (Knox 1995) where θ FWHM is the full-width-at-half-maximum (FWHM) of the beam, w T and w E are the inverse square of the detector noise levels for temperature and polarisation, respectively.CMB lensing noise is reconstructed through minimum-variance estimator for the lensing potential using quicklens code;8 see Okamoto & Hu (2003).
For the Planck-like experiment, we use an effective sensitivity in order to reproduce the Planck 2018 results for the ΛCDM model (Planck Collaboration: Aghanim et al. 2020b) avoiding the complexity of the real-data likelihood (Planck Collaboration: Aghanim et al. 2020c).We use noise specifications corresponding to in-flight performances of the 143 GHz channel of the high-frequency instrument (HFI) (Planck Collaboration: Aghanim et al. 2020a) with a sky fraction f sky = 0.7 and a multipole range for temperature and polarisation from ℓ min = 2 to ℓ max = 1500.The E-mode polarisation noise is inflated by a factor of 8 to reproduce the uncertainty of the optical depth parameter τ; see Bermejo-Climent et al. (2021).Finally, CMB lensing is obtained by combining the 143 GHz and 217 GHz HFI channels assuming a conservative multipole range of 8-400.For the CMB, we also consider the optical depth at reionisation τ and then we marginalise over it before combining the CMB Fisher matrix with the Euclid ones.
For the SO-like experiment, we use noise curves provided by the SO collaboration in Ade et al. (2019) taking into account residuals and noises from component separation. 9We use spectra from ℓ min = 40 to ℓ max = 3000 for temperature and temperature-polarisation cross-correlation, and ℓ max = 5000 for E-mode polarisation.CMB lensing and temperature-lensing cross-correlation spectra cover the multipole range of 2-3000.The sky fraction considered is f sky = 0.4.We complement the SO-like information with the Planck-like large-scale data in the multipole range of 2-40 in both temperature and polarisation.
For CMB-S4, we use noise sensitivities of 1 µK arcmin in temperature and √ 2 µK arcmin in polarisation, with resolution of θ FWHM = 1 arcmin.We assume data over the same multipole ranges of SO and with the same sky coverage.Notes.We show results for LIN and LOG models in the pessimistic and optimistic settings in the case, relative to their corresponding fiducial values (A X = 0.01, ω X = 10, using Euclid (GC sp +WL+GC ph +XC) in combination with the CMB.
In Table 4, we show the marginalised uncertainties on the primordial feature parameters, percentage relative to the corresponding fiducial values, for the optimistic and pessimistic Eu-clid + CMB combination.Uncertainties on the primordial feature amplitude improve by 15%-50% depending on the choice of Euclid settings and on the CMB experiment considered, for a frequency value of ω X = 10.

Conclusions
We have explored the constraining power of the future Euclid space mission for oscillatory primordial features in the primordial power spectrum (PPS) of density perturbations.We have considered two templates with undamped oscillations, i.e. with constant amplitude all over the PPS, one linearly and one logarithmically spaced in k space.We have used a common baseline set of parameters for both the templates with amplitude A X = 0.01, frequency ω X = 10, and phase ϕ X = 0, where X = {lin, log}.
Following previous studies for ΛCDM and simple extensions in Euclid Collaboration: Blanchard et al. (2020), we have calculated Fisher-matrix-based uncertainties from Euclid primary probes, i.e. spectroscopic galaxy clustering (GC sp ) and the combination of photometric galaxy clustering (GC ph ) and photometric weak lensing (WL); see Sect. 3. We have considered two sets of Euclid specifications: a pessimistic setting with k max = 0.25 h Mpc −1 for GC sp , ℓ max = 1500 for WL, ℓ max = 750 for GC ph and XC, and a cut in redshift of z < 0.9 to GC ph ; an optimistic setting with k max = 0.30 h Mpc −1 for GC sp , ℓ max = 5000 for WL, and ℓ max = 3000 for GC ph and XC.
We have modelled the nonlinear matter power spectrum with time-sliced perturbation theory predictions at leading order (Vasudevan et al. 2019;Beutler et al. 2019;Ballardini et al. 2020;Ballardini & Finelli 2022) calibrated on a set of N-body simulations from COLA (Tassev et al. 2013;Howlett et al. 2015) and GADGET (Springel et al. 2001;Springel 2005;Ballardini et al. 2020;Ballardini & Finelli 2022) for GC sp ; see Sect.4.2.For modelling of the photometric probes we rely on HMCODE (Mead et al. 2016) to describe the broadband small-scale behaviour and on PT to describe the smearing of BAO and primordial features; see Sect.4.2.
With the full set of probes, i.e.GC sp +WL+GC ph +XC, we have found that Euclid alone is able to constrain the amplitude of a primordial oscillatory feature with frequency ω X = 10 to A X = 0.010 ± 0.002 at a 68.3% CL both in the pessimistic and optimistic settings; see Sect. 5 and Table 2.When we consider single Euclid probes these uncertainties depend strongly on the frequency (see Fig. 5) and are weakened for low and high frequencies.However, we have found robust constraints of 0.002-0.003from the full combination of Euclid probes over the frequency range of (1, 10 2.1 ).
We have then considered further information available from Euclid measurements including a numerical reconstruction of nonlinear spectroscopic galaxy clustering (Beutler et al. 2019;Li et al. 2022) and the galaxy clustering bispectrum (Karagiannis et al. 2018); see Sect.4.4 and Sect.6, respectively.We have found A lin = 0.0100 ± 0.0014 (A log = 0.0100 ± 0.0015) in the pessimistic setting and A lin = 0.0100 ± 0.0012 (A log = 0.0100 ± 0.0012) in the optimistic setting, both at 68.3% CL, from GC sp (rec)+WL+GC ph +XC.This shows that reconstruction can potentially tighten the constraints on the primordial feature amplitude.Including the galaxy clustering bispectrum, we further tighten the uncertainties on the primordial feature amplitude to A lin = 0.0100 ± 0.0009 (A log = 0.0100 ± 0.0009) in the pessimistic setting and A lin = 0.0100 ± 0.0008 (A log = 0.0100 ± 0.0008) in the optimistic setting, both at 68.3% CL, from GC sp (rec)+WL+GC ph +XC+BS.
Finally, we have studied the combination of Euclid probes and the expected information from CMB experiments by adding (without including the cross-correlation) forecasted results from Planck, Simons Observatory (SO), and CMB-S4 following Euclid Collaboration: Ilić et al. (2022); see Sect. 7. Combining the Fisher matrix information from a SO-like experiment (complemented with a Planck-like one at low-ℓ) with the Euclid primary probes, we have found A lin = 0.0100 ± 0.0014 (A log = 0.0100 ± 0.0011) in the pessimistic setting and A lin = 0.0100 ± 0.0013 (A log = 0.0100 ± 0.0010) in the optimistic setting, both at 68.3% CL.
We summarise in Fig. 6 the comparison of Euclid measurements in combination with other non-primary sources of information.Our tightest uncertainties, combining all the sources of information expected from Euclid in combination with future CMB experiments, e.g SO complemented with Planck at low-ℓ, correspond to A lin = 0.0100 ± 0.0008 at a 68.3% CL and to A log = 0.0100 ± 0.0008 for GC sp (PS rec + BS)+WL+GC ph +XC+SO-like for both the pessimistic and the optimistic settings.
Our results highlight the power of Euclid measurements to constrain primordial feature signals helping us to reach a more complete picture of the physics of the early Universe and to improve over the current bounds on the primordial feature amplitude parameter A X .In addition, the expected observations from Euclid will allow us to scrutinise the primordial interpretation of some of the anomalies in the CMB temperature and polarisation angular power spectra corresponding to the following best-fits of Planck data (Planck Collaboration: Aghanim et al.

Fig. 2 .
Fig.2.Ratio with respect to the ΛCDM case of the nonlinear matter power spectra in Fourier space at redshift z = 0 calculated with the approximate N-body method COLA and with the full N-body code GADGET-3 for the linear oscillations with A lin = 0.03 and ω lin= 10  (Ballardini et al. 2020).

Fig. 3 .
Fig.3.Ratio of IR resummed matter power spectrum at LO (blue) and NLO (orange) obtained for linear (top panel) and logarithmic (bottom panel) oscillations to the one obtained with a power-law PPS at redshift z = 0 when the frequency ω X = 10 and the IR separation scale k S = ϵk is varied (with ϵ ∈ [0.3, 0.7]).We also show the results of the linear theory (grey) and the ones obtained from N-body simulations (green).
Fig.4.Fisher-forecast marginalised two-dimensional contours and one-dimensional probability distribution functions from Euclid on the primordial feature parameters for the LIN model with ω lin = 10 (top panels) and the LOG model with ω log = 10 (bottom panels).Left (right) panels correspond to the pessimistic (optimistic) setting for GC sp (green), WL+GC ph +XC (orange), and their combination (blue).

Fig. 5 .
Fig.5.Marginalised uncertainties on A X as a function of the primordial frequency ω X for the LIN model (top panel) and the LOG model (bottom panel).We show uncertainties for GC sp (blue), WL+GC ph +XC (orange), and their combination (green).Dashed (solid) curves correspond to the pessimistic (optimistic) setting.

Fig. 6 .
Fig.6.Fisher-forecast marginalised two-dimensional contours and one-dimensional probability distribution functions from Euclid on the primordial feature parameters for the LIN model with ω lin = 10 (top panels) and the LOG model with ω log = 10 (bottom panels).Left (right) panels correspond to the pessimistic (optimistic) setting for GC sp +WL+GC ph +XC (red), GC sp +WL+GC ph +XC with numerical reconstruction of GC sp (green), GC sp (rec)+WL+GC ph +XC in combination with Euclid GC sp bispectrum (orange), and their combination plus SO-like CMB (blue).

Fig. B. 1 .
Fig. B.1.Ratio of the linear (black solid curves), nonlinear (blue solid curves), and reconstructed (orange dashed curves) density fields for the 20 selected feature models to the one obtained with a power-law PPS.Each of the rows is respectively for the redshifts marked on the right.The linear O w (k) are measured from the initial conditions at z = 9, the nonlinear O w (k) are from the output snapshots of DM, and the reconstructed O w (k) are from the reconstructed density fields.

Table 1 .
Best-fit values of ζ(z)for unreconstructed (reconstructed) wiggle spectra.performedreconstructionfromhalo catalogues, though the results are noisier and not shown here; the DM reconstruction result can be considered as an ideal scenario, which serves as a limit of what can be achieved from halo or galaxy reconstruction with increasing tracer number density.We present the reconstruction from mock galaxy catalogues in redshift space based on the same set of simulations in a separate work.Fig.B.1 shows that reconstruction can help recover the weakened wiggles down to k ≃ 1 h Mpc −1 .Note that, in order to avoid artificial features in the P(k) of the models with high-frequency wiggles, we have binned P(k) in k bins of width 0.001 h Mpc −1 .To estimate the best-fit values of ζ(z), the Gaussian function is applied to fit the envelopes of the unreconstructed and reconstructed wiggled spectra, O also

Table 2 .
Fisher-forecast 68.3% CL marginalised uncertainties on cosmological and primordial feature parameters, relative to their corresponding fiducial values.Notes.We show results for LIN and LOG models in the pessimistic and optimistic settings, using Euclid observations GC sp , WL+GC ph +XC, and their combination.

Table 3 .
Fisher-forecast 68.3% CL marginalised uncertainties on primordial feature parameters, relative to their corresponding fiducial values, using Euclid GC sp bispectrum.

Table 4 .
Fisher-matrix-based forecasted 68.3% CL marginalised uncertainties using Euclid in combination with the CMB.