© The Authors 2024

1. Introduction

Future galaxy surveys are the new frontier in the study of initial conditions of the Universe. They are expected to improve significantly 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, 1983; Albrecht & Steinhardt 1982; Hawking et al. 1982) 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 single-field slow-roll inflation (Laureijs et al. 2011; Amendola et al. 2018; Euclid Collaboration 2020).

The goal of this paper is to summarise the status of, motivations for, and challenges in searching for oscillatory features in the primordial fluctuations (see Chluba et al. 2015; Achúcarro et al. 2022, for reviews) with large-scale structure (LSS) observations expected from Euclid. Searches for primordial features based on the cosmic microwave background (CMB) angular power spectra (Wang & Mathews 2002; Adams et al. 2001; Peiris et al. 2003; Mukherjee & Wang 2003; Covi et al. 2006; Hamann et al. 2007; Meerburg et al. 2012, 2014; Planck Collaboration XXII 2014; Benetti 2013; Miranda & Hu 2014; Easther & Flauger 2014; Chen & Namjoo 2014; Achúcarro et al. 2014; Hazra et al. 2014a,b, 2016; Hu & Torrado 2015; Planck Collaboration XX 2016; Gruppuso & Sagnotti 2015; Gruppuso et al. 2016; Torrado et al. 2017; Planck Collaboration X 2020; Zeng et al. 2019; Cañas-Herrera et al. 2021; Braglia et al. 2021, 2022b, 2022a; Naik et al. 2022; Hamann & Wons 2022) and bispectra (Planck Collaboration XXIV 2014; Fergusson et al. 2015; Planck Collaboration XVII 2016; Meerburg et al. 2016; Planck Collaboration X 2020; Planck Collaboration IX 2020) have so far not detected any statistically significant signal, setting constraints on deviations from a pure power-law primordial power spectrum (PPS) at the few percent level (depending on the methodology applied) for wavenumbers in the range 0.005 < k/(h Mpc⁻¹)< 0.2. However, it is important to stress that there are interesting models of primordial features that are able to reproduce some of the anomalous features observed in the CMB temperature angular power spectrum at a marginal statistical significance of 99.7%. These include the dip in power in the multipole range ℓ ∼ 20–40 and an oscillatory pattern at the intermediate scales ℓ ∼ 700–800. Models with primordial features have also been studied in the context of the anomalous CMB lensing smoothing excess, as well as the H₀ and S₈ tensions (Planck Collaboration X 2020; Domènech & Kamionkowski 2019; Liu & Huang 2020; Domènech et al. 2020; Keeley et al. 2020; Hazra et al. 2022; Antony et al. 2023; Ballardini & Finelli 2022).

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,b), Ballardini et al. (2016, 2018, 2020), Xu et al. (2016), Ansari Fard & Baghram (2018), Palma et al. (2018), Ballardini (2019), Beutler et al. (2019), Debono et al. (2020), 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.

2. 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 power-law PPS, 𝒫_ℛ, 0(k), with the following parametrisation

$$\begin{array}{c}\hfill {\mathcal{P}}_{\mathcal{R},0}(k)=A_{\mathrm{s}}{\left(\frac{k}{k_{\ast }}\right)}^{n_{\mathrm{s}}-1},\end{array}$$(1)

where A_s and n_s are the amplitude and the spectral index of the comoving curvature perturbations ℛ on superhorizon scales at the pivot scale k_* = 0.05 Mpc⁻¹.

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¹.

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,

$$\begin{array}{c}\hfill S={\displaystyle \int }{\mathrm{d}}^{4}x{a}^3{M}_{\mathrm{Pl}}^2{ϵ}_1[\frac{{\dot{\mathcal{R}}}^2}{c_{\mathrm{s}}^2}-\frac{{({\partial }_i\mathcal{R})}^2}{a^2}]+\cdots ,\end{array}$$(2)

where a(t) is the scale factor, c_s(t) is the speed of sound for curvature perturbations, and ${ϵ}_1(t)\equiv -\dot{H}/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.

Equation (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 (ϵ₁ 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 time-independence 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 ϵ₁(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 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, ϵ₁ 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 ϵ₁ and c_s may be completely different in other models.

3. 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 important 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.

3.1. Spectroscopic galaxy clustering

Following Euclid Collaboration (2020), we define the observed galaxy power spectrum as

$$\begin{array}{cc}\hfill P_{\mathrm{obs}}(k_{\mathrm{ref}},{\mu }_{\mathrm{ref}};z)=& \frac{1}{q_{\perp }^2{q}_{\Vert }}\frac{1}{1+k^2{\mu }^2{f}^2(z){\sigma }_{\mathrm{p}}^2(z)}\hfill \\ \hfill & \times {[b(z){\sigma }_8(z)+f(z){\sigma }_8(z){\mu }^2]}^2\hfill \\ \hfill & \times \frac{P^{\mathrm{IR}\mathrm{res},\mathrm{LO}}(k,\mu ;z)}{{\sigma }_8^2(z)}{F}_z(k,\mu ;z)\hfill \\ \hfill & +P_{\mathrm{s}}(z),\hfill \end{array}$$(3)

where μ is the cosine of the angle of the wave mode with respect to the line of sight pointing into the direction $\widehat{r}$. Here b(z) is the linear clustering bias, f(z) is the growth rate, and σ₈(z) is the root mean square linearly evolved density fluctuations in spheres of 8 h⁻¹ Mpc².

The 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

$$\begin{array}{c}\hfill q_{\perp }(z)=\frac{D_{\mathrm{A}}(z)}{D_{\mathrm{A},\mathrm{ref}}(z)}\text{and}{q}_{\Vert }(z)=\frac{H_{\mathrm{ref}}(z)}{H(z)},\end{array}$$(4)

where the subscript ref refers to the fiducial cosmology. This leads also to a rescaling of the wavevector components as

$$\begin{array}{c}\hfill k_{\perp }=\frac{k_{\perp ,\mathrm{ref}}}{q_{\perp }}\text{and}{k}_{\Vert }=\frac{k_{\Vert ,\mathrm{ref}}}{q_{\Vert }}.\end{array}$$(5)

Using Eqs. (4) and (5) we can convert the known reference cosmology (k_ref and μ_ref) to the true unknown cosmology (k and μ),

$$\begin{array}{cc}& k(k_{\mathrm{ref}},{\mu }_{\mathrm{ref}})=k_{\mathrm{ref}}\alpha ({\mu }_{\mathrm{ref}}),\hfill \end{array}$$(6)

$$\begin{array}{cc}\hfill & \mu ({\mu }_{\mathrm{ref}})=\frac{{\mu }_{\mathrm{ref}}}{\alpha ({\mu }_{\mathrm{ref}})q_{\Vert }},\hfill \end{array}$$(7)

where

$$\begin{array}{c}\hfill \alpha (\mu )={[\frac{{\mu }^2}{q_{\Vert }^2}+\frac{(1-{\mu }^2)}{q_{\perp }^2}]}^{1/2}.\end{array}$$(8)

The relations above can be used to define the effect of the choice of reference cosmology on the observed power spectrum via

$$\begin{array}{c}\hfill P_{\mathrm{obs}}(k_{\mathrm{ref}},{\mu }_{\mathrm{ref}};z)=\frac{1}{q_{\perp }^2{q}_{\Vert }}{P}_{\mathrm{g}}(k(k_{\mathrm{ref}}),\mu ({\mu }_{\mathrm{ref}});z).\end{array}$$(9)

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 redshift-dependent 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

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{p}}^2(z)=\frac{1}{6{\pi }^2}{\displaystyle {\int }_0^{k_{\max }}}\mathrm{d}k{P}_{\mathrm{lin}}(k,z).\end{array}$$(10)

The total galaxy power spectrum in Eq. (3) includes errors on the measurement of the redshift through the exponential factor

$$\begin{array}{c}\hfill F_z(k,\mu ;z)={\text{e}}^{-k^2{\mu }^2{\sigma }_r^2(z)},\end{array}$$(11)

where ${\sigma }_r^2(z)=c(1+z){\sigma }_{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

$$\begin{array}{cc}\hfill F_{\alpha \beta }(z_i)& =\frac{1}{8{\pi }^2}{\displaystyle {\int }_{-1}^1}\mathrm{d}\mu {\displaystyle {\int }_{k_{\min }}^{k_{\max }}}{k}^2\mathrm{d}k{V}_{\mathrm{eff}}(z_i,k)\hfill \\ \hfill & \times \frac{\partial ln{P}_{\mathrm{obs}}(k,\mu ;z_i)}{\partial {\theta }_{\alpha }}\frac{\partial ln{P}_{\mathrm{obs}}(k,\mu ;z_i)}{\partial {\theta }_{\beta }},\hfill \end{array}$$(12)

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

$$\begin{array}{c}\hfill V_{\mathrm{eff}}(k,\mu ;z)=V_{\mathrm{s}}{\left[\frac{n(z)P_{\mathrm{obs}}(k,\mu ;z)}{1+n(z)P_{\mathrm{obs}}(k,\mu ;z)}\right]}^2,\end{array}$$(13)

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

$$\begin{array}{c}\hfill \mathrm{Cov}(k,k^{\prime })\approx \frac{2{(2\pi )}^3}{V_{\mathrm{eff}}(k,\mu ;z)}{P}_{\mathrm{obs}}^2(k,\mu ;z){\delta }^{\mathrm{D}}(k-k^{\prime }),\end{array}$$(14)

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),

$$\begin{array}{c}\hfill {\widehat{P}}_{\mathrm{obs}}(k,\mu ;z)=\frac{1}{\mathrm{\Delta }k}{\displaystyle {\int }_{k-\mathrm{\Delta }k/2}^{k+\mathrm{\Delta }k/2}}\mathrm{d}{k}^{\prime }{P}_{\mathrm{obs}}(k^{\prime },\mu ;z).\end{array}$$(15)

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.0031 h Mpc⁻¹ for a cubic volume for the four redshift bins. We assume a bin width of Δk = 0.004 h Mpc⁻¹ with k_min = 0.002 h Mpc⁻¹. Finally, we can rewrite the Fisher matrix of Eq. (12) for the discrete number of averaged bandpower for the redshift bin z_i as

$$\begin{array}{cc}\hfill F_{\alpha \beta }(z_i)=& \frac{1}{8{\pi }^2}{\displaystyle {\int }_{-1}^1}\mathrm{d}\mu \sum _{k=k_{\min }}^{k_{\max }}{k}^2\mathrm{\Delta }k{V}_{\mathrm{eff}}(z_i,k)\hfill \\ \hfill & \times \frac{\partial ln{\widehat{P}}_{\mathrm{obs}}(k,\mu ;z_i)}{\partial {\theta }_{\alpha }}\frac{\partial ln{\widehat{P}}_{\mathrm{obs}}(k,\mu ;z_i)}{\partial {\theta }_{\beta }}.\hfill \end{array}$$(16)

The total spectroscopic Fisher matrix is then calculated by summing over the redshift bins as $F_{\alpha \beta }^{\mathrm{sp}}={\displaystyle \sum _i}{F}_{\alpha \beta }(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 redshift-dependent parameters bσ₈(z_i) and P_s(z_i) for each of the four redshift bins.

3.2. 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 (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_{ij,\ell }^{\mathit{XY}}$ between probe X in the i-th redshift bin and probe Y in the j-th redshift bin, where the probes X and Y are L for weak lensing or G for photometric galaxy clustering; $C_{ij,\ell }^{\mathit{XY}}$ 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

$$\begin{array}{c}\hfill C_{ij,\ell }^{\mathit{XY}}=c{\displaystyle {\int }_{z_{\min }}^{z_{\max }}}\mathrm{d}z\frac{W_i^X(z)W_j^Y(z)}{H(z)r^2(z)}{P}_{\mathrm{NL}}(k_{\ell },z),\end{array}$$(17)

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

$$\begin{array}{cc}\hfill W_i^{\mathrm{G}}(z)=& b_i(z)n_i(z)\frac{H(z)}{c},\hfill \end{array}$$(18)

$$\begin{array}{cc}\hfill W_i^{\mathrm{L}}(z)=& \frac{3}{2}{\mathrm{\Omega }}_{\mathrm{m},0}{\left(\frac{H_0}{c}\right)}^2(1+z)r(z){\displaystyle {\int }_z^{z_{\max }}}\mathrm{d}{z}^{\prime }{n}_i(z^{\prime })[1-\frac{r(z)}{r(z^{\prime })}]\hfill \\ \hfill & +W_i^{\mathrm{IA}}(z),\hfill \end{array}$$(19)

where b_i(z) is the galaxy bias in the i-th redshift bin, and $W_i^{\text{IA}}(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

$$\begin{array}{c}\hfill n_i(z)=\frac{{\int }_{z_i^{-}}^{z_i^{+}}\mathrm{d}{z}^{\prime }n(z)p_{\mathrm{ph}}(z^{\prime }|z)}{{\int }_{z_{\min }}^{z_{\max }}\mathrm{d}z{\displaystyle {\int }_{z_i^{-}}^{z_i^{+}}}\mathrm{d}{z}^{\prime }n(z)p_{\mathrm{ph}}(z^{\prime }|z)},\end{array}$$(20)

where $(z_i^{-},z_i^{+})$ are the edges of the ith 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 (2020).

The IA contribution is computed following the extended nonlinear alignment (eNLA) model adopted in Euclid Collaboration (2020) so that the corresponding window function is

$$\begin{array}{c}\hfill W_i^{\mathrm{IA}}(z)=-\frac{{\mathcal{A}}_{\mathrm{IA}}{\mathcal{C}}_{\mathrm{IA}}{\mathrm{\Omega }}_{\mathrm{m},0}{\mathcal{F}}_{\mathrm{IA}}(z)}{D(z)}{n}_i(z)\frac{H(z)}{c},\end{array}$$(21)

where

$$\begin{array}{c}\hfill {\mathcal{F}}_{\mathrm{IA}}(z)={(1+z)}^{{\eta }_{\mathrm{IA}}}{\left[\frac{⟨L⟩(z)}{L_{\star }(z)}\right]}^{{\beta }_{\mathrm{IA}}},\end{array}$$(22)

with ⟨L⟩(z) the redshift-dependent mean, and L_⋆(z) the characteristic luminosity of source galaxies as computed from the luminosity function. 𝒜_IA, β_IA and η_IA are the nuisance parameters of the model, and 𝒞_IA is a constant accounting for dimensional units.

We consider a Gaussian-only covariance whose elements are given by

$$\begin{array}{cc}\hfill \text{Cov}[C_{ij,\ell }^{\mathrm{AB}},C_{mn,{\ell }^{\prime }}^{\mathrm{CD}}]=& \frac{{\delta }_{\ell {\ell }^{\prime }}^{\mathrm{K}}}{(2\ell +1)f_{\mathrm{sky}}\mathrm{\Delta }\ell }\hfill \\ \hfill & \times \{[C_{im,\ell }^{\mathrm{AC}}+{\mathcal{N}}_{im,\ell }^{\mathrm{AC}}][C_{jn,{\ell }^{\prime }}^{\mathrm{BD}}+{\mathcal{N}}_{jn,{\ell }^{\prime }}^{\mathrm{BD}}]\hfill \\ \hfill & +[C_{in,\ell }^{\mathrm{AD}}+{\mathcal{N}}_{in,\ell }^{\mathrm{AD}}][C_{jm,{\ell }^{\prime }}^{\mathrm{BC}}+{\mathcal{N}}_{jm,{\ell }^{\prime }}^{\mathrm{BC}}]\},\hfill \end{array}$$(23)

where the upper- and lower-case Latin indices run over L and G (all tomographic bins), ${\delta }_{\mathcal{l}{\mathcal{l}}^{\prime }}^{\text{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,

$$\begin{array}{cc}\hfill {\mathcal{N}}_{ij,\ell }^{\mathrm{LL}}& =\frac{{\delta }_{\mathit{ij}}^{\mathrm{K}}}{{\overline{n}}_i}{\sigma }_{ϵ}^2,{\mathcal{N}}_{ij,\ell }^{\mathrm{GG}}=\frac{{\delta }_{\mathit{ij}}^{\mathrm{K}}}{{\overline{n}}_i},{\mathcal{N}}_{ij,\ell }^{\mathrm{GL}}=0,\hfill \end{array}$$(24)

where σ_ϵ is the variance of observed ellipticities.

For evaluating the Fisher matrix $F_{\alpha \beta }^{\mathrm{ph}}$ for the observed galaxy power spectrum, we use

$$\begin{array}{c}\hfill F_{\alpha \beta }^{\mathrm{ph}}={\displaystyle \sum _{\ell ={\ell }_{\min }}^{{\ell }_{\max }}\sum _{ij,mn}}\frac{\partial C_{ij,\ell }^{\mathrm{AB}}}{\partial {\theta }_{\alpha }}{\text{Cov}}^{-1}[C_{ij,\ell }^{\mathrm{AB}},C_{mn,\ell }^{\mathrm{CD}}]\frac{\partial C_{mn,\ell }^{\mathrm{CD}}}{\partial {\theta }_{\beta }}.\end{array}$$(25)

3.3. Survey specifications

Following Euclid Collaboration (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⁻⁴.

For the photometric probes, the sources are split into 10 equi-populated redshift bins whose limits are obtained from the redshift distribution

$$\begin{array}{c}\hfill n(z)\propto {\left(\frac{z}{z_0}\right)}^2\text{exp}[-{\left(\frac{z}{z_0}\right)}^{3/2}],\end{array}$$(26)

with $z_0=0.9/\sqrt{2}$ and the normalisation set by the requirement that the surface density of galaxies is ${\overline{n}}_g=30{\mathrm{arcmin}}^{-2}$. This is then convolved with the sum of two Gaussians to account for the effect of photometric redshift (see Euclid Collaboration 2020, for details). The galaxy bias is assumed to be constant within each redshift bin, with fiducial values $b_i=\sqrt{1+{\overline{z}}_i}$, where ${\overline{z}}_i$ 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 2022b).

4. 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.

4.1. Models

We consider a superimposed pattern of oscillations as

$$\begin{array}{c}\hfill {\mathcal{P}}_{\mathcal{R}}(k)={\mathcal{P}}_{\mathcal{R},0}(k)[1+\delta P^{\mathrm{X}}(k)],\end{array}$$(27)

where 𝒫_ℛ, 0(k) is the standard power-law PPS of the comoving curvature perturbations ℛ on superhorizon scales, given in Eq. (1). We study the following templates with superimposed oscillations on the PPS

$$\begin{array}{c}\hfill \delta P^{\mathrm{X}}(k)={\mathcal{A}}_{\mathrm{X}}sin({\omega }_{\mathrm{X}}{\mathrm{\Xi }}_{\mathrm{X}}+2\pi {\varphi }_{\mathrm{X}}),\end{array}$$(28)

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 VI 2020): ω_b = 0.02237, ω_c = 0.1200, H₀ = 67.36 km s⁻¹ Mpc⁻¹, ln(10¹⁰A_s)=3.044, n_s = 0.9649, σ₈ = 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:

   $$\begin{array}{c}\hfill {\mathrm{\Theta }}_{\mathrm{lin}}=\{{\mathcal{A}}_{\mathrm{lin}}=0.01,{\omega }_{\mathrm{lin}}=10,{\varphi }_{\mathrm{lin}}=0\}.\end{array}$$(29)

2. Logarithmic oscillations:

   $$\begin{array}{c}\hfill {\mathrm{\Theta }}_{\log }=\{{\mathcal{A}}_{\log }=0.01,{\omega }_{\log }=10,{\varphi }_{\log }=0\}.\end{array}$$(30)

We show the effect of the feature parameters on the PPS in Fig. 1.

Fig. 1. Primordial power spectrum of curvature perturbations for linear (LIN) and logarithmic (LOG) feature models with fiducial parameter values given in Eqs. (29) and (30).

4.2. Modelling nonlinear scales in perturbation theory in the presence of primordial oscillatory features

The linearly propagated matter power spectrum is given by

$$\begin{array}{c}\hfill P_{\mathrm{lin}}(k,z)=M^2(k,z)P_{\mathrm{\Phi }}(k),\end{array}$$(31)

where $M(k,z)=2k^2{c}^{2}T(k)D(z)/(3{\mathrm{\Omega }}_{\text{m}}{H}_0^2)$, 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π²)𝒫_ℛ(k)/k³.

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,

$$\begin{array}{c}\hfill P_{\mathrm{lin}}(k,z)=D^2(z)[P_{\mathrm{nw}}(k)+P_{\mathrm{w}}(k)],\end{array}$$(32)

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

$$\begin{array}{c}\hfill P_{\mathrm{w}}(k)\equiv P_{\mathrm{nw}}[\delta P_{\mathrm{w}}^{\mathrm{BAO}}(k)+\delta P_{\mathrm{w}}^{\mathrm{X}}(k)+\delta P_{\mathrm{w}}^{\mathrm{BAO}}(k)\delta P_{\mathrm{w}}^{\mathrm{X}}(k)].\end{array}$$(33)

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 𝒜_BAO × 𝒜_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)

$$\begin{array}{cc}& P^{\mathrm{IR}\mathrm{res},\mathrm{LO}+\mathrm{NLO}}(k,z)=D^2(z)P_{\mathrm{nw}}(k)\hfill \\ \hfill & \times \{1+[1+k^2{D}^2(z){\mathrm{\Sigma }}_{\mathrm{BAO}}^2]{\mathrm{e}}^{-k^2{D}^2(z){\mathrm{\Sigma }}_{\mathrm{BAO}}^2}\delta P_{\mathrm{w}}^{\mathrm{BAO}}(k)\hfill \\ \hfill & +[1+k^2{D}^2(z){\mathrm{\Sigma }}_{\mathrm{lin}}^2]{\mathrm{e}}^{-k^2{D}^2(z){\mathrm{\Sigma }}_{\mathrm{lin}}^2}\delta P_{\mathrm{w}}^{\mathrm{lin}}(k)\}\hfill \\ \hfill & +D^4(z)P^{1-\mathrm{loop}}\left[P^{\mathrm{IR}\mathrm{res},\mathrm{LO}}(k)\right],\hfill \end{array}$$(34)

and for logarithmic oscillations as (Vasudevan et al. 2019; Beutler et al. 2019)

$$\begin{array}{cc}& P^{\mathrm{IR}\mathrm{res},\mathrm{LO}+\mathrm{NLO}}(k,z)=D^2(z)P_{\mathrm{nw}}(k)\hfill \\ \hfill & \times \{1+[1+k^2{D}^2(z){\mathrm{\Sigma }}_{\mathrm{BAO}}^2]{\mathrm{e}}^{-k^2{D}^2(z){\mathrm{\Sigma }}_{\mathrm{BAO}}^2}\delta P_{\mathrm{w}}^{\mathrm{BAO}}(k)\hfill \\ \hfill & +[1+k^2{D}^2(z){\mathrm{\Sigma }}_{\log }^2]{\mathrm{e}}^{-k^2{D}^2(z){\mathrm{\Sigma }}_{\log }^2}[cos\left(k^2{D}^2(z){\widehat{\mathrm{\Sigma }}}_{\log }^2\right)\delta P_{\mathrm{w}}^{\log }(k)\hfill \\ \hfill & -sin\left(k^2{D}^2(z){\widehat{\mathrm{\Sigma }}}_{\log }^2\right)\frac{\mathrm{d}\delta P_{\mathrm{w}}^{\log }(k)}{{\omega }_{\log }\mathrm{d}lnk}]\}\hfill \\ \hfill & +D^4(z)P^{1-\mathrm{loop}}\left[P^{\mathrm{IR}\mathrm{res},\mathrm{LO}}(k)\right].\hfill \end{array}$$(35)

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₂₂ + 2P₁₃ with

$$\begin{array}{cc}& P_{22}(k)={\displaystyle \int \frac{{\mathrm{d}}^{3}q}{4{\pi }^3}}{F}_2^2(\mathbf{q},\mathbf{k}-\mathbf{q})P_{\mathrm{lin}}(q)P_{\mathrm{lin}}(|\mathbf{k}-\mathbf{q}|),\hfill \end{array}$$(36)

$$\begin{array}{cc}\hfill & P_{13}(k)=3P_{\mathrm{lin}}(k){\displaystyle \int \frac{{\mathrm{d}}^{3}q}{4{\pi }^3}}{F}_3(\mathbf{q},-\mathbf{q},\mathbf{k})P_{\mathrm{lin}}(q),\hfill \end{array}$$(37)

while evaluating the loop integrals P₂₂ and P₁₃ 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)³. The IR resummed power spectrum at leading order is given by

$$\begin{array}{cc}& P^{\mathrm{IR}\mathrm{res},\mathrm{LO}}(k,z)=D^2(z)P_{\mathrm{nw}}(k)\hfill \\ \hfill & \times [1+{\mathrm{e}}^{-k^2{D}^2(z){\mathrm{\Sigma }}_{\mathrm{BAO}}^2}\delta P_{\mathrm{w}}^{\mathrm{BAO}}(k)+{\mathrm{e}}^{-k^2{D}^2(z){\mathrm{\Sigma }}_{\mathrm{lin}}^2}\delta P_{\mathrm{w}}^{\mathrm{lin}}(k)]\hfill \end{array}$$(38)

for the linear oscillations (Blas et al. 2016b; Beutler et al. 2019) and by

$$\begin{array}{cc}\hfill P^{\mathrm{IR}\mathrm{res},\mathrm{LO}}(k,z)& =D^2(z)P_{\mathrm{nw}}(k)\hfill \\ \hfill & \times [1+{\mathrm{e}}^{-k^2{D}^2(z){\mathrm{\Sigma }}_{\mathrm{BAO}}^2}\delta P_{\mathrm{w}}^{\mathrm{BAO}}(k)\hfill \\ \hfill & +{\mathrm{e}}^{-k^2{D}^2(z){\mathrm{\Sigma }}_{\log }^2}cos\left(k^2{D}^2(z){\widehat{\mathrm{\Sigma }}}_{\log }^2\right)\delta P_{\mathrm{w}}^{\log }(k)\hfill \\ \hfill & -{\mathrm{e}}^{-k^2{D}^2(z){\mathrm{\Sigma }}_{\log }^2}sin\left(k^2{D}^2(z){\widehat{\mathrm{\Sigma }}}_{\log }^2\right)\frac{\mathrm{d}\delta P_{\mathrm{w}}^{\log }(k)}{{\omega }_{\log }\mathrm{d}lnk}]\hfill \end{array}$$(39)

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

$$\begin{array}{cc}& {\mathrm{\Sigma }}_{\mathrm{BAO}}^2(k_{\mathrm{S}})\equiv {\displaystyle {\int }_0^{k_{\mathrm{S}}}}\frac{\mathrm{d}q}{6{\pi }^2}{P}_{\mathrm{nw}}(q)[1-j_0\left(q{r}_{\mathrm{s}}\right)+2j_2\left(q{r}_{\mathrm{s}}\right)],\hfill \end{array}$$(40)

$$\begin{array}{cc}\hfill & {\mathrm{\Sigma }}_{\mathrm{lin}}^2({\omega }_{\mathrm{lin}},k_{\mathrm{S}})\equiv \hfill \\ \hfill & {\displaystyle {\int }_0^{k_{\mathrm{S}}}}\frac{\mathrm{d}q}{6{\pi }^2}{P}_{\mathrm{nw}}(q)[1-j_0\left(q\frac{{\omega }_{\mathrm{lin}}}{k_{\ast }}\right)+2j_2\left(q\frac{{\omega }_{\mathrm{lin}}}{k_{\ast }}\right)],\hfill \end{array}$$(41)

$$\begin{array}{cc}\hfill & {\mathrm{\Sigma }}_{\log }^2(k,{\omega }_{\log },k_{\mathrm{S}})\equiv {\displaystyle {\int }_0^{k_S}}\frac{\mathrm{d}q}{4{\pi }^2}{P}_{\mathrm{nw}}(q)\hfill \\ \hfill & \times {\displaystyle {\int }_{-1}^1}\mathrm{d}\mu {\mu }^2\{1-cos[{\omega }_{\log }ln(1-\frac{q\mu }{k})]\},\hfill \end{array}$$(42)

$$\begin{array}{cc}\hfill & {\widehat{\mathrm{\Sigma }}}_{\log }^2(k,{\omega }_{\log },k_{\mathrm{S}})\equiv \hfill \\ \hfill & -{\displaystyle {\int }_0^{k_S}}\frac{\mathrm{d}q}{4{\pi }^2}{P}_{\mathrm{nw}}(q){\displaystyle {\int }_{-1}^1}\mathrm{d}\mu {\mu }^{2}sin[{\omega }_{\log }ln(1-\frac{q\mu }{k})],\hfill \end{array}$$(43)

where j_n are spherical Bessel functions, r_s ≃ 147 Mpc is the scale setting the period of the BAO (Planck Collaboration VI 2020), 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) D(z)Σ_X → D(z)Σ_X(μ) with

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\mathrm{X}}^2(k,\mu ,z)={\mathrm{\Sigma }}_{\mathrm{X}}^2(z)\{1-{\mu }^2+{\mu }^2{[1+f(k,z)]}^2\}.\end{array}$$(44)

We note that while in the next subsection, Sect. 4.3, we present both PT results at LO and NLO, that is Eqs. (38) and (39)–(34) and (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

$$\begin{array}{c}\hfill P_{\mathrm{NL}}(k,z)=\frac{P_{\mathrm{NL}}^{\mathrm{\Lambda }\mathrm{CDM}}(k,z)}{P_{\mathrm{lin}}^{\mathrm{\Lambda }\mathrm{CDM}}(k,z)}{P}^{\mathrm{IR}\mathrm{res},\mathrm{LO}}(k,z).\end{array}$$(45)

4.3. 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, 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)⁴. We 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³ dark matter particles in a comoving box with side length of 1024 h⁻¹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)⁵, for 𝒜_lin = 0.03 and ω_lin = 10.

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 𝒜lin = 0.03 and ωlin = 10 (Ballardini et al. 2020).

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.

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 kS = ϵ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).

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₁₀ω_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₁₀ω_lin > 0.6 and between 1% and 2% for log₁₀ω_lin ≤ 0.6. At NLO, we find an agreement better than 1% for all the frequencies and differences < 0.5% for log₁₀ω_lin > 0.4. For the template with undamped logarithmic oscillations, we find a similar agreement for the frequencies log₁₀ω_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 N-body 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.

4.4. 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 primordial 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⁻¹)≲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⁻¹ or s ≳ 15–20 h⁻¹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

$$\begin{array}{c}\hfill \rho (\mathit{x}){\mathrm{d}}^3\mathit{x}=\rho (\mathit{q}){\mathrm{d}}^3\mathit{q}\approx \overline{\rho }{\mathrm{d}}^3\mathit{q},\end{array}$$(46)

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. $\rho (\mathit{q})\simeq \overline{\rho }$.

We define the displacement field Ψ between the final and initial positions of a particle as Ψ(x)=x − q, which can be rewritten as

$$\begin{array}{c}\hfill {\mathbf{\nabla }}_{\mathit{x}}\mathrm{\Theta }(\mathit{x})\equiv \mathit{q}=\mathit{x}-\mathbf{\Psi }(\mathit{x}),\end{array}$$(47)

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⁶. Combining Eqs. (47) and (46) gives

$$\begin{array}{c}\hfill det[{\mathrm{\nabla }}^i{\mathrm{\nabla }}_j\mathrm{\Theta }(\mathbf{x})]=\frac{\rho (\mathbf{x})}{\overline{\rho }}\equiv 1+\delta (\mathbf{x}),\end{array}$$(48)

where $\delta (\mathit{x})=\rho (\mathit{x})/\overline{\rho }-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

$$\begin{array}{c}\hfill {\delta }_{\mathrm{r}}=-{\mathrm{\nabla }}_{\mathit{q}}·\mathbf{\Psi }(\mathit{q}),\end{array}$$(49)

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 Tessellation Field Estimator code (DTFE; Schaap & van de Weygaert 2000; van de Weygaert & Schaap 2009; Cautun & van de Weygaert 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_{\mathrm{w}}^{\mathrm{undamped}}$ by a Gaussian function

$$\begin{array}{c}\hfill O_{\mathrm{w}}^{\mathrm{damped}}(k,z)=O_{\mathrm{w}}^{\mathrm{undamped}}(k,z)exp[-\frac{k^2\zeta {(z)}^2}{2}],\end{array}$$(50)

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 also performed reconstruction from halo 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. Figure B.1 shows that reconstruction can help recover the weakened wiggles down to k ≃ 1 h Mpc⁻¹. 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⁻¹.

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_{\mathrm{w}}^{\mathrm{damped}}(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_{\mathrm{w}}^{\mathrm{damped}}$ and $O_{\mathrm{w}}^{\mathrm{undamped}}$ 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 low-frequency featured models due to the very few wiggles in the k range of fitting, k = (0.05 − 1) h Mpc⁻¹, 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.

5. 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,

$$\begin{array}{c}\hfill {\mathrm{\Theta }}_{\mathrm{final}}=\{{\mathrm{\Omega }}_{\mathrm{m},0},{\mathrm{\Omega }}_{\mathrm{b},0},h,n_{\mathrm{s}},{\sigma }_8,{\mathcal{A}}_X,{\omega }_X,{\varphi }_X\},\end{array}$$(51)

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⁻¹, 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⁻¹, 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 𝒜_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 𝒜_lin = 0.0100 ± 0.0023 (±0.0023) and 𝒜_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 𝒜_lin = 0.0100 ± 0.0021 (±0.0018) and 𝒜_log = 0.0100 ± 0.0022 (±0.0018) at a 68.3% CL for the pessimistic (optimistic) setting.

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 GCsp (green), WL+GCph+XC (orange), and their combination (blue).

For GC_sp, bounds on the primordial amplitude parameter 𝒜_X strongly depend on the primordial frequency value ω_X; see Huang et al. (2012), Ballardini et al. (2016, 2020), Slosar et al. (2019), Beutler et al. (2019). 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 high-frequency 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₁₀ω_lin ∼ 0.87, are degraded. For the combination of photometric probes, bounds on the primordial amplitude parameter 𝒜_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₁₀ω_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.003 on the primordial feature amplitude for the entire range of frequencies.

Fig. 5. Marginalised uncertainties on 𝒜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 GCsp (blue), WL+GCph+XC (orange), and their combination (green). Dashed (solid) curves correspond to the pessimistic (optimistic) setting.

The nonlinear reconstruction described in Sect. 4.4 applied to GC_sp reduces significantly the uncertainties on the feature amplitude, going from σ(𝒜_lin) = 0.0023 (0.0023) at a 68.3% CL for the pessimistic (optimistic) setting to σ(𝒜_lin) = 0.0014 (0.0013), and going from σ(𝒜_log) = 0.0025 (0.0025) to σ(𝒜_log) = 0.0016 (0.0012). Combining with the photometric probes, we find σ(𝒜_lin) = 0.0014 (0.0012) at a 68.3% CL and σ(𝒜_log) = 0.0015 (0.0012) for the pessimistic (optimistic) setting, compared with σ(𝒜_lin) = 0.0021 (0.0018) at a 68.3% CL and σ(𝒜_log) = 0.0022 (0.0018) for the same settings without reconstruction. We also note that in the analysis here k_max is set to be 0.25–0.30 h Mpc⁻¹, while in Li et al. (2022) k_max = 0.50 h Mpc⁻¹ is used.

6. 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_{\mathrm{NL}}^{\mathrm{lin}}$, the phase ${\varphi }_{\text{lin}}^B$, and the frequency ω_lin, and it is given by (Chen et al. 2007, 2008)

$$\begin{array}{c}\hfill B_{\mathrm{\Phi }}^{\mathrm{lin}}(k_1,k_2,k_3)=f_{\mathrm{NL}}^{\mathrm{lin}}\frac{6A^2}{k_1^2{k}_2^2{k}_3^2}sin[{\omega }_{\mathrm{lin}}\frac{K}{k_{\ast }}+2\pi {\varphi }_{\mathrm{lin}}^B],\end{array}$$(52)

where K = k₁ + k₂ + k₃ and $A=9/25(2{\pi }^2)k_{*}^{1-n_{\text{s}}}{A}_{\text{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 logarithmically-spaced 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)

$$\begin{array}{c}\hfill B_{\mathrm{\Phi }}^{\log }(k_1,k_2,k_3)=f_{\mathrm{NL}}^{\log }\frac{6A^2}{k_1^2{k}_2^2{k}_3^2}sin[{\omega }_{\log }ln\left(\frac{K}{k_{\ast }}\right)+2\pi {\varphi }_{\log }^B],\end{array}$$(53)

where ${\varphi }_{\text{log}}^B$ and $f_{\mathrm{NL}}^{\log }$ 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

$$\begin{array}{cc}\hfill P_g(\mathit{k},z)& =D^P(\mathit{k},z)Z_1^2(\mathit{k},z)P^{\mathrm{lin}}(k,z)+P_{\epsilon }(z),\hfill \end{array}$$(54)

$$\begin{array}{c}\hfill \\ \hfill & B_g({\mathit{k}}_1,{\mathit{k}}_2,{\mathit{k}}_3,z)=D^B({\mathit{k}}_1,{\mathit{k}}_2,{\mathit{k}}_3,z)\hfill \\ \hfill & \times [Z_1({\mathit{k}}_1,z)Z_1({\mathit{k}}_2,z)Z_1({\mathit{k}}_3,z)B_I({\mathit{k}}_1,{\mathit{k}}_2,{\mathit{k}}_3,z)\hfill \\ \hfill & +2Z_1({\mathit{k}}_1,z)Z_1({\mathit{k}}_2,z)Z_2({\mathit{k}}_1,{\mathit{k}}_2,z)P_{\mathrm{lin}}(k_1,z)P_{\mathrm{lin}}(k_2,z)+2\text{perm}]\hfill \\ \hfill & +2P_{\epsilon {\epsilon }_{\delta }}[Z_1({\mathit{k}}_1,z)P_{\mathrm{lin}}(k_1,z)+2\text{perm}]+B_{\epsilon },\hfill \end{array}$$(55)

where the linearly propagated primordial oscillatory bispectrum is

$$\begin{array}{c}\hfill B_I({\mathit{k}}_1,{\mathit{k}}_2,{\mathit{k}}_3,z)=M(k_1,z)M(k_2,z)M(k_3,z)B_{\mathrm{\Phi }}^X(k_1,k_2,k_3,z)\end{array}$$(56)

with $B_{\mathrm{\Phi }}^X$ given by Eqs. (52) and (53) for the linear and logarithmic oscillations, respectively. The redshift space kernels are given by

$$\begin{array}{cc}\hfill Z_1({\mathit{k}}_i,z)& =b_1(z)+f{\mu }_i^2,\hfill \end{array}$$(57)

$$\begin{array}{cc}\hfill Z_2({\mathit{k}}_i,{\mathit{k}}_j,z)& =b_1(z)F_2({\mathit{k}}_i,{\mathit{k}}_j)+f(z){\mu }_{\mathit{ij}}^2{G}_2({\mathit{k}}_i,{\mathit{k}}_j)+\frac{b_2(z)}{2}\hfill \\ \hfill & +\frac{b_{s^2}(z)}{2}{S}_2({\mathit{k}}_i,{\mathit{k}}_j)+\frac{f(z){\mu }_{\mathit{ij}}{k}_{\mathit{ij}}}{2}[\frac{{\mu }_i}{k_i}{Z}_1({\mathit{k}}_j,z)+\frac{{\mu }_j}{k_j}{Z}_1({\mathit{k}}_i,z)],\hfill \end{array}$$(58)

where ${\mu }_i={\mathit{k}}_i·\widehat{z}/k_i$ is the cosine of the angle between the wave vector k_i and the line of sight, μ = μ₁ + …μ_n, and k = |k₁ + …k_n|. Specifically μ₁₂ = (μ₁k₁ + μ₂k₂)/k₁₂ and $k_{12}^2={({\mathit{k}}_1+{\mathit{k}}_2)}^2$. The kernels F₂(k_i, k_j) and G₂(k_i, k_j) are the second-order symmetric PT kernels (Bernardeau et al. 2002), while $S_2({\mathit{k}}_1,{\mathit{k}}_2)={({\mathit{k}}_1\cdot {\mathit{k}}_2)}^2/(k_1^2{k}_2^2)-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²(z)= − 4[b₁(z)−1]/7. The redshift-space bispectrum is characterised by five variables: three to define the triangle shape (the sides k₁, k₂, k₃) and two that characterise the orientation of the triangle relative to the line of sight, i.e. B_g(k₁, k₂, k₃, z)=B_g(k₁, k₂, k₃, μ₁, ϕ, z). 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_{\mathrm{NL}}^{\mathrm{X}}$), shown in Barreira (2020), does not affect our forecasts.

The multiplicative factors D^P(k, z) and D^B(k₁, k₂, k₃, z), which incorporate errors on redshift measurements and FOG effect, are $D^P(\mathit{k},z)=\exp [-k\mu ({\sigma }_{\text{p}}^2(z)+{\sigma }_{\text{r}}^2(z))]$ and $D^B({\mathit{k}}_1,{\mathit{k}}_2,{\mathit{k}}_3,z)=\exp [-(k_1^2{\mu }_1^2+k_2^2{\mu }_2^2+k_3^2{\mu }_3^2)({\sigma }_{\text{p}}^2(z)+{\sigma }_{\text{r}}^2(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_{\epsilon }=1/{\overline{n}}_{\mathrm{g}}$, $P_{\epsilon {\epsilon }_{\delta }}=b_1/(2{\overline{n}}_{\mathrm{g}})$, and $B_{\epsilon }=1/{\overline{n}}_{\mathrm{g}}^2$.

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)

$$\begin{array}{cc}\hfill B_{\mathrm{g}}^{\mathrm{obs}}(k_{\mathrm{ref},1}& ,k_{\mathrm{ref},2},k_{\mathrm{ref},3},{\mu }_{\mathrm{ref},1},\varphi ;z)\hfill \\ \hfill & =\frac{1}{q_{\Vert }^2{q}_{\perp }^4}{B}_{\mathrm{g}}(k(k_{\mathrm{ref},1}),k(k_{\mathrm{ref},2}),k(k_{\mathrm{ref},3}),\mu ({\mu }_{\mathrm{ref},1}),\varphi ;z).\hfill \end{array}$$(59)

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

$$\begin{array}{cc}\hfill F_{\alpha \beta }^B(z_i)& =\frac{1}{4\pi }{\displaystyle \sum _{k_1,k_2,k_3}}{\displaystyle {\int }_{-1}^1}\mathrm{d}{\mu }_1{\displaystyle {\int }_0^{2\pi }}\mathrm{d}\varphi \frac{1}{\mathrm{\Delta }{B}^2({\mathit{k}}_1,{\mathit{k}}_2,{\mathit{k}}_3,z_i)}\hfill \\ \hfill & \times \frac{\partial B_{\mathrm{g}}^{\mathrm{obs}}({\mathit{k}}_1,{\mathit{k}}_2,{\mathit{k}}_3,z_i)}{\partial p_{\alpha }}\frac{\partial B_{\mathrm{g}}^{\mathrm{obs}}({\mathit{k}}_1,{\mathit{k}}_2,{\mathit{k}}_3,z_i)}{\partial p_{\beta }},\hfill \end{array}$$(60)