# Euclid: The search for primordial features

M Ballardini, Y Akrami, F Finelli, D Karagiannis, B Li, Y Li, Z Sakr, D Sapone, A Achúcarro, M Baldi, et al.

► To cite this version:

M Ballardini, Y Akrami, F Finelli, D Karagiannis, B Li, et al. Euclid: The search for primordial features. Astron.Astrophys., 2024, 683, pp.A220. 10.1051/0004-6361/202348162. hal-04237840

HAL Id: hal-04237840

https://hal.science/hal-04237840

Submitted on 19 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Euclid: The search for primordial features

M Ballardini, Y Akrami, F Finelli, D Karagiannis, B Li, Y Li, Z Sakr, D Sapone, A Achúcarro, M Baldi, et al.

► To cite this version:

M Ballardini, Y Akrami, F Finelli, D Karagiannis, B Li, et al. Euclid: The search for primordial features. Astronomy & Astrophysics - A&A, 2024, 683, pp.A220. 10.1051/0004-6361/202348162. hal-04237840

HAL Id: hal-04237840

https://hal.science/hal-04237840v1

Submitted on 19 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Euclid: The search for primordial features*

M. Ballardini<sup>1,2,3</sup> ◦, Y. Akrami<sup>4,5</sup>, F. Finelli<sup>3,6</sup>, D. Karagiannis<sup>7</sup>, B. Li<sup>8</sup>, Y. Li<sup>9</sup>, Z. Sakr<sup>10,11,12</sup>, D. Sapone<sup>13</sup>, A. Achúcarro<sup>14,15</sup>, M. Baldi<sup>16,3,17</sup>, N. Bartolo<sup>18,19,20</sup>, G. Cañas-Herrera<sup>21,14</sup>, S. Casas<sup>22</sup>, R. Murgia<sup>23,24</sup>, H.-A. Winther<sup>25</sup>, M. Viel<sup>26,27,28,29</sup>, A. Andrews<sup>3</sup>, J. Jasche<sup>30,31</sup>, G. Lavaux<sup>32</sup> ◦, D. K. Hazra<sup>3,33,34</sup>, D. Paoletti<sup>3,6</sup>, J. Valiviita<sup>35,36</sup>, A. Amara<sup>37</sup>, S. Andreon<sup>38</sup>, N. Auricchio<sup>3</sup>, P. Battaglia<sup>3</sup>, D. Bonino<sup>39</sup>, E. Branchini<sup>40,41</sup>, M. Brescia<sup>42,43,44</sup>, J. Brinchmann<sup>45</sup>, S. Camera<sup>46,47,39</sup>, V. Capobianco<sup>39</sup>, C. Carbone<sup>48</sup>, J. Carretero<sup>49,50</sup>, M. Castellano<sup>51</sup>, S. Cavuoti<sup>43,44</sup>, A. Cimatti<sup>52</sup>, G. Congedo<sup>53</sup>, L. Conversi<sup>54,55</sup>, Y. Copin<sup>56</sup>, L. Corcione<sup>39</sup>, F. Courbin<sup>57</sup>, H. M. Courtois<sup>58</sup>, A. Da Silva<sup>59,60</sup>, H. Degaudenzi<sup>61</sup>, F. Dubath<sup>61</sup>, X. Dupac<sup>55</sup>, M. Farina<sup>62</sup>, S. Farrens<sup>63</sup>, M. Frailis<sup>27</sup>, E. Franceschi<sup>3</sup>, M. Fumana<sup>48</sup>, S. Galeotta<sup>27</sup>, B. Gillis<sup>53</sup>, C. Giocoli<sup>3,17</sup>, A. Grazian<sup>20</sup>, F. Grupp<sup>64,65</sup>, S. V. H. Haugan<sup>25</sup>, W. Holmes<sup>66</sup>, F. Hormuth<sup>67</sup>, A. Hornstrup<sup>68,69</sup>, P. Hudelot<sup>32</sup>, K. Jahnke<sup>70</sup>, S. Kermiche<sup>71</sup>, A. Kiessling<sup>66</sup>, M. Kunz<sup>72</sup>, H. Kurki-Suonio<sup>35,36</sup>, P. B. Lilje<sup>25</sup>, V. Lindholm<sup>35,36</sup>, I. Lloro<sup>73</sup>, E. Maiorano<sup>3</sup>, O. Mansutti<sup>27</sup>, O. Marggraf<sup>74</sup>, N. Martinet<sup>75</sup>, F. Marulli<sup>76,3,17</sup>, R. Massey<sup>8</sup>, E. Medinaceli<sup>3</sup>, S. Mei<sup>77</sup>, Y. Mellier<sup>31,32</sup>, M. Meneghetti<sup>3,17</sup>, E. Merlin<sup>51</sup>, G. Meylan<sup>57</sup>, M. Moresco<sup>76,3</sup>, L. Moscardini<sup>76,3,17</sup>, E. Munari<sup>27</sup>, S.-M. Niemi<sup>21</sup>, C. Padilla<sup>49</sup>, S. Paltani<sup>61</sup>, F. Pasian<sup>27</sup>, K. Pedersen<sup>78</sup>, W. J. Percival<sup>79,80,81</sup>, V. Pettorino<sup>82</sup>, S. Pires<sup>63</sup>, G. Polenta<sup>83</sup>, M. Poncet<sup>84</sup>, L. A. Popa<sup>85</sup>, L. Pozzetti<sup>3</sup>, F. Raison<sup>64</sup>, A. Renzi<sup>18,19</sup>, J. Rhodes<sup>66</sup>, G. Riccio<sup>43</sup>, E. Romelli<sup>27</sup>, M. Roncarelli<sup>3</sup>, R. Saglia<sup>86,64</sup>, B. Sartoris<sup>86,27</sup>, T. Schrabback<sup>87</sup>, A. Secroun<sup>71</sup>, G. Seidel<sup>70</sup>, S. Serrano<sup>88,89,90</sup>, C. Sirignano<sup>18,19</sup>, G. Sirri<sup>17</sup>, L. Stanco<sup>19</sup>, J.-L. Starck<sup>91</sup>, C. Surace<sup>75</sup>, P. Tallada-Crespi<sup>92,50</sup>, A. N. Taylor<sup>53</sup>, I. Tereno<sup>59,93</sup>, R. Toledo-Moreo<sup>94</sup>, F. Torradeflot<sup>50,92</sup>, I. Tutusaus<sup>11</sup>, E. A. Valentijn<sup>95</sup>, L. Valenziano<sup>3,6</sup>, T. Vassallo<sup>86,27</sup>, A. Veropalumbo<sup>38,41</sup>, Y. Wang<sup>96</sup>, J. Weller<sup>86,64</sup>, G. Zamorani<sup>3</sup>, J. Zoubian<sup>71</sup>, and V. Scottez<sup>31,97</sup>

*(Affiliations can be found after the references)*

Received 5 October 2023 / Accepted 21 December 2023

## ABSTRACT

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 described by $1 + \mathcal{A}_X \sin(\omega_X \Xi_X + 2\pi\phi_X)$, one linearly spaced in $k$ space with $\Xi_{\text{lin}} \equiv k/k_*$ where $k_*$ is $0.05 \text{ Mpc}^{-1}$ and the other logarithmically spaced in $k$ space with $\Xi_{\text{log}} \equiv \ln(k/k_*)$. We note that $\mathcal{A}_X$ is the amplitude of the primordial feature, $\omega_X$ is the dimensionless frequency, and $\phi_X$ is the normalised phase, where $X = \{\text{lin, log}\}$. We provide forecasts from spectroscopic and photometric primary Euclid probes on the standard cosmological parameters $\Omega_{m,0}$, $\Omega_{b,0}$, $h$, $n_s$, and $\sigma_8$, and the primordial feature parameters $\mathcal{A}_X$, $\omega_X$, and $\phi_X$. We focus on the uncertainties of the primordial feature amplitude $\mathcal{A}_X$ and on the capability of Euclid to detect primordial features at a given frequency. We also study a nonlinear density reconstruction method in order to retrieve the oscillatory signals in the primordial power spectrum, which are damped on small scales in the late-time Universe due to cosmic structure formation. Finally, we also include the expected measurements from Euclid's galaxy-clustering bispectrum and from observations of the cosmic microwave background (CMB). We forecast uncertainties in estimated values of the cosmological parameters with a Fisher matrix method applied to spectroscopic galaxy clustering (GC<sub>sp</sub>), weak lensing (WL), photometric galaxy clustering (GC<sub>ph</sub>), the cross correlation (XC) between GC<sub>ph</sub> and WL, the spectroscopic galaxy clustering bispectrum, the CMB temperature and E-mode polarisation, the temperature-polarisation cross correlation, and CMB weak lensing. We consider two sets of specifications for the Euclid probes (pessimistic and optimistic) and three different CMB experiment configurations, that is, Planck, Simons Observatory (SO), and CMB Stage-4 (CMB-S4). We find the following percentage relative errors in the feature amplitude with Euclid primary probes: for the linear (logarithmic) feature model, with a fiducial value of $\mathcal{A}_X = 0.01$, $\omega_X = 10$, and $\phi_X = 0$: 21% (22%) in the pessimistic settings and 18% (18%) in the optimistic settings at a 68.3% confidence level (CL) using GC<sub>sp</sub>+WL+GC<sub>ph</sub>+XC. While the uncertainties on the feature amplitude are strongly dependent on the frequency value when single Euclid probes are considered, we find robust constraints on $\mathcal{A}_X$ from the combination of spectroscopic and photometric measurements over the frequency range of (1, $10^{2.1}$). Due to the inclusion of numerical reconstruction, the GC<sub>sp</sub> bispectrum, SO-like CMB reduces the uncertainty on the primordial feature amplitude by 32%–48%, 50%–65%, and 15%–50%, respectively. Combining all the sources of information explored expected from Euclid in combination with the future SO-like CMB experiment, we forecast $\mathcal{A}_{\text{lin}} \approx 0.010 \pm 0.001$ at a 68.3% CL and $\mathcal{A}_{\text{log}} \approx 0.010 \pm 0.001$ for GC<sub>sp</sub>(PS rec + BS)+WL+GC<sub>ph</sub>+XC+SO-like for both the optimistic and pessimistic settings over the frequency range (1, $10^{2.1}$).

**Key words.** gravitation – gravitational lensing: weak – cosmological parameters – early Universe – large-scale structure of Universe

<footnote>* This paper is published on behalf of the Euclid Consortium.</footnote>
# Euclid: The search for primordial features*

M. Ballardini<sup>1,2,3</sup> ◦, Y. Akrami<sup>4,5</sup>, F. Finelli<sup>3,6</sup>, D. Karagiannis<sup>7</sup>, B. Li<sup>8</sup>, Y. Li<sup>9</sup>, Z. Sakr<sup>10,11,12</sup>, D. Sapone<sup>13</sup>, A. Achúcarro<sup>14,15</sup>, M. Baldi<sup>16,3,17</sup>, N. Bartolo<sup>18,19,20</sup>, G. Cañas-Herrera<sup>21,14</sup>, S. Casas<sup>22</sup>, R. Murgia<sup>23,24</sup>, H.-A. Winther<sup>25</sup>, M. Viel<sup>26,27,28,29</sup>, A. Andrews<sup>3</sup>, J. Jasche<sup>30,31</sup>, G. Lavaux<sup>32</sup> ◦, D. K. Hazra<sup>3,33,34</sup>, D. Paoletti<sup>3,6</sup>, J. Valiviita<sup>35,36</sup>, A. Amara<sup>37</sup>, S. Andreon<sup>38</sup>, N. Auricchio<sup>3</sup>, P. Battaglia<sup>3</sup>, D. Bonino<sup>39</sup>, E. Branchini<sup>40,41</sup>, M. Brescia<sup>42,43,44</sup>, J. Brinchmann<sup>45</sup>, S. Camera<sup>46,47,39</sup>, V. Capobianco<sup>39</sup>, C. Carbone<sup>48</sup>, J. Carretero<sup>49,50</sup>, M. Castellano<sup>51</sup>, S. Cavuoti<sup>43,44</sup>, A. Cimatti<sup>52</sup>, G. Congedo<sup>53</sup>, L. Conversi<sup>54,55</sup>, Y. Copin<sup>56</sup>, L. Corcione<sup>39</sup>, F. Courbin<sup>57</sup>, H. M. Courtois<sup>58</sup>, A. Da Silva<sup>59,60</sup>, H. Degaudenzi<sup>61</sup>, F. Dubath<sup>61</sup>, X. Dupac<sup>55</sup>, M. Farina<sup>62</sup>, S. Farrens<sup>63</sup>, M. Frailis<sup>27</sup>, E. Franceschi<sup>3</sup>, M. Fumana<sup>48</sup>, S. Galeotta<sup>27</sup>, B. Gillis<sup>53</sup>, C. Giocoli<sup>3,17</sup>, A. Grazian<sup>20</sup>, F. Grupp<sup>64,65</sup>, S. V. H. Haugan<sup>25</sup>, W. Holmes<sup>66</sup>, F. Hormuth<sup>67</sup>, A. Hornstrup<sup>68,69</sup>, P. Hudelot<sup>32</sup>, K. Jahnke<sup>70</sup>, S. Kermiche<sup>71</sup>, A. Kiessling<sup>66</sup>, M. Kunz<sup>72</sup>, H. Kurki-Suonio<sup>35,36</sup>, P. B. Lilje<sup>25</sup>, V. Lindholm<sup>35,36</sup>, I. Lloro<sup>73</sup>, E. Maiorano<sup>3</sup>, O. Mansutti<sup>27</sup>, O. Marggraf<sup>74</sup>, N. Martinet<sup>75</sup>, F. Marulli<sup>76,3,17</sup>, R. Massey<sup>8</sup>, E. Medinaceli<sup>3</sup>, S. Mei<sup>77</sup>, Y. Mellier<sup>31,32</sup>, M. Meneghetti<sup>3,17</sup>, E. Merlin<sup>51</sup>, G. Meylan<sup>57</sup>, M. Moresco<sup>76,3</sup>, L. Moscardini<sup>76,3,17</sup>, E. Munari<sup>27</sup>, S.-M. Niemi<sup>21</sup>, C. Padilla<sup>49</sup>, S. Paltani<sup>61</sup>, F. Pasian<sup>27</sup>, K. Pedersen<sup>78</sup>, W. J. Percival<sup>79,80,81</sup>, V. Pettorino<sup>82</sup>, S. Pires<sup>63</sup>, G. Polenta<sup>83</sup>, M. Poncet<sup>84</sup>, L. A. Popa<sup>85</sup>, L. Pozzetti<sup>3</sup>, F. Raison<sup>64</sup>, A. Renzi<sup>18,19</sup>, J. Rhodes<sup>66</sup>, G. Riccio<sup>43</sup>, E. Romelli<sup>27</sup>, M. Roncarelli<sup>3</sup>, R. Saglia<sup>86,64</sup>, B. Sartoris<sup>86,27</sup>, T. Schrabback<sup>87</sup>, A. Secroun<sup>71</sup>, G. Seidel<sup>70</sup>, S. Serrano<sup>88,89,90</sup>, C. Sirignano<sup>18,19</sup>, G. Sirri<sup>17</sup>, L. Stanco<sup>19</sup>, J.-L. Starck<sup>91</sup>, C. Surace<sup>75</sup>, P. Tallada-Crespi<sup>92,50</sup>, A. N. Taylor<sup>53</sup>, I. Tereno<sup>59,93</sup>, R. Toledo-Moreo<sup>94</sup>, F. Torradeflot<sup>50,92</sup>, I. Tutusaus<sup>11</sup>, E. A. Valentijn<sup>95</sup>, L. Valenziano<sup>3,6</sup>, T. Vassallo<sup>86,27</sup>, A. Veropalumbo<sup>38,41</sup>, Y. Wang<sup>96</sup>, J. Weller<sup>86,64</sup>, G. Zamorani<sup>3</sup>, J. Zoubian<sup>71</sup>, and V. Scottez<sup>31,97</sup>

*(Affiliations can be found after the references)*

Received 5 October 2023 / Accepted 21 December 2023

## ABSTRACT

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 described by $1 + \mathcal{A}_X \sin(\omega_X \Xi_X + 2\pi\phi_X)$, one linearly spaced in $k$ space with $\Xi_{\text{lin}} \equiv k/k_*$ where $k_*$ is $0.05 \text{ Mpc}^{-1}$ and the other logarithmically spaced in $k$ space with $\Xi_{\text{log}} \equiv \ln(k/k_*)$. We note that $\mathcal{A}_X$ is the amplitude of the primordial feature, $\omega_X$ is the dimensionless frequency, and $\phi_X$ is the normalised phase, where $X = \{\text{lin, log}\}$. We provide forecasts from spectroscopic and photometric primary Euclid probes on the standard cosmological parameters $\Omega_{m,0}$, $\Omega_{b,0}$, $h$, $n_s$, and $\sigma_8$, and the primordial feature parameters $\mathcal{A}_X$, $\omega_X$, and $\phi_X$. We focus on the uncertainties of the primordial feature amplitude $\mathcal{A}_X$ and on the capability of Euclid to detect primordial features at a given frequency. We also study a nonlinear density reconstruction method in order to retrieve the oscillatory signals in the primordial power spectrum, which are damped on small scales in the late-time Universe due to cosmic structure formation. Finally, we also include the expected measurements from Euclid's galaxy-clustering bispectrum and from observations of the cosmic microwave background (CMB). We forecast uncertainties in estimated values of the cosmological parameters with a Fisher matrix method applied to spectroscopic galaxy clustering (GC<sub>sp</sub>), weak lensing (WL), photometric galaxy clustering (GC<sub>ph</sub>), the cross correlation (XC) between GC<sub>ph</sub> and WL, the spectroscopic galaxy clustering bispectrum, the CMB temperature and E-mode polarisation, the temperature-polarisation cross correlation, and CMB weak lensing. We consider two sets of specifications for the Euclid probes (pessimistic and optimistic) and three different CMB experiment configurations, that is, Planck, Simons Observatory (SO), and CMB Stage-4 (CMB-S4). We find the following percentage relative errors in the feature amplitude with Euclid primary probes: for the linear (logarithmic) feature model, with a fiducial value of $\mathcal{A}_X = 0.01$, $\omega_X = 10$, and $\phi_X = 0$: 21% (22%) in the pessimistic settings and 18% (18%) in the optimistic settings at a 68.3% confidence level (CL) using GC<sub>sp</sub>+WL+GC<sub>ph</sub>+XC. While the uncertainties on the feature amplitude are strongly dependent on the frequency value when single Euclid probes are considered, we find robust constraints on $\mathcal{A}_X$ from the combination of spectroscopic and photometric measurements over the frequency range of (1, 10<sup>2.1</sup>). Due to the inclusion of numerical reconstruction, the GC<sub>sp</sub> bispectrum, SO-like CMB reduces the uncertainty on the primordial feature amplitude by 32%–48%, 50%–65%, and 15%–50%, respectively. Combining all the sources of information explored expected from Euclid in combination with the future SO-like CMB experiment, we forecast $\mathcal{A}_{\text{lin}} \approx 0.010 \pm 0.001$ at a 68.3% CL and $\mathcal{A}_{\text{log}} \approx 0.010 \pm 0.001$ for GC<sub>sp</sub>(PS rec + BS)+WL+GC<sub>ph</sub>+XC+SO-like for both the optimistic and pessimistic settings over the frequency range (1, 10<sup>2.1</sup>).

**Key words.** gravitation – gravitational lensing: weak – cosmological parameters – early Universe – large-scale structure of Universe

<footnote>* This paper is published on behalf of the Euclid Consortium.</footnote>
# 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 $\Omega_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,a; 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^{-1}) < 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 $\ell \sim 20-40$ and an oscillatory pattern at the intermediate scales $\ell \sim 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_0$ and $S_8$ 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 comple-

mentary 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, $P_{R,0}(k)$, with the following parametrisation

$$P_{R,0}(k) = A_s \left(\frac{k}{k_*}\right)^{n_s-1}, \qquad (1)$$

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 \text{ Mpc}^{-1}$.

The simplest class of models that can produce such a spectrum is known as canonical single-field slow-roll inflation. In
# 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 $\Omega_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,a; 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 \text{ Mpc}^{-1}) < 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 $\ell \sim 20-40$ and an oscillatory pattern at the intermediate scales $\ell \sim 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_0$ and $S_8$ 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 comple-

mentary 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, $P_{R,0}(k)$, with the following parametrisation

$$P_{R,0}(k) = A_s \left(\frac{k}{k_*}\right)^{n_s-1}, \qquad (1)$$

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 \text{ 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<fnref n="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,

$$S = \int d^4 x a^3 M_{\text{Pl}}^2 \epsilon_1 \left[ \frac{\dot{R}^2}{c_s^2} - \frac{(\partial_i R)^2}{a^2} \right] + \dots, \quad (2)$$

where $a(t)$ is the scale factor, $c_s(t)$ is the speed of sound for curvature perturbations, and $\epsilon_1(t) \equiv -\dot{H}/H^2$ is the first slow-roll parameter; $H \equiv \dot{a}/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 ($\epsilon_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 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 $\epsilon_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 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, $\epsilon_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 $\epsilon_1$ 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\alpha$ 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

<footnote n="1">¹ These models are based additionally on other simplifying assumptions about the inflaton field, such as canonical kinetic terms, minimal coupling to gravity and the Bunch-Davies quantum vacuum for the perturbations.</footnote>
these models a single scalar field – the inflaton – is responsible for both the primordial perturbations and the background evolution during inflation<fnref n="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,

$$S = \int d^4 x a^3 M_{\text{Pl}}^2 \epsilon_1 \left[ \frac{\dot{R}^2}{c_s^2} - \frac{(\partial_i R)^2}{a^2} \right] + \dots, \quad (2)$$

where $a(t)$ is the scale factor, $c_s(t)$ is the speed of sound for curvature perturbations, and $\epsilon_1(t) \equiv -\dot{H}/H^2$ is the first slow-roll parameter; $H \equiv \dot{a}/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 ($\epsilon_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 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 $\epsilon_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 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, $\epsilon_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 $\epsilon_1$ 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\alpha$ 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

<footnote n="1">¹ These models are based additionally on other simplifying assumptions about the inflaton field, such as canonical kinetic terms, minimal coupling to gravity and the Bunch-Davies quantum vacuum for the perturbations.</footnote>
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

$$P_{\text{obs}}(k_{\text{ref}}, \mu_{\text{ref}}; z) = \frac{1}{q_{\perp}^2 q_{\parallel}} \frac{1}{1 + k^2 \mu^2 f^2(z) \sigma_{\text{p}}^2(z)} \times \left[b(z)\sigma_8(z) + f(z)\sigma_8(z)\mu^2\right]^2 \times \frac{P^{\text{IR res, LO}}(k, \mu; z)}{\sigma_8^2(z)} F_z(k, \mu; z) + P_s(z), \quad (3)$$

where $\mu$ is the cosine of the angle of the wave mode with respect to the line of sight pointing into the direction $\hat{r}$. Here $b(z)$ is the linear clustering bias, $f(z)$ is the growth rate, and $\sigma_8(z)$ is the root mean square linearly evolved density fluctuations in spheres of $8 h^{-1} \text{ Mpc}^2$.

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

$$q_{\perp}(z) = \frac{D_A(z)}{D_{A, \text{ref}}(z)} \quad \text{and} \quad q_{\parallel}(z) = \frac{H_{\text{ref}}(z)}{H(z)}, \quad (4)$$

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

$$k_{\perp} = \frac{k_{\perp, \text{ref}}}{q_{\perp}} \quad \text{and} \quad k_{\parallel} = \frac{k_{\parallel, \text{ref}}}{q_{\parallel}}. \quad (5)$$

Using Eqs. (4) and (5) we can convert the known reference cosmology ($k_{\text{ref}}$ and $\mu_{\text{ref}}$) to the true unknown cosmology ($k$ and $\mu$),

$$k(k_{\text{ref}}, \mu_{\text{ref}}) = k_{\text{ref}} \alpha(\mu_{\text{ref}}), \quad (6)$$

$$\mu(\mu_{\text{ref}}) = \frac{\mu_{\text{ref}}}{\alpha(\mu_{\text{ref}})q_{\parallel}}, \quad (7)$$

where

$$\alpha(\mu) = \left[ \frac{\mu^2}{q_{\parallel}^2} + \frac{(1-\mu^2)}{q_{\perp}^2} \right]^{1/2}. \quad (8)$$

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

$$P_{\text{obs}}(k_{\text{ref}}, \mu_{\text{ref}}; z) = \frac{1}{q_{\perp}^2 q_{\parallel}^2} P_g(k(k_{\text{ref}}, \mu(\mu_{\text{ref}}); z)). \quad (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^{\text{IR res, LO}}$ in Sect. 4.2. The pairwise velocity dispersion, $\sigma_p(z)$, is evaluated from the linear matter power spectrum as

$$\sigma_p^2(z) = \frac{1}{6\pi^2} \int_0^{k_{\text{max}}} dk P_{\text{lin}}(k, z). \quad (10)$$

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

$$F_z(k, \mu; z) = e^{-k^2 \mu^2 \sigma_r^2(z)}, \quad (11)$$

where $\sigma_r^2(z) = c(1+z)\sigma_{0,z}/H(z)$ with $\sigma_{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 (GCsp) observable for one redshift bin $z_i$ is

$$F_{\alpha\beta}(z_i) = \frac{1}{8\pi^2} \int_{-1}^{1} d\mu \int_{k_{\text{min}}}^{k_{\text{max}}} k^2 dk V_{\text{eff}}(z_i, k) \times \frac{\partial \ln P_{\text{obs}}(k, \mu; z_i)}{\partial \theta_{\alpha}} \frac{\partial \ln P_{\text{obs}}(k, \mu; z_i)}{\partial \theta_{\beta}}, \quad (12)$$

where the derivatives are evaluated at the parameter values of the fiducial model and $V_{\text{eff}}$ is the effective volume of the survey, given by

$$V_{\text{eff}}(k, \mu; z) = V_s \left[ \frac{n(z) P_{\text{obs}}(k, \mu; z)}{1 + n(z) P_{\text{obs}}(k, \mu; z)} \right]^2, \quad (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

$$\text{Cov}(k, k') \approx \frac{2(2\pi)^3}{V_{\text{eff}}(k, \mu; z)} P_{\text{obs}}^2(k, \mu; z) \delta^D(k-k'), \quad (14)$$

where $\delta^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 $\Delta k$ with a top-hat window function as in Huang et al. (2012), Ballardini et al. (2016),

$$\hat{P}_{\text{obs}}(k, \mu; z) = \frac{1}{\Delta k} \int_{k-\Delta k/2}^{k+\Delta k/2} dk' P_{\text{obs}}(k', \mu; z). \quad (15)$$

The size of the k-bin $\Delta k$ should be equal or larger than the effective fundamental frequency defined as $k_{\text{eff}} \equiv 2\pi/V^{1/3}$ for an ideal cubic volume in order to guarantee uncorrelated measurements of the observed power spectrum. The values of $k_{\text{eff}}$ span in the

<footnote n="2">² h is the dimensionless Hubble parameter defined as $H_0 = 100 h \text{ km s}^{-1} \text{ Mpc}^{-1}$.</footnote>
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

$$P_{\text{obs}}(k_{\text{ref}}, \mu_{\text{ref}}; z) = \frac{1}{q_{\perp}^2 q_{\parallel}} \frac{1}{1 + k^2 \mu^2 f^2(z) \sigma_{\text{p}}^2(z)} \times \left[ b(z) \sigma_8(z) + f(z) \sigma_8(z) \mu^2 \right]^2 \times \frac{P^{\text{IR res, LO}}(k, \mu; z)}{\sigma_8^2(z)} F_z(k, \mu; z) + P_s(z), \quad (3)$$

where $\mu$ is the cosine of the angle of the wave mode with respect to the line of sight pointing into the direction $\hat{r}$. Here $b(z)$ is the linear clustering bias, $f(z)$ is the growth rate, and $\sigma_8(z)$ is the root mean square linearly evolved density fluctuations in spheres of $8 h^{-1} \text{ Mpc}^2$.

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

$$q_{\perp}(z) = \frac{D_A(z)}{D_{A, \text{ref}}(z)} \quad \text{and} \quad q_{\parallel}(z) = \frac{H_{\text{ref}}(z)}{H(z)}, \quad (4)$$

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

$$k_{\perp} = \frac{k_{\perp, \text{ref}}}{q_{\perp}} \quad \text{and} \quad k_{\parallel} = \frac{k_{\parallel, \text{ref}}}{q_{\parallel}}. \quad (5)$$

Using Eqs. (4) and (5) we can convert the known reference cosmology ($k_{\text{ref}}$ and $\mu_{\text{ref}}$) to the true unknown cosmology ($k$ and $\mu$),

$$k(k_{\text{ref}}, \mu_{\text{ref}}) = k_{\text{ref}} \alpha(\mu_{\text{ref}}), \quad (6)$$

$$\mu(\mu_{\text{ref}}) = \frac{\mu_{\text{ref}}}{\alpha(\mu_{\text{ref}})q_{\parallel}}, \quad (7)$$

where

$$\alpha(\mu) = \left[ \frac{\mu^2}{q_{\parallel}^2} + \frac{(1 - \mu^2)}{q_{\perp}^2} \right]^{1/2}. \quad (8)$$

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

$$P_{\text{obs}}(k_{\text{ref}}, \mu_{\text{ref}}; z) = \frac{1}{q_{\perp}^2 q_{\parallel}^2} P_g(k(k_{\text{ref}}, \mu(\mu_{\text{ref}}); z)). \quad (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^{\text{IR res, LO}}$ in Sect. 4.2. The pairwise velocity dispersion, $\sigma_p(z)$, is evaluated from the linear matter power spectrum as

$$\sigma_p^2(z) = \frac{1}{6\pi^2} \int_0^{k_{\text{max}}} dk P_{\text{lin}}(k, z). \quad (10)$$

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

$$F_z(k, \mu; z) = e^{-k^2 \mu^2 \sigma_r^2(z)}, \quad (11)$$

where $\sigma_r^2(z) = c(1+z)\sigma_{0,z}/H(z)$ with $\sigma_{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 (GCsp) observable for one redshift bin $z_i$ is

$$F_{\alpha\beta}(z_i) = \frac{1}{8\pi^2} \int_{-1}^{1} d\mu \int_{k_{\text{min}}}^{k_{\text{max}}} k^2 dk V_{\text{eff}}(z_i, k) \times \frac{\partial \ln P_{\text{obs}}(k, \mu; z_i)}{\partial \theta_{\alpha}} \frac{\partial \ln P_{\text{obs}}(k, \mu; z_i)}{\partial \theta_{\beta}}, \quad (12)$$

where the derivatives are evaluated at the parameter values of the fiducial model and $V_{\text{eff}}$ is the effective volume of the survey, given by

$$V_{\text{eff}}(k, \mu; z) = V_s \left[ \frac{n(z) P_{\text{obs}}(k, \mu; z)}{1 + n(z) P_{\text{obs}}(k, \mu; z)} \right]^2, \quad (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

$$\text{Cov}(k, k') \approx \frac{2(2\pi)^3}{V_{\text{eff}}(k, \mu; z)} P_{\text{obs}}^2(k, \mu; z) \delta^{\text{D}}(k - k'), \quad (14)$$

where $\delta^{\text{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 $\Delta k$ with a top-hat window function as in Huang et al. (2012), Ballardini et al. (2016),

$$\hat{P}_{\text{obs}}(k, \mu; z) = \frac{1}{\Delta k} \int_{k-\Delta k/2}^{k+\Delta k/2} dk' P_{\text{obs}}(k', \mu; z). \quad (15)$$

The size of the k-bin $\Delta k$ should be equal or larger than the effective fundamental frequency defined as $k_{\text{eff}} \equiv 2\pi/V^{1/3}$ for an ideal cubic volume in order to guarantee uncorrelated measurements of the observed power spectrum. The values of $k_{\text{eff}}$ span in the

<footnote n="2">² h is the dimensionless Hubble parameter defined as $H_0 = 100 h \text{ km s}^{-1} \text{ Mpc}^{-1}$.</footnote>
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ᵢ as

$$F_{\alpha\beta}(z_i) = \frac{1}{8\pi^2} \int_{-1}^{1} d\mu \sum_{k=k_{\min}}^{k_{\max}} k^2 \Delta k V_{\text{eff}}(z_i, k) \\ \times \frac{\partial \ln \hat{P}_{\text{obs}}(k, \mu; z_i) \partial \ln \hat{P}_{\text{obs}}(k, \mu; z_i)}{\partial \theta_{\alpha} \partial \theta_{\beta}}. \quad (16)$$

The total spectroscopic Fisher matrix is then calculated by summing over the redshift bins as $F_{\alpha\beta}^{\text{sp}} = \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<sub>sp</sub> Fisher matrix over two redshift-dependent parameters $b\sigma_8(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_{i,j,l}^{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_{i,j,l}^{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

$$C_{i,j,l}^{XY} = c \int_{z_{\min}}^{z_{\max}} dz \frac{W_i^X(z) W_j^Y(z)}{H(z) r^2(z)} P_{\text{NL}}(k_l, z), \quad (17)$$

with $k_l = (l + 1/2)/r(z)$, $r(z)$ the comoving distance to redshift z, and $P_{\text{NL}}(k_l, z)$ the nonlinear power spectrum of matter density fluctuations at wave number $k_l$ and redshift z. The GC<sub>ph</sub> and WL window functions read

$$W_i^G(z) = b_i(z) n_i(z) \frac{H(z)}{c}, \quad (18)$$

$$W_i^L(z) = \frac{3}{2} \Omega_{m,0} \left(\frac{H_0}{c}\right)^2 (1+z) r(z) \int_z^{z_{\max}} dz' n_i(z') \left[1 - \frac{r(z)}{r(z')}\right] + W_i^{\text{IA}}(z), \quad (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

$$n_i(z) = \frac{\int_{z_i^-}^{z_i^+} dz' n(z) p_{\text{ph}}(z'|z)}{\int_{z_{\min}}^{z_{\max}} dz \int_{z_i^-}^{z_i^+} dz' n(z) p_{\text{ph}}(z'|z)}, \quad (20)$$

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_{\text{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

$$W_i^{\text{IA}}(z) = - \frac{\mathcal{A}_{\text{IA}} C_{\text{IA}} \Omega_{\text{m},0} \mathcal{F}_{\text{IA}}(z)}{D(z)} n_i(z) \frac{H(z)}{c}, \quad (21)$$

where

$$\mathcal{F}_{\text{IA}}(z) = (1+z)^{\eta_{\text{IA}}} \left[ \frac{\langle L \rangle(z)}{L_{\star}(z)} \right]^{\beta_{\text{IA}}}, \quad (22)$$

with $\langle L \rangle(z)$ the redshift-dependent mean, and $L_{\star}(z)$ the characteristic luminosity of source galaxies as computed from the luminosity function. $\mathcal{A}_{\text{IA}}$, $\beta_{\text{IA}}$ and $\eta_{\text{IA}}$ are the nuisance parameters of the model, and $C_{\text{IA}}$ is a constant accounting for dimensional units.

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

$$\text{Cov}[C_{i,j,l}^{\text{AB}}, C_{mn,l'}^{\text{CD}}] = \frac{\delta_{l'l}^{K}}{(2\ell + 1) f_{\text{sky}} \Delta \ell} \\ \times \left[ \left[ C_{im,l}^{\text{AC}} + N_{im,l}^{\text{AC}} \right] \left[ C_{jn,l'}^{\text{BD}} + N_{jn,l'}^{\text{BD}} \right] \right. \\ \left. + \left[ C_{in,l}^{\text{AD}} + N_{in,l}^{\text{AD}} \right] \left[ C_{jm,l'}^{\text{BC}} + N_{jm,l'}^{\text{BC}} \right] \right], \quad (23)$$

where the upper- and lower-case Latin indices run over L and G (all tomographic bins), $\delta_{l'l}^{K}$ is the Kronecker delta coming from the lack of correlation between different multipoles $(\ell, l')$, $f_{\text{sky}}$ is the survey's sky fraction, and $\Delta\ell$ denotes the width of the logarithmic equi-spaced multipole bins. We consider a white noise,

$$N_{i,j,l}^{\text{LL}} = \frac{\delta_{ij}^{K}}{\bar{n}_{i}} \sigma_{\epsilon}^{2}, \quad N_{i,j,l}^{\text{GG}} = \frac{\delta_{ij}^{K}}{\bar{n}_{i}}, \quad N_{i,j,l}^{\text{GL}} = 0, \quad (24)$$

where $\sigma_{\epsilon}$ is the variance of observed ellipticities.

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

$$F_{\alpha\beta}^{\text{ph}} = \sum_{\ell=\ell_{\min}}^{\ell_{\max}} \sum_{i,j,mn} \frac{\partial C_{i,j,l}^{\text{AB}}}{\partial \theta_{\alpha}} \text{Cov}^{-1} [C_{i,j,l}^{\text{AB}}, C_{mn,l'}^{\text{CD}}] \frac{\partial C_{mn,l'}^{\text{CD}}}{\partial \theta_{\beta}}. \quad (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 $\Delta z = 0.2$ for the first three bins and $\Delta z = 0.3$ for the last bin. The $H\alpha$ 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} $\times 10^{-4}$.

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

$$n(z) \propto \left(\frac{z}{z_0}\right)^2 \exp\left[-\left(\frac{z}{z_0}\right)^{3/2}\right], \quad (26)$$

with $z_0 = 0.9/\sqrt{2}$ and the normalisation set by the requirement that the surface density of galaxies is $\bar{n}_g = 30 \text{ 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
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ᵢ as

$$F_{\alpha\beta}(z_i) = \frac{1}{8\pi^2} \int_{-1}^{1} d\mu \sum_{k=k_{\min}}^{k_{\max}} k^2 \Delta k V_{\text{eff}}(z_i, k) \\ \times \frac{\partial \ln \hat{P}_{\text{obs}}(k, \mu; z_i) \partial \ln \hat{P}_{\text{obs}}(k, \mu; z_i)}{\partial \theta_{\alpha} \partial \theta_{\beta}}. \quad (16)$$

The total spectroscopic Fisher matrix is then calculated by summing over the redshift bins as $F_{\alpha\beta}^{\text{sp}} = \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<sub>sp</sub> Fisher matrix over two redshift-dependent parameters $b\sigma_8(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}^{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}^{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

$$C_{ij,\ell}^{XY} = c \int_{z_{\min}}^{z_{\max}} dz \frac{W_i^X(z) W_j^Y(z)}{H(z) r^2(z)} P_{\text{NL}}(k_\ell, z), \quad (17)$$

with $k_\ell = (\ell + 1/2)/r(z)$, $r(z)$ the comoving distance to redshift z, and $P_{\text{NL}}(k_\ell, z)$ the nonlinear power spectrum of matter density fluctuations at wave number $k_\ell$ and redshift z. The GC<sub>ph</sub> and WL window functions read

$$W_i^G(z) = b_i(z) n_i(z) \frac{H(z)}{c}, \quad (18)$$

$$W_i^L(z) = \frac{3}{2} \Omega_{m,0} \left(\frac{H_0}{c}\right)^2 (1+z) r(z) \int_z^{z_{\max}} dz' n_i(z') \left[1 - \frac{r(z)}{r(z')}\right] + W_i^{\text{IA}}(z), \quad (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

$$n_i(z) = \frac{\int_{z_i^-}^{z_i^+} dz' n(z) p_{\text{ph}}(z'|z)}{\int_{z_{\min}}^{z_{\max}} dz \int_{z_i^-}^{z_i^+} dz' n(z) p_{\text{ph}}(z'|z)}, \quad (20)$$

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_{\text{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

$$W_i^{\text{IA}}(z) = - \frac{\mathcal{A}_{\text{IA}} C_{\text{IA}} \Omega_{m,0} \mathcal{F}_{\text{IA}}(z)}{D(z)} n_i(z) \frac{H(z)}{c}, \quad (21)$$

where

$$\mathcal{F}_{\text{IA}}(z) = (1+z)^{\eta_{\text{IA}}} \left[ \frac{\langle L \rangle(z)}{L_{\star}(z)} \right]^{\beta_{\text{IA}}}, \quad (22)$$

with $\langle L \rangle(z)$ the redshift-dependent mean, and $L_{\star}(z)$ the characteristic luminosity of source galaxies as computed from the luminosity function. $\mathcal{A}_{\text{IA}}$, $\beta_{\text{IA}}$ and $\eta_{\text{IA}}$ are the nuisance parameters of the model, and $C_{\text{IA}}$ is a constant accounting for dimensional units.

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

$$\text{Cov}[C_{ij,\ell}^{\text{AB}}, C_{mn,\ell'}^{\text{CD}}] = \frac{\delta_{\ell\ell'}^{K}}{(2\ell + 1) f_{\text{sky}} \Delta\ell} \\ \times \left\{ \left[C_{im,\ell}^{\text{AC}} + N_{im,\ell}^{\text{AC}}\right] \left[C_{jn,\ell'}^{\text{BD}} + N_{jn,\ell'}^{\text{BD}}\right] \right. \\ \left. + \left[C_{in,\ell}^{\text{AD}} + N_{in,\ell}^{\text{AD}}\right] \left[C_{jm,\ell'}^{\text{BC}} + N_{jm,\ell'}^{\text{BC}}\right] \right\}, \quad (23)$$

where the upper- and lower-case Latin indices run over L and G (all tomographic bins), $\delta_{\ell\ell'}^K$ is the Kronecker delta coming from the lack of correlation between different multipoles $(\ell, \ell')$, $f_{\text{sky}}$ is the survey's sky fraction, and $\Delta\ell$ denotes the width of the logarithmic equi-spaced multipole bins. We consider a white noise,

$$N_{ij,\ell}^{\text{LL}} = \frac{\delta_{ij}^{K}}{\bar{n}_{i}} \sigma_{\epsilon}^{2}, \quad N_{ij,\ell}^{\text{GG}} = \frac{\delta_{ij}^{K}}{\bar{n}_{i}}, \quad N_{ij,\ell}^{\text{GL}} = 0, \quad (24)$$

where $\sigma_\epsilon$ is the variance of observed ellipticities.

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

$$F_{\alpha\beta}^{\text{ph}} = \sum_{\ell=\ell_{\min}}^{\ell_{\max}} \sum_{i,j,mn} \frac{\partial C_{ij,\ell}^{\text{AB}}}{\partial \theta_{\alpha}} \text{Cov}^{-1} [C_{ij,\ell}^{\text{AB}}, C_{mn,\ell}^{\text{CD}}] \frac{\partial C_{mn,\ell}^{\text{CD}}}{\partial \theta_{\beta}}. \quad (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 $\Delta z = 0.2$ for the first three bins and $\Delta z = 0.3$ for the last bin. The $H\alpha$ 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} $\times 10^{-4}$.

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

$$n(z) \propto \left(\frac{z}{z_0}\right)^2 \exp\left[-\left(\frac{z}{z_0}\right)^{3/2}\right], \quad (26)$$

with $z_0 = 0.9/\sqrt{2}$ and the normalisation set by the requirement that the surface density of galaxies is $\bar{n}_g = 30 \text{ 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 + \bar{z}_i}$, where $\bar{z}_i$ is the bin centre. We consider 100 logarithmic equi-spaced multipole bins with $\sigma_\epsilon = 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} \approx 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

$$ \mathcal{P}_{\mathcal{R}}(k) = \mathcal{P}_{\mathcal{R},0}(k) [1 + \delta P^X(k)], \quad (27) $$

where $\mathcal{P}_{\mathcal{R},0}(k)$ is the standard power-law PPS of the comoving curvature perturbations $\mathcal{R}$ on superhorizon scales, given in Eq. (1). We study the following templates with superimposed oscillations on the PPS

$$ \delta P^X(k) = \mathcal{A}_X \sin(\omega_X \Xi_X + 2\pi\phi_X), \quad (28) $$

where $X = \{\text{lin, log}\}$ and $\Xi_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): $\omega_b = 0.02237$, $\omega_c = 0.1200$, $H_0 = 67.36 \text{ km s}^{-1} \text{ Mpc}^{-1}$, $\ln(10^10 A_s) = 3.044$, $n_s = 0.9649$, $\sigma_8 = 0.8107$, and $m_\nu = 0.06 \text{ eV}$ (with one massive neutrino and two massless ones). In what follows we specify the fiducial values of the model parameters.

### 1. Linear oscillations:

$$ \Theta_{\text{lin}} = \{\mathcal{A}_{\text{lin}} = 0.01, \omega_{\text{lin}} = 10, \phi_{\text{lin}} = 0\}. \quad (29) $$

### 2. Logarithmic oscillations:

$$ \Theta_{\text{log}} = \{\mathcal{A}_{\text{log}} = 0.01, \omega_{\text{log}} = 10, \phi_{\text{log}} = 0\}. \quad (30) $$

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

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

The linearly propagated matter power spectrum is given by

$$ P_{\text{lin}}(k, z) = M^2(k, z) P_{\Phi}(k), \quad (31) $$

where $M(k, z) = 2k^2 c^2 T(k) D(z) / (3\Omega_m H_0^2)$, with $T(k)$ the matter transfer function normalised to unity at large scales (i.e. $k \to 0$) and $c$ the speed of light. The gravitational potential power spectrum will then be derived from the two oscillatory models considered: $P_\Phi = 9/25(2\pi^2)\mathcal{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

<figure><img src="image_2.png" /><figcaption>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).</figcaption></figure>

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,

$$ P_{\text{lin}}(k, z) = D^2(z) [P_{\text{nw}}(k) + P_{\text{w}}(k)], \quad (32) $$

where $P_{\text{nw}}$ is the no-wiggle power spectrum, and the oscillatory part $P_{\text{w}}$ describes both the BAO feature and the primordial oscillating feature as

$$ P_{\text{w}}(k) = P_{\text{nw}} \left[ \delta P_{\text{w}}^{\text{BAO}}(k) + \delta P_{\text{w}}^{\text{X}}(k) + \delta P_{\text{w}}^{\text{BAO}}(k) \delta P_{\text{w}}^{\text{X}}(k) \right]. \quad (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 $A_{\text{BAO}} \times \mathcal{A}_{\text{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)

$$ P^{\text{IR res, LO+NLO}}(k, z) = D^2(z) P_{\text{nw}}(k) \times \left\{ \begin{aligned} & 1 + \left[1 + k^2 D^2(z) \Sigma_{\text{BAO}}^2\right] e^{-k^2 D^2(z) \Sigma_{\text{BAO}}^2} \delta P_{\text{w}}^{\text{BAO}}(k) \\ & + \left[1 + k^2 D^2(z) \Sigma_{\text{lin}}^2\right] e^{-k^2 D^2(z) \Sigma_{\text{lin}}^2} \delta P_{\text{w}}^{\text{lin}}(k) \end{aligned} \right\} + D^4(z) P^{1-\text{loop}} \left[ P^{\text{IR res, LO}}(k) \right], \quad (34) $$

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

$$ P^{\text{IR res, LO+NLO}}(k, z) = D^2(z) P_{\text{nw}}(k) \times \left\{ \begin{aligned} & 1 + \left[1 + k^2 D^2(z) \Sigma_{\text{BAO}}^2\right] e^{-k^2 D^2(z) \Sigma_{\text{BAO}}^2} \delta P_{\text{w}}^{\text{BAO}}(k) \\ & + \left[1 + k^2 D^2(z) \Sigma_{\text{log}}^2\right] e^{-k^2 D^2(z) \Sigma_{\text{log}}^2} \left[ \cos(k^2 D^2(z) \hat{\Sigma}_{\text{log}}^2) \delta P_{\text{w}}^{\text{log}}(k) \right. \\ & \qquad \left. - \sin(k^2 D^2(z) \hat{\Sigma}_{\text{log}}^2) \frac{d}{\omega_{\text{log}}} \frac{d}{\ln k} \right] \\ & + D^4(z) P^{1-\text{loop}} \left[ P^{\text{IR res, LO}}(k) \right]. \end{aligned} \right\} \quad (35) $$

$P^{1-\text{loop}}$ is the standard one-loop result, but computed with the leading-order (LO) IR resummed power spectrum. We take the
each redshift bin, with fiducial values $b_i = \sqrt{1 + \bar{z}_i}$, where $\bar{z}_i$ is the bin centre. We consider 100 logarithmic equi-spaced multipole bins with $\sigma_\epsilon = 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} \approx 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

$$ \mathcal{P}_R(k) = \mathcal{P}_{R,0}(k) [1 + \delta P^X(k)], \quad (27) $$

where $\mathcal{P}_{R,0}(k)$ is the standard power-law PPS of the comoving curvature perturbations $\mathcal{R}$ on superhorizon scales, given in Eq. (1). We study the following templates with superimposed oscillations on the PPS

$$ \delta P^X(k) = \mathcal{A}_X \sin(\omega_X \Xi_X + 2\pi\phi_X), \quad (28) $$

where $X = \{\text{lin, log}\}$ and $\Xi_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): $\omega_b = 0.02237$, $\omega_c = 0.1200$, $H_0 = 67.36 \text{ km s}^{-1} \text{ Mpc}^{-1}$, $\ln(10^10 A_s) = 3.044$, $n_s = 0.9649$, $\sigma_8 = 0.8107$, and $m_\nu = 0.06 \text{ eV}$ (with one massive neutrino and two massless ones). In what follows we specify the fiducial values of the model parameters.

### 1. Linear oscillations:

$$ \Theta_{\text{lin}} = \{\mathcal{A}_{\text{lin}} = 0.01, \omega_{\text{lin}} = 10, \phi_{\text{lin}} = 0\}. \quad (29) $$

### 2. Logarithmic oscillations:

$$ \Theta_{\text{log}} = \{\mathcal{A}_{\text{log}} = 0.01, \omega_{\text{log}} = 10, \phi_{\text{log}} = 0\}. \quad (30) $$

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

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

The linearly propagated matter power spectrum is given by

$$ P_{\text{lin}}(k, z) = M^2(k, z) P_{\Phi}(k), \quad (31) $$

where $M(k, z) = 2k^2 c^2 T(k) D(z) / (3\Omega_m H_0^2)$, with $T(k)$ the matter transfer function normalised to unity at large scales (i.e. $k \to 0$) and $c$ the speed of light. The gravitational potential power spectrum will then be derived from the two oscillatory models considered: $P_\Phi = 9/25(2\pi^2)\mathcal{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

<figure><img src="image_2.png" /><figcaption>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).</figcaption></figure>

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,

$$ P_{\text{lin}}(k, z) = D^2(z) [P_{\text{nw}}(k) + P_{\text{w}}(k)], \quad (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

$$ P_w(k) = P_{nw} \left[ \delta P_w^{\text{BAO}}(k) + \delta P_w^X(k) + \delta P_w^{\text{BAO}}(k) \delta P_w^X(k) \right]. \quad (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 $A_{\text{BAO}} \times \mathcal{A}_{\text{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)

$$ P^{\text{IR res, LO+NLO}}(k, z) = D^2(z) P_{\text{nw}}(k) \times \left\{ \begin{aligned} & 1 + \left[1 + k^2 D^2(z) \Sigma_{\text{BAO}}^2\right] e^{-k^2 D^2(z) \Sigma_{\text{BAO}}^2} \delta P_w^{\text{BAO}}(k) \\ & + \left[1 + k^2 D^2(z) \Sigma_{\text{lin}}^2\right] e^{-k^2 D^2(z) \Sigma_{\text{lin}}^2} \delta P_w^{\text{lin}}(k) \end{aligned} \right\} + D^4(z) P^{1-\text{loop}} \left[ P^{\text{IR res, LO}}(k) \right], \quad (34) $$

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

$$ P^{\text{IR res, LO+NLO}}(k, z) = D^2(z) P_{\text{nw}}(k) \times \left\{ \begin{aligned} & 1 + \left[1 + k^2 D^2(z) \Sigma_{\text{BAO}}^2\right] e^{-k^2 D^2(z) \Sigma_{\text{BAO}}^2} \delta P_w^{\text{BAO}}(k) \\ & + \left[1 + k^2 D^2(z) \Sigma_{\text{log}}^2\right] e^{-k^2 D^2(z) \Sigma_{\text{log}}^2} \left[ \cos(k^2 D^2(z) \hat{\Sigma}_{\text{log}}^2) \delta P_w^{\text{log}}(k) \right. \\ & \qquad \left. - \sin(k^2 D^2(z) \hat{\Sigma}_{\text{log}}^2) \frac{d}{\omega_{\text{log}}} \frac{d}{\ln k} \right] \\ & + D^4(z) P^{1-\text{loop}} \left[ P^{\text{IR res, LO}}(k) \right]. \end{aligned} \right\} \quad (35) $$

$P^{1-\text{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-\text{loop}} = P_{22} + 2P_{13}$ with

$$P_{22}(k) = \int \frac{d^3 q}{4\pi^3} F_2^2(\mathbf{q}, \mathbf{k} - \mathbf{q}) P_{\text{lin}}(q) P_{\text{lin}}(|\mathbf{k} - \mathbf{q}|), \quad (36)$$

$$P_{13}(k) = 3P_{\text{lin}}(k) \int \frac{d^3 q}{4\pi^3} F_3(\mathbf{q}, -\mathbf{q}, \mathbf{k}) P_{\text{lin}}(q), \quad (37)$$

while evaluating the loop integrals $P_{22}$ and $P_{13}$ with the input spectrum $P^{\text{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)$^3$. The IR resummed power spectrum at leading order is given by

$$P^{\text{IR res, LO}}(k, z) = D^2(z) P_{\text{nw}}(k) \times \left[ 1 + e^{-k^2 D^2(z) \Sigma_{\text{BAO}}^2} \delta P_{\text{w}}^{\text{BAO}}(k) + e^{-k^2 D^2(z) \Sigma_{\text{lin}}^2} \delta P_{\text{w}}^{\text{lin}}(k) \right] \quad (38)$$

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

$$P^{\text{IR res, LO}}(k, z) = D^2(z) P_{\text{nw}}(k) \times \left[ 1 + e^{-k^2 D^2(z) \Sigma_{\text{BAO}}^2} \delta P_{\text{w}}^{\text{BAO}}(k) + e^{-k^2 D^2(z) \Sigma_{\text{log}}^2} \cos(k^2 D^2(z) \hat{\Sigma}_{\text{log}}^2) \delta P_{\text{w}}^{\text{log}}(k) - e^{-k^2 D^2(z) \Sigma_{\text{log}}^2} \sin(k^2 D^2(z) \hat{\Sigma}_{\text{log}}^2) \frac{d \delta P_{\text{w}}^{\text{log}}(k)}{\omega_{\text{log}} d \ln k} \right] \quad (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

$$\Sigma_{\text{BAO}}^2(k_S) \equiv \int_0^{k_S} \frac{dq}{6\pi^2} P_{\text{nw}}(q) [1 - j_0(q r_S) + 2 j_2(q r_S)], \quad (40)$$

$$\Sigma_{\text{lin}}^2(\omega_{\text{lin}}, k_S) \equiv \int_0^{k_S} \frac{dq}{6\pi^2} P_{\text{nw}}(q) \left[ 1 - j_0\left(q \frac{\omega_{\text{lin}}}{k_*}\right) + 2j_2\left(q \frac{\omega_{\text{lin}}}{k_*}\right) \right], \quad (41)$$

$$\Sigma_{\text{log}}^2(k, \omega_{\text{log}}, k_S) \equiv \int_0^{k_S} \frac{dq}{4\pi^2} P_{\text{nw}}(q) \times \int_{-1}^{1} d\mu \mu^2 \left\{ 1 - \cos\left[\omega_{\text{log}} \ln\left(1 - \frac{q\mu}{k}\right)\right] \right\}, \quad (42)$$

$$\hat{\Sigma}_{\text{log}}^2(k, \omega_{\text{log}}, k_S) \equiv - \int_0^{k_S} \frac{dq}{4\pi^2} P_{\text{nw}}(q) \int_{-1}^{1} d\mu \mu^2 \sin\left[\omega_{\text{log}} \ln\left(1 - \frac{q\mu}{k}\right)\right], \quad (43)$$

where $j_n$ are spherical Bessel functions, $r_S \approx 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 \ll k$, we assume $k_S = \epsilon k$ with $\epsilon \in [0.3, 0.7]$ (Baldauf et al. 2015; Vasudevan et al. 2019).

In redshift space, the damping factor becomes also dependent on the cosine $\mu$ 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)\Sigma_X \rightarrow D(z)\Sigma_X(\mu)$ with

$$\Sigma_X^2(k, \mu, z) = \Sigma_X^2(z) \{1 - \mu^2 + \mu^2 [1 + f(k, z)]^2\}. \quad (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 $\Lambda$CDM nonlinear matter power spectrum is dressed with the primordial oscillations according to

$$P_{\text{NL}}(k, z) = \frac{P_{\text{NL}}^{\Lambda\text{CDM}}(k, z)}{P_{\text{lin}}^{\Lambda\text{CDM}}(k, z)} P^{\text{IR res, LO}}(k, z). \quad (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)$^4$. 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 $\omega_X \in \{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 $\mathcal{A}_{\text{lin}} = 0.03$ and $\omega_{\text{lin}} = 10$.

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.

<footnote>$^3$ https://github.com/JoeMcEwen/FAST-PT</footnote>

<footnote>$^4$ https://github.com/CullanHowlett/l-picola</footnote>

<footnote>$^5$ https://wwwmpa.mpa-garching.mpg.de/gadget</footnote>
usual expression $P^{1-\text{loop}} = P_{22} + 2P_{13}$ with

$$P_{22}(k) = \int \frac{d^3 q}{4\pi^3} F_2^2(\mathbf{q}, \mathbf{k} - \mathbf{q}) P_{\text{lin}}(q) P_{\text{lin}}(|\mathbf{k} - \mathbf{q}|), \quad (36)$$

$$P_{13}(k) = 3P_{\text{lin}}(k) \int \frac{d^3 q}{4\pi^3} F_3(\mathbf{q}, -\mathbf{q}, \mathbf{k}) P_{\text{lin}}(q), \quad (37)$$

while evaluating the loop integrals $P_{22}$ and $P_{13}$ with the input spectrum $P^{\text{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)$^3$. The IR resummed power spectrum at leading order is given by

$$P^{\text{IR res, LO}}(k, z) = D^2(z) P_{\text{nw}}(k) \times \left[ 1 + e^{-k^2 D^2(z) \Sigma_{\text{BAO}}^2} \delta P_{\text{w}}^{\text{BAO}}(k) + e^{-k^2 D^2(z) \Sigma_{-\text{lin}}^2} \delta P_{\text{w}}^{\text{lin}}(k) \right] \quad (38)$$

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

$$P^{\text{IR res, LO}}(k, z) = D^2(z) P_{\text{nw}}(k) \times \left[ 1 + e^{-k^2 D^2(z) \Sigma_{\text{BAO}}^2} \delta P_{\text{w}}^{\text{BAO}}(k) + e^{-k^2 D^2(z) \Sigma_{\text{log}}^2} \cos(k^2 D^2(z) \hat{\Sigma}_{\text{log}}^2) \delta P_{\text{w}}^{\text{log}}(k) - e^{-k^2 D^2(z) \Sigma_{\text{log}}^2} \sin(k^2 D^2(z) \hat{\Sigma}_{\text{log}}^2) \frac{d \delta P_{\text{w}}^{\text{log}}(k)}{\omega_{\text{log}} d \ln k} \right] \quad (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

$$\Sigma_{\text{BAO}}^2(k_S) = \int_0^{k_S} \frac{dq}{6\pi^2} P_{\text{nw}}(q) [1 - j_0(qr_S) + 2j_2(qr_S)], \quad (40)$$

$$\Sigma_{\text{lin}}^2(\omega_{\text{lin}}, k_S) = \int_0^{k_S} \frac{dq}{6\pi^2} P_{\text{nw}}(q) \left[ 1 - j_0\left(q \frac{\omega_{\text{lin}}}{k_*}\right) + 2j_2\left(q \frac{\omega_{\text{lin}}}{k_*}\right) \right], \quad (41)$$

$$\Sigma_{\text{log}}^2(k, \omega_{\text{log}}, k_S) = \int_0^{k_S} \frac{dq}{4\pi^2} P_{\text{nw}}(q) \times \int_{-1}^{1} d\mu\mu^2 \left\{ 1 - \cos\left[\omega_{\text{log}} \ln\left(1 - \frac{q\mu}{k}\right)\right] \right\}, \quad (42)$$

$$\hat{\Sigma}_{\text{log}}^2(k, \omega_{\text{log}}, k_S) = - \int_0^{k_S} \frac{dq}{4\pi^2} P_{\text{nw}}(q) \int_{-1}^{1} d\mu\mu^2 \sin\left[\omega_{\text{log}} \ln\left(1 - \frac{q\mu}{k}\right)\right], \quad (43)$$

where $j_n$ are spherical Bessel functions, $r_S \approx 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 \ll k$, we assume $k_S = \epsilon k$ with $\epsilon \in [0.3, 0.7]$ (Baldauf et al. 2015; Vasudevan et al. 2019).

In redshift space, the damping factor becomes also dependent on the cosine $\mu$ 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)\Sigma_X \rightarrow D(z)\Sigma_X(\mu)$ with

$$\Sigma_X^2(k, \mu, z) = \Sigma_X^2(z) \{1 - \mu^2 + \mu^2 [1 + f(k, z)]^2\}. \quad (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

$$P_{\text{NL}}(k, z) = \frac{P_{\text{NL}}^{\Lambda\text{CDM}}(k, z)}{P_{\text{lin}}^{\Lambda\text{CDM}}(k, z)} P^{\text{IR res, LO}}(k, z). \quad (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)$^4$. 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 $\omega_X \in \{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} \text{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 $\mathcal{A}_{\text{lin}} = 0.03$ and $\omega_{\text{lin}} = 10$.

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.

<footnote>$^3$ https://github.com/JoeMcEwen/FAST-PT</footnote>

<footnote>$^4$ https://github.com/CullanHowlett/l-picola</footnote>

<footnote>$^5$ https://wwwmpa.mpa-garching.mpg.de/gadget</footnote>
<figure><img src="image_3.png" /><figcaption>**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 $\mathcal{A}_{\text{lin}} = 0.03$ and $\omega_{\text{lin}} = 10$ (Ballardini et al. 2020).</figcaption></figure>

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} \omega_{\text{lin}} \in (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} \omega_{\text{lin}} > 0.6$ and between 1% and 2% for $\log_{10} \omega_{\text{lin}} \le 0.6$. At NLO, we find an agreement better than 1% for all the frequencies and differences <0.5% for $\log_{10} \omega_{\text{lin}} > 0.4$. For the template with undamped logarithmic oscillations, we find a similar agreement for the frequencies $\log_{10} \omega_{\text{log}} \ge 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

<figure><img src="image_4.png" /><figcaption>**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 $\omega_X = 10$ and the IR separation scale $k_S = \epsilon k$ is varied (with $\epsilon \in [0.3, 0.7]$). We also show the results of the linear theory (grey) and the ones obtained from N-body simulations (green).</figcaption></figure>

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 \le k/(h \, \text{Mpc}^{-1}) \le 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 \le 0.2 h \, \text{Mpc}^{-1}$ or $s \ge 15–20 h^{-1} \text{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
<figure><img src="image_3.png" /><figcaption>**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 $\mathcal{A}_{\text{lin}} = 0.03$ and $\omega_{\text{lin}} = 10$ (Ballardini et al. 2020).</figcaption></figure>

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} \omega_{\text{lin}} \in (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} \omega_{\text{lin}} > 0.6$ and between 1% and 2% for $\log_{10} \omega_{\text{lin}} \le 0.6$. At NLO, we find an agreement better than 1% for all the frequencies and differences <0.5% for $\log_{10} \omega_{\text{lin}} > 0.4$. For the template with undamped logarithmic oscillations, we find a similar agreement for the frequencies $\log_{10} \omega_{\text{log}} \ge 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

<figure><img src="image_4.png" /><figcaption>**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 $\omega_X = 10$ and the IR separation scale $k_S = \epsilon k$ is varied (with $\epsilon \in [0.3, 0.7]$). We also show the results of the linear theory (grey) and the ones obtained from N-body simulations (green).</figcaption></figure>

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 \le k/(h \, \text{Mpc}^{-1}) \le 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 \le 0.2 h \, \text{Mpc}^{-1}$ or $s \ge 15–20 h^{-1} \text{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

$$ \rho(\mathbf{x}) \, d^3\mathbf{x} = \rho(\mathbf{q}) \, d^3\mathbf{q} \approx \bar{\rho} \, d^3\mathbf{q}, \qquad (46) $$

where $\rho(\mathbf{q})$ and $\rho(\mathbf{x})$ are, respectively, the initial density field in some infinitesimal volume element $d\mathbf{q}$, and the density field in the same volume element at time $t$, which now has been moved, and possibly deformed, and described by $d\mathbf{x}$. As the density field in the early Universe is homogeneous to a very good approximation, one can take $\rho(\mathbf{q})$ to be a constant, i.e. $\rho(\mathbf{q}) \approx \bar{\rho}$.

We define the displacement field $\mathbf{\Psi}$ between the final and initial positions of a particle as $\mathbf{\Psi}(\mathbf{x}) = \mathbf{x} - \mathbf{q}$, which can be rewritten as

$$ \nabla_{\mathbf{x}} \Theta(\mathbf{x}) \equiv \mathbf{q} = \mathbf{x} - \mathbf{\Psi}(\mathbf{x}), \qquad (47) $$

with $\Theta(\mathbf{x})$ being the displacement potential. Here we have assumed the displacement field to be curl-free, i.e. $\nabla \times \mathbf{\Psi} = \mathbf{0}$, which is valid only on large scales$^6$. Combining Eqs. (47) and (46) gives

$$ \det[\nabla^i \nabla_j \Theta(\mathbf{x})] = \frac{\rho(\mathbf{x})}{\bar{\rho}} \equiv 1 + \delta(\mathbf{x}), \qquad (48) $$

where $\delta(\mathbf{x}) = \rho(\mathbf{x})/\bar{\rho} - 1$ is the density contrast at time $t$, ‘det’ is the determinant of a matrix, here the Hessian of $\Theta(\mathbf{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 $\Theta$. This was later extended by Birkin et al. (2019) to cases where $\delta(\mathbf{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 $\Theta(\mathbf{x})$ is solved, one can obtain $\mathbf{\Psi}(\mathbf{x})$ and the reconstructed density field is given by

$$ \delta_r = -\nabla_q \cdot \mathbf{\Psi}(\mathbf{q}), \qquad (49) $$

where we express $\mathbf{\Psi}$ in terms of the Lagrangian coordinate $\mathbf{q}$, since the divergence $\nabla_q$ therein is with respect to $\mathbf{q}$. This calculation means that we need to have $\mathbf{\Psi}(\mathbf{q})$ on a regular $\mathbf{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 $\mathbf{\Psi}(\mathbf{x})$ to the target $\mathbf{q}$-grid.

The quantity $\delta_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 $\Lambda$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_w^{\text{undamped}}$ by a Gaussian function

$$ O_w^{\text{damped}}(k, z) = O_w^{\text{undamped}}(k, z) \exp \left[ -\frac{k^2 \zeta(z)^2}{2} \right], \qquad (50) $$

<table><caption>**Table 1.** Best-fit values of $\zeta(z)$ for unreconstructed (reconstructed) wiggle spectra.</caption><thead><tr><th>z</th><th>Unrec</th><th>Rec</th><th>Efficiency</th></tr></thead><tbody><tr><td>0.0</td><td>8.16</td><td>2.04</td><td>4.00</td></tr><tr><td>0.5</td><td>6.49</td><td>1.69</td><td>3.84</td></tr><tr><td>1.0</td><td>5.23</td><td>1.41</td><td>3.71</td></tr><tr><td>1.5</td><td>4.32</td><td>1.21</td><td>3.57</td></tr></tbody></table>

**Notes.** Columns respectively denote (1) redshift, (2) unreconstructed wiggles, (3) reconstructed wiggles, and (4) reconstruction efficiency (defined as $\zeta_{\text{unrec}}/\zeta_{\text{rec}}$).

where $\zeta(z)$ is the parameter quantifying the level of damping. In the linearly evolved density field with no damping, one has $\zeta = 0$. For the evolved nonlinear matter or galaxy field, $\zeta > 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 $\zeta$.

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 \approx 1 h \text{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 \text{Mpc}^{-1}$.

To estimate the best-fit values of $\zeta(z)$, the Gaussian function is applied to fit the envelopes of the unreconstructed and reconstructed wiggled spectra, $O_w^{\text{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 $\zeta(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_w^{\text{damped}}$ and $O_w^{\text{undamped}}$ of the different models, and finding them to agree very well. As a result, to get $\zeta(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 \text{Mpc}^{-1}$, we select only the 10 high-frequency feature models to fit $\zeta(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

<footnote>$^6$ This is similar to the assumption of no particle trajectory crossing, which must break down on small scales. Together, these assumptions mean that the reconstruction method described here, similar to other reconstruction methods, should only be expected to work for relatively large scales, and its performance progressively degrades if one goes to smaller scales.</footnote>
mapping is unique and can be obtained by solving the mass conservation equation

$$ \rho(\mathbf{x}) \, d^3\mathbf{x} = \rho(\mathbf{q}) \, d^3\mathbf{q} \approx \bar{\rho} \, d^3\mathbf{q}, \qquad (46) $$

where $\rho(\mathbf{q})$ and $\rho(\mathbf{x})$ are, respectively, the initial density field in some infinitesimal volume element $d\mathbf{q}$, and the density field in the same volume element at time $t$, which now has been moved, and possibly deformed, and described by $d\mathbf{x}$. As the density field in the early Universe is homogeneous to a very good approximation, one can take $\rho(\mathbf{q})$ to be a constant, i.e. $\rho(\mathbf{q}) \approx \bar{\rho}$.

We define the displacement field $\mathbf{\Psi}$ between the final and initial positions of a particle as $\mathbf{\Psi}(\mathbf{x}) = \mathbf{x} - \mathbf{q}$, which can be rewritten as

$$ \nabla_{\mathbf{x}} \Theta(\mathbf{x}) \equiv \mathbf{q} = \mathbf{x} - \mathbf{\Psi}(\mathbf{x}), \qquad (47) $$

with $\Theta(\mathbf{x})$ being the displacement potential. Here we have assumed the displacement field to be curl-free, i.e. $\nabla \times \mathbf{\Psi} = \mathbf{0}$, which is valid only on large scales<fnref n="6" />. Combining Eqs. (47) and (46) gives

$$ \det[\nabla^i \nabla_j \Theta(\mathbf{x})] = \frac{\rho(\mathbf{x})}{\bar{\rho}} \equiv 1 + \delta(\mathbf{x}), \qquad (48) $$

where $\delta(\mathbf{x}) = \rho(\mathbf{x})/\bar{\rho} - 1$ is the density contrast at time $t$, ‘det’ is the determinant of a matrix, here the Hessian of $\Theta(\mathbf{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 $\Theta$. This was later extended by Birkin et al. (2019) to cases where $\delta(\mathbf{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 $\Theta(\mathbf{x})$ is solved, one can obtain $\mathbf{\Psi}(\mathbf{x})$ and the reconstructed density field is given by

$$ \delta_r = -\nabla_q \cdot \mathbf{\Psi}(\mathbf{q}), \qquad (49) $$

where we express $\mathbf{\Psi}$ in terms of the Lagrangian coordinate $\mathbf{q}$, since the divergence $\nabla_q$ therein is with respect to $\mathbf{q}$. This calculation means that we need to have $\mathbf{\Psi}(\mathbf{q})$ on a regular $\mathbf{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 $\mathbf{\Psi}(\mathbf{x})$ to the target $\mathbf{q}$-grid.

The quantity $\delta_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 $\Lambda$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_w^{\text{undamped}}$ by a Gaussian function

$$ O_w^{\text{damped}}(k, z) = O_w^{\text{undamped}}(k, z) \exp \left[ -\frac{k^2 \zeta(z)^2}{2} \right], \qquad (50) $$

<table><caption>**Table 1.** Best-fit values of $\zeta(z)$ for unreconstructed (reconstructed) wiggle spectra.</caption><thead><tr><th>z</th><th>Unrec</th><th>Rec</th><th>Efficiency</th></tr></thead><tbody><tr><td>0.0</td><td>8.16</td><td>2.04</td><td>4.00</td></tr><tr><td>0.5</td><td>6.49</td><td>1.69</td><td>3.84</td></tr><tr><td>1.0</td><td>5.23</td><td>1.41</td><td>3.71</td></tr><tr><td>1.5</td><td>4.32</td><td>1.21</td><td>3.57</td></tr></tbody></table>

**Notes.** Columns respectively denote (1) redshift, (2) unreconstructed wiggles, (3) reconstructed wiggles, and (4) reconstruction efficiency (defined as $\zeta_{\text{unrec}}/\zeta_{\text{rec}}$).

where $\zeta(z)$ is the parameter quantifying the level of damping. In the linearly evolved density field with no damping, one has $\zeta = 0$. For the evolved nonlinear matter or galaxy field, $\zeta > 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 $\zeta$.

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 \approx 1 h \text{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 \text{Mpc}^{-1}$.

To estimate the best-fit values of $\zeta(z)$, the Gaussian function is applied to fit the envelopes of the unreconstructed and reconstructed wiggled spectra, $O_w^{\text{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 $\zeta(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_w^{\text{damped}}$ and $O_w^{\text{undamped}}$ of the different models, and finding them to agree very well. As a result, to get $\zeta(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 \text{Mpc}^{-1}$, we select only the 10 high-frequency feature models to fit $\zeta(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

<footnote n="6">⁶ This is similar to the assumption of no particle trajectory crossing, which must break down on small scales. Together, these assumptions mean that the reconstruction method described here, similar to other reconstruction methods, should only be expected to work for relatively large scales, and its performance progressively degrades if one goes to smaller scales.</footnote>
<table><caption>**Table 2.** Fisher-forecast 68.3% CL marginalised uncertainties on cosmological and primordial feature parameters, relative to their corresponding fiducial values.</caption><thead><tr><th></th><th colspan="8">LIN ω<sub>lin</sub> = 10</th></tr><tr><th></th><th>Ω<sub>m,0</sub></th><th>Ω<sub>b,0</sub></th><th>h</th><th>n<sub>s</sub></th><th>σ<sub>8</sub></th><th>A<sub>lin</sub></th><th>ω<sub>lin</sub></th><th>φ<sub>lin</sub></th></tr></thead><tbody><tr><td colspan="9"><strong>Pessimistic setting</strong></td></tr><tr><td>GC<sub>sp</sub> (k<sub>max</sub> = 0.25 h Mpc<sup>-1</sup>)</td><td>1.2%</td><td>1.7%</td><td>1.4%</td><td>1.5%</td><td>0.83%</td><td>23%</td><td>4.1%</td><td>0.11</td></tr><tr><td>WL+GC<sub>ph</sub>+XC</td><td>0.83%</td><td>5.7%</td><td>4.0%</td><td>1.6%</td><td>0.39%</td><td>32%</td><td>5.7%</td><td>0.12</td></tr><tr><td>GC<sub>sp</sub>+WL+GC<sub>ph</sub>+XC(z &lt; 0.9)</td><td>0.61%</td><td>1.3%</td><td>0.59%</td><td>0.70%</td><td>0.28%</td><td>21%</td><td>2.9%</td><td>0.083</td></tr><tr><td colspan="9"><strong>Optimistic setting</strong></td></tr><tr><td>GC<sub>sp</sub> (k<sub>max</sub> = 0.3 h Mpc<sup>-1</sup>)</td><td>1.1%</td><td>1.6%</td><td>1.1%</td><td>1.3%</td><td>0.76%</td><td>23%</td><td>3.6%</td><td>0.10</td></tr><tr><td>WL+GC<sub>ph</sub>+XC</td><td>0.27%</td><td>4.4%</td><td>2.5%</td><td>0.66%</td><td>0.13%</td><td>31%</td><td>2.8%</td><td>0.072</td></tr><tr><td>GC<sub>sp</sub>+WL+GC<sub>ph</sub>+XC</td><td>0.23%</td><td>1.1%</td><td>0.45%</td><td>0.28%</td><td>0.11%</td><td>18%</td><td>1.2%</td><td>0.047</td></tr></tbody><thead><tr><th></th><th colspan="8">LOG ω<sub>log</sub> = 10</th></tr><tr><th></th><th>Ω<sub>m,0</sub></th><th>Ω<sub>b,0</sub></th><th>h</th><th>n<sub>s</sub></th><th>σ<sub>8</sub></th><th>A<sub>log</sub></th><th>ω<sub>log</sub></th><th>φ<sub>log</sub></th></tr></thead><tbody><tr><td colspan="9"><strong>Pessimistic setting</strong></td></tr><tr><td>GC<sub>sp</sub> (k<sub>max</sub> = 0.25 h Mpc<sup>-1</sup>)</td><td>1.1%</td><td>1.8%</td><td>1.4%</td><td>1.5%</td><td>0.82%</td><td>25%</td><td>4.7%</td><td>0.053</td></tr><tr><td>WL+GC<sub>ph</sub>+XC</td><td>0.83%</td><td>5.7%</td><td>4.0%</td><td>1.6%</td><td>0.38%</td><td>32%</td><td>5.7%</td><td>0.12</td></tr><tr><td>GC<sub>sp</sub>+WL+GC<sub>ph</sub>+XC(z &lt; 0.9)</td><td>0.61%</td><td>1.4%</td><td>0.59%</td><td>0.71%</td><td>0.28%</td><td>22%</td><td>2.4%</td><td>0.044</td></tr><tr><td colspan="9"><strong>Optimistic setting</strong></td></tr><tr><td>GC<sub>sp</sub> (k<sub>max</sub> = 0.3 h Mpc<sup>-1</sup>)</td><td>1.1%</td><td>1.8%</td><td>1.1%</td><td>1.3%</td><td>0.76%</td><td>25%</td><td>4.7%</td><td>0.052</td></tr><tr><td>WL+GC<sub>ph</sub>+XC</td><td>0.27%</td><td>4.4%</td><td>2.5%</td><td>0.66%</td><td>0.13%</td><td>31%</td><td>2.8%</td><td>0.072</td></tr><tr><td>GC<sub>sp</sub>+WL+GC<sub>ph</sub>+XC</td><td>0.23%</td><td>1.1%</td><td>0.46%</td><td>0.28%</td><td>0.11%</td><td>18%</td><td>1.1%</td><td>0.035</td></tr></tbody></table>

**Notes.** We show results for LIN and LOG models in the pessimistic and optimistic settings, using Euclid observations GC<sub>sp</sub>, WL+GC<sub>ph</sub>+XC, and their combination.

flat ΛCDM,

$$ \Theta_{\text{final}} = \{\Omega_{\text{m,0}}, \Omega_{\text{b,0}}, h, n_s, \sigma_8, \mathcal{A}_X, \omega_X, \phi_X\}, \quad (51) $$

after marginalisation over spectroscopic and photometric nuisance parameters. We compute the spectroscopic galaxy clustering (GC<sub>sp</sub>), the photometric galaxy clustering (GC<sub>ph</sub>), 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<sub>sp</sub> we have k<sub>max</sub> = 0.25 h Mpc<sup>-1</sup>, for GC<sub>ph</sub> and XC we have l<sub>max</sub> = 750, and for WL we have l<sub>max</sub> = 1500. In addition, we impose a cut in redshift of z < 0.9 to GC<sub>ph</sub> to limit any possible cross-correlation between spectroscopic and photometric data. For the optimistic setting, we extend GC<sub>sp</sub> to k<sub>max</sub> = 0.30 h Mpc<sup>-1</sup>, GC<sub>ph</sub> and XC to l<sub>max</sub> = 3000, and WL to l<sub>max</sub> = 5000, without imposing any redshift cut to GC<sub>ph</sub>.

We start by presenting results for a fiducial scenario with A<sub>X</sub> = 0.01, ω<sub>X</sub> = 10, and φ<sub>X</sub> = 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<sub>sp</sub> measurements resulting in A<sub>lin</sub> = 0.0100 ± 0.0023 (±0.0023) and A<sub>log</sub> = 0.0100 ± 0.0025 (±0.0025) at a 68.3% CL for the pessimistic (optimistic) setting. Combining GC<sub>sp</sub> with the combination of photometric information, i.e. WL+GC<sub>ph</sub>+XC, we find A<sub>lin</sub> = 0.0100 ± 0.0021 (±0.0018) and A<sub>log</sub> = 0.0100 ± 0.0022 (±0.0018) at a 68.3% CL for the pessimistic (optimistic) setting.

For GC<sub>sp</sub>, bounds on the primordial amplitude parameter A<sub>X</sub> strongly depend on the primordial frequency value ω<sub>X</sub>; 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<sub>10</sub> ω<sub>lin</sub> ~ 0.87, are degraded. For the combination of photometric probes, bounds on the primordial amplitude parameter A<sub>X</sub> 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<sub>sp</sub>. 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<sub>10</sub> ω<sub>X</sub> 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.

The nonlinear reconstruction described in Sect. 4.4 applied to GC<sub>sp</sub> reduces significantly the uncertainties on the feature amplitude, going from σ(A<sub>lin</sub>) = 0.0023 (0.0023) at a 68.3% CL for the pessimistic (optimistic) setting to σ(A<sub>lin</sub>) = 0.0014 (0.0013), and going from σ(A<sub>log</sub>) = 0.0025 (0.0025) to σ(A<sub>log</sub>) = 0.0016 (0.0012). Combining with the photometric probes, we find σ(A<sub>lin</sub>) = 0.0014 (0.0012) at a 68.3% CL and σ(A<sub>log</sub>) = 0.0015 (0.0012) for the pessimistic (optimistic) setting, compared with σ(A<sub>lin</sub>) = 0.0021 (0.0018) at a 68.3% CL and σ(A<sub>log</sub>) = 0.0022 (0.0018) for the same settings without reconstruction. We also note that in the analysis here k<sub>max</sub> is set to be 0.25–0.30 h Mpc<sup>-1</sup>, while in Li et al. (2022) k<sub>max</sub> = 0.50 h Mpc<sup>-1</sup> is used.
<table><caption>**Table 2.** Fisher-forecast 68.3% CL marginalised uncertainties on cosmological and primordial feature parameters, relative to their corresponding fiducial values.</caption><thead><tr><th rowspan="2"></th><th colspan="7">LIN ω<sub>lin</sub> = 10</th></tr><tr><th>Ω<sub>m,0</sub></th><th>Ω<sub>b,0</sub></th><th>h</th><th>n<sub>s</sub></th><th>σ<sub>8</sub></th><th>A<sub>lin</sub></th><th>ω<sub>lin</sub></th><th>φ<sub>lin</sub></th></tr></thead><tbody><tr><td colspan="9"><strong>Pessimistic setting</strong></td></tr><tr><td>GC<sub>sp</sub> (k<sub>max</sub> = 0.25 h Mpc<sup>-1</sup>)</td><td>1.2%</td><td>1.7%</td><td>1.4%</td><td>1.5%</td><td>0.83%</td><td>23%</td><td>4.1%</td><td>0.11</td></tr><tr><td>WL+GC<sub>ph</sub>+XC</td><td>0.83%</td><td>5.7%</td><td>4.0%</td><td>1.6%</td><td>0.39%</td><td>32%</td><td>5.7%</td><td>0.12</td></tr><tr><td>GC<sub>sp</sub>+WL+GC<sub>ph</sub>+XC(z &lt; 0.9)</td><td>0.61%</td><td>1.3%</td><td>0.59%</td><td>0.70%</td><td>0.28%</td><td>21%</td><td>2.9%</td><td>0.083</td></tr><tr><td colspan="9"><strong>Optimistic setting</strong></td></tr><tr><td>GC<sub>sp</sub> (k<sub>max</sub> = 0.3 h Mpc<sup>-1</sup>)</td><td>1.1%</td><td>1.6%</td><td>1.1%</td><td>1.3%</td><td>0.76%</td><td>23%</td><td>3.6%</td><td>0.10</td></tr><tr><td>WL+GC<sub>ph</sub>+XC</td><td>0.27%</td><td>4.4%</td><td>2.5%</td><td>0.66%</td><td>0.13%</td><td>31%</td><td>2.8%</td><td>0.072</td></tr><tr><td>GC<sub>sp</sub>+WL+GC<sub>ph</sub>+XC</td><td>0.23%</td><td>1.1%</td><td>0.45%</td><td>0.28%</td><td>0.11%</td><td>18%</td><td>1.2%</td><td>0.047</td></tr></tbody><thead><tr><th rowspan="2"></th><th colspan="7">LOG ω<sub>log</sub> = 10</th></tr><tr><th>Ω<sub>m,0</sub></th><th>Ω<sub>b,0</sub></th><th>h</th><th>n<sub>s</sub></th><th>σ<sub>8</sub></th><th>A<sub>log</sub></th><th>ω<sub>log</sub></th><th>φ<sub>log</sub></th></tr></thead><tbody><tr><td colspan="9"><strong>Pessimistic setting</strong></td></tr><tr><td>GC<sub>sp</sub> (k<sub>max</sub> = 0.25 h Mpc<sup>-1</sup>)</td><td>1.1%</td><td>1.8%</td><td>1.4%</td><td>1.5%</td><td>0.82%</td><td>25%</td><td>4.7%</td><td>0.053</td></tr><tr><td>WL+GC<sub>ph</sub>+XC</td><td>0.83%</td><td>5.7%</td><td>4.0%</td><td>1.6%</td><td>0.38%</td><td>32%</td><td>5.7%</td><td>0.12</td></tr><tr><td>GC<sub>sp</sub>+WL+GC<sub>ph</sub>+XC(z &lt; 0.9)</td><td>0.61%</td><td>1.4%</td><td>0.59%</td><td>0.71%</td><td>0.28%</td><td>22%</td><td>2.4%</td><td>0.044</td></tr><tr><td colspan="9"><strong>Optimistic setting</strong></td></tr><tr><td>GC<sub>sp</sub> (k<sub>max</sub> = 0.3 h Mpc<sup>-1</sup>)</td><td>1.1%</td><td>1.8%</td><td>1.1%</td><td>1.3%</td><td>0.76%</td><td>25%</td><td>4.7%</td><td>0.052</td></tr><tr><td>WL+GC<sub>ph</sub>+XC</td><td>0.27%</td><td>4.4%</td><td>2.5%</td><td>0.66%</td><td>0.13%</td><td>31%</td><td>2.8%</td><td>0.072</td></tr><tr><td>GC<sub>sp</sub>+WL+GC<sub>ph</sub>+XC</td><td>0.23%</td><td>1.1%</td><td>0.46%</td><td>0.28%</td><td>0.11%</td><td>18%</td><td>1.1%</td><td>0.035</td></tr></tbody></table>

**Notes.** We show results for LIN and LOG models in the pessimistic and optimistic settings, using Euclid observations GC<sub>sp</sub>, WL+GC<sub>ph</sub>+XC, and their combination.

flat ΛCDM,

$$ \Theta_{\text{final}} = \{\Omega_{\text{m,0}}, \Omega_{\text{b,0}}, h, n_s, \sigma_8, \mathcal{A}_X, \omega_X, \phi_X\}, \quad (51) $$

after marginalisation over spectroscopic and photometric nuisance parameters. We compute the spectroscopic galaxy clustering (GC<sub>sp</sub>), the photometric galaxy clustering (GC<sub>ph</sub>), 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<sub>sp</sub> we have k<sub>max</sub> = 0.25 h Mpc<sup>-1</sup>, for GC<sub>ph</sub> and XC we have l<sub>max</sub> = 750, and for WL we have l<sub>max</sub> = 1500. In addition, we impose a cut in redshift of z < 0.9 to GC<sub>ph</sub> to limit any possible cross-correlation between spectroscopic and photometric data. For the optimistic setting, we extend GC<sub>sp</sub> to k<sub>max</sub> = 0.30 h Mpc<sup>-1</sup>, GC<sub>ph</sub> and XC to l<sub>max</sub> = 3000, and WL to l<sub>max</sub> = 5000, without imposing any redshift cut to GC<sub>ph</sub>.

We start by presenting results for a fiducial scenario with A<sub>X</sub> = 0.01, ω<sub>X</sub> = 10, and φ<sub>X</sub> = 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<sub>sp</sub> measurements resulting in A<sub>lin</sub> = 0.0100 ± 0.0023 (±0.0023) and A<sub>log</sub> = 0.0100 ± 0.0025 (±0.0025) at a 68.3% CL for the pessimistic (optimistic) setting. Combining GC<sub>sp</sub> with the combination of photometric information, i.e. WL+GC<sub>ph</sub>+XC, we find A<sub>lin</sub> = 0.0100 ± 0.0021 (±0.0018) and A<sub>log</sub> = 0.0100 ± 0.0022 (±0.0018) at a 68.3% CL for the pessimistic (optimistic) setting.

For GC<sub>sp</sub>, bounds on the primordial amplitude parameter A<sub>X</sub> strongly depend on the primordial frequency value ω<sub>X</sub>; 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<sub>10</sub> ω<sub>lin</sub> ~ 0.87, are degraded. For the combination of photometric probes, bounds on the primordial amplitude parameter A<sub>X</sub> 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<sub>sp</sub>. 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<sub>10</sub> ω<sub>X</sub> 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.

The nonlinear reconstruction described in Sect. 4.4 applied to GC<sub>sp</sub> reduces significantly the uncertainties on the feature amplitude, going from σ(A<sub>lin</sub>) = 0.0023 (0.0023) at a 68.3% CL for the pessimistic (optimistic) setting to σ(A<sub>lin</sub>) = 0.0014 (0.0013), and going from σ(A<sub>log</sub>) = 0.0025 (0.0025) to σ(A<sub>log</sub>) = 0.0016 (0.0012). Combining with the photometric probes, we find σ(A<sub>lin</sub>) = 0.0014 (0.0012) at a 68.3% CL and σ(A<sub>log</sub>) = 0.0015 (0.0012) for the pessimistic (optimistic) setting, compared with σ(A<sub>lin</sub>) = 0.0021 (0.0018) at a 68.3% CL and σ(A<sub>log</sub>) = 0.0022 (0.0018) for the same settings without reconstruction. We also note that in the analysis here k<sub>max</sub> is set to be 0.25–0.30 h Mpc<sup>-1</sup>, while in Li et al. (2022) k<sub>max</sub> = 0.50 h Mpc<sup>-1</sup> is used.
<figure><img src="image_5.png" /><figcaption>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 $\omega_{\text{lin}} = 10$ (top panels) and the LOG model with $\omega_{\text{log}} = 10$ (bottom panels). Left (right) panels correspond to the pessimistic (optimistic) setting for GC<sub>sp</sub> (green), WL+GC<sub>ph</sub>+XC (orange), and their combination (blue).</figcaption></figure>

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

$$B_{\Phi}^{\text{lin}}(k_1, k_2, k_3) = f_{\text{NL}}^{\text{lin}} \frac{6A^2}{k_1^2 k_2^2 k_3^2} \sin \left[ \omega_{\text{lin}} \frac{K}{k_*} + 2\pi \phi_{\text{lin}}^B \right], \quad (52)$$

where $K = k_1 + k_2 + k_3$ and $A = 9/25(2\pi^2)k_*^{1-n_s}A_s$ is the normalisation parameter of $P_{\Phi} = 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)

$$B_{\Phi}^{\text{log}}(k_1, k_2, k_3) = f_{\text{NL}}^{\text{log}} \frac{6A^2}{k_1^2 k_2^2 k_3^2} \sin \left[ \omega_{\text{log}} \ln \left( \frac{K}{k_*} \right) + 2\pi \phi_{\text{log}}^B \right], \quad (53)$$
<figure><img src="image_5.png" /><figcaption>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 $\omega_{\text{lin}} = 10$ (top panels) and the LOG model with $\omega_{\text{log}} = 10$ (bottom panels). Left (right) panels correspond to the pessimistic (optimistic) setting for GC<sub>sp</sub> (green), WL+GC<sub>ph</sub>+XC (orange), and their combination (blue).</figcaption></figure>

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

$$B_{\Phi}^{\text{lin}}(k_1, k_2, k_3) = f_{\text{NL}}^{\text{lin}} \frac{6A^2}{k_1^2 k_2^2 k_3^2} \sin\left[\omega_{\text{lin}} \frac{K}{k_*} + 2\pi\phi_{\text{lin}}^B\right], \quad (52)$$

where $K = k_1 + k_2 + k_3$ and $A = 9/25(2\pi^2)k_*^{1-n_s}A_s$ is the normalisation parameter of $P_{\Phi} = 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)

$$B_{\Phi}^{\text{log}}(k_1, k_2, k_3) = f_{\text{NL}}^{\text{log}} \frac{6A^2}{k_1^2 k_2^2 k_3^2} \sin\left[\omega_{\text{log}} \ln\left(\frac{K}{k_*}\right) + 2\pi\phi_{\text{log}}^B\right], \quad (53)$$
<figure><img src="image_6.png" /><figcaption>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.</figcaption></figure>

where $\phi_{\log}^B$ and $f_{\text{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

$$P_g(k, z) = D^P(k, z)Z_1^2(k, z)P^{\text{lin}}(k, z) + P_\varepsilon(z), \quad (54)$$

$$B_g(k_1, k_2, k_3, z) = D^B(k_1, k_2, k_3, z) \times \left[ Z_1(k_1, z)Z_1(k_2, z)Z_1(k_3, z)B_I(k_1, k_2, k_3, z) + 2Z_1(k_1, z)Z_1(k_2, z)Z_2(k_1, k_2, z)P_{\text{lin}}(k_1, z)P_{\text{lin}}(k_2, z) + 2 \text{perm} \right] + 2P_{\varepsilon\varepsilon\delta}[Z_1(k_1, z)P_{\text{lin}}(k_1, z) + 2 \text{perm}] + B_\varepsilon, \quad (55)$$

where the linearly propagated primordial oscillatory bispectrum is

$$B_I(k_1, k_2, k_3, z) = M(k_1, z)M(k_2, z)M(k_3, z)B_\Phi^X(k_1, k_2, k_3, z) \quad (56)$$

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

given by

$$Z_1(k_i, z) = b_1(z) + f\mu_i^2, \quad (57)$$

$$Z_2(k_i, k_j, z) = b_1(z)F_2(k_i, k_j) + f(z)\mu_i^2 G_2(k_i, k_j) + \frac{b_2(z)}{2} + \frac{b_{s2}(z)}{2}S_2(k_i, k_j) + \frac{f(z)\mu_{ij}k_{ij}}{2} \left[ \frac{\mu_i}{k_i}Z_1(k_j, z) + \frac{\mu_j}{k_j}Z_1(k_i, z) \right], \quad (58)$$

where $\mu_i = k_i \cdot \hat{z}/k_i$ is the cosine of the angle between the wave vector $k_i$ and the line of sight, $\mu = \mu_1 + \dots + \mu_n$, and $k = |k_1 + \dots + k_n|$. Specifically $\mu_{12} = (\mu_1 k_1 + \mu_2 k_2)/k_{12}$ and $k_{12}^2 = (k_1 + k_2)^2$. 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), while $S_2(k_1, k_2) = (k_1 \cdot 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_{s2}(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. $B_g(k_1, k_2, k_3, z) = B_g(k_1, k_2, k_3, \mu_1, \phi, 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 upcom- ing 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 degen- eracy between the scale-dependent bias coefficient and the pri- mordial bispectrum amplitude ($f_{\text{NL}}^X$), 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\mu(\sigma_p^2(z) + \sigma_f^2(z))]$ and $D^B(k_1, k_2, k_3, z) = \exp[-(k_1^2\mu_1^2 + k_2^2\mu_2^2 + k_3^2\mu_3^2)(\sigma_p^2(z) + \sigma_f^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_\varepsilon = 1/\bar{n}_g$, $P_{\varepsilon\varepsilon\delta} = b_1/(2\bar{n}_g)$, and $B_\varepsilon = 1/\bar{n}_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)

$$B_g^{\text{obs}}(k_{\text{ref},1}, k_{\text{ref},2}, k_{\text{ref},3}, \mu_{\text{ref},1}, \phi; z) = \frac{1}{q_{\parallel}^2 q_{\perp}^4} B_g(k_{\text{ref},1}), k_{\text{ref},2}), k_{\text{ref},3}), \mu(\mu_{\text{ref},1}), \phi; z). \quad (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
<figure><img src="image_6.png" /><figcaption>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.</figcaption></figure>

where $\phi_{\log}^B$ and $f_{\text{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

$$P_g(k, z) = D^p(k, z)Z_1^2(k, z)P^{\text{lin}}(k, z) + P_\varepsilon(z), \quad (54)$$

$$B_g(k_1, k_2, k_3, z) = D^B(k_1, k_2, k_3, z) \times \left[ Z_1(k_1, z)Z_1(k_2, z)Z_1(k_3, z)B_I(k_1, k_2, k_3, z) + 2Z_1(k_1, z)Z_1(k_2, z)Z_2(k_1, k_2, z)P_{\text{lin}}(k_1, z)P_{\text{lin}}(k_2, z) + 2 \text{perm} \right] + 2P_{\varepsilon\varepsilon\delta}[Z_1(k_1, z)P_{\text{lin}}(k_1, z) + 2 \text{perm}] + B_\varepsilon, \quad (55)$$

where the linearly propagated primordial oscillatory bispectrum is

$$B_I(k_1, k_2, k_3, z) = M(k_1, z)M(k_2, z)M(k_3, z)B_\Phi^X(k_1, k_2, k_3, z) \quad (56)$$

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

given by

$$Z_1(k_i, z) = b_1(z) + f\mu_i^2, \quad (57)$$

$$Z_2(k_i, k_j, z) = b_1(z)F_2(k_i, k_j) + f(z)\mu_i^2 G_2(k_i, k_j) + \frac{b_2(z)}{2} + \frac{b_{s2}(z)}{2}S_2(k_i, k_j) + \frac{f(z)\mu_{ij}k_{ij}}{2} \left[ \frac{\mu_i}{k_i}Z_1(k_j, z) + \frac{\mu_j}{k_j}Z_1(k_i, z) \right], \quad (58)$$

where $\mu_i = k_i \cdot \hat{z}/k_i$ is the cosine of the angle between the wave vector $k_i$ and the line of sight, $\mu = \mu_1 + \dots + \mu_n$, and $k = |k_1 + \dots + k_n|$. Specifically $\mu_{12} = (\mu_1 k_1 + \mu_2 k_2)/k_{12}$ and $k_{12}^2 = (k_1 + k_2)^2$. 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), while $S_2(k_1, k_2) = (k_1 \cdot 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_{s2}(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. $B_g(k_1, k_2, k_3, z) = B_g(k_1, k_2, k_3, \mu_1, \phi, 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 upcom- ing 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 degen- eracy between the scale-dependent bias coefficient and the pri- mordial bispectrum amplitude ($f_{\text{NL}}^X$), 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\mu(\sigma_p^2(z) + \sigma_f^2(z))]$ and $D^B(k_1, k_2, k_3, z) = \exp[-(k_1^2\mu_1^2 + k_2^2\mu_2^2 + k_3^2\mu_3^2)(\sigma_p^2(z) + \sigma_f^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_\varepsilon = 1/\bar{n}_g$, $P_{\varepsilon\varepsilon\delta} = b_1/(2\bar{n}_g)$, and $B_\varepsilon = 1/\bar{n}_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)

$$B_g^{\text{obs}}(k_{\text{ref},1}, k_{\text{ref},2}, k_{\text{ref},3}, \mu_{\text{ref},1}, \phi; z) = \frac{1}{q_{\parallel}^2 q_{\perp}^4} B_g(k_{\text{ref},1}), k_{\text{ref},2}), k_{\text{ref},3}), \mu(\mu_{\text{ref},1}), \phi; z). \quad (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

$$F_{\alpha\beta}^{B}(z_i) = \frac{1}{4\pi} \sum_{k_1,k_2,k_3} \int_{-1}^{1} d\mu_1 \int_{0}^{2\pi} d\phi \frac{1}{\Delta B^2(k_1, k_2, k_3, z_i)} \times \frac{\partial B_{g}^{\text{obs}}(k_1, k_2, k_3, z_i) \partial B_{g}^{\text{obs}}(k_1, k_2, k_3, z_i)}{\partial p_{\alpha} \partial p_{\beta}}, \quad (60)$$

where the sum over triangles has $k_{min} \le k_3 \le k_2 \le k_1 \le k_{max}$, and $k_1, k_2$ and $k_3$ satisfy the triangle inequality. The bin size $\Delta k$ is taken to be the fundamental frequency of the survey, $k_f = 2\pi/L$, where for simplicity we approximate the survey volume as a cube, $L = V_{\text{survey}}^{1/3}$. The $k$ mode is binned with a bin size of $\Delta k$, and here we consider $\Delta 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<fnref n="7" />. The variance for the bispectrum is then (Sefusatti et al. 2007)

$$\Delta B^2(k_1, k_2, k_3, z) = s_{123} \pi k_f(z)^3 \frac{P_g(k_1, z) P_g(k_2, z) P_g(k_3, z)}{k_1 k_2 k_3 \Delta k(z)^3}, \quad (61)$$

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

$$\Theta_{\text{final}}^{B} = \{\Omega_{m,0}, \Omega_{b,0}, h, n_s, \sigma_8, A_X, \omega_X, \phi_X, f_{\text{NL}}^{X}, \phi_{\text{NL}}^{B}\}. \quad (62)$$

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 $\Theta_{\text{final}}^B$. The forecasted marginalised $1\sigma$ 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_{\text{lin}}(k_i)P_{\text{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_{\text{NL}}^X$ and $\phi_X^B$ feature parameters is the primordial component $B_I$ of the galaxy bispectrum, justifying the less rigorous forecasts presented here.

## 7. 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 (LiteBIRD Collaboration 2023), will be able to reduce the uncertainties on the PPS, in particular in the presence of oscillatory

<table><caption>**Table 3.** Fisher-forecast 68.3% CL marginalised uncertainties on primordial feature parameters, relative to their corresponding fiducial values, using Euclid GC<sub>sp</sub> bispectrum.</caption><thead><tr><th></th><th>P</th><th>B</th><th>P+B</th></tr></thead><tbody><tr><td colspan="4">Pessimistic setting (k<sub>max</sub> = 0.25 h Mpc<sup>-1</sup>)</td></tr><tr><td>A<sub>lin</sub></td><td>23%</td><td>16%</td><td>13%</td></tr><tr><td>ω<sub>lin</sub></td><td>4.1%</td><td>2.1%</td><td>1.7%</td></tr><tr><td>φ<sub>lin</sub></td><td>0.11</td><td>—</td><td>0.049</td></tr><tr><td>f<sub>NL</sub><sup>lin</sup></td><td>—</td><td>31</td><td>31</td></tr><tr><td>φ<sub>lin</sub><sup>B</sup></td><td>—</td><td>4.9</td><td>4.9</td></tr><tr><td>A<sub>log</sub></td><td>25%</td><td>17%</td><td>13%</td></tr><tr><td>ω<sub>log</sub></td><td>4.7%</td><td>2.9%</td><td>2.4%</td></tr><tr><td>φ<sub>log</sub></td><td>0.053</td><td>—</td><td>0.026</td></tr><tr><td>f<sub>NL</sub><sup>log</sup></td><td>—</td><td>33</td><td>32</td></tr><tr><td>φ<sub>log</sub><sup>B</sup></td><td>—</td><td>5.2</td><td>5.1</td></tr><tr><td colspan="4">Optimistic setting (k<sub>max</sub> = 0.3 h Mpc<sup>-1</sup>)</td></tr><tr><td>A<sub>lin</sub></td><td>23%</td><td>12%</td><td>11%</td></tr><tr><td>ω<sub>lin</sub></td><td>3.6%</td><td>1.4%</td><td>1.2%</td></tr><tr><td>φ<sub>lin</sub><sup>B</sup></td><td>0.10</td><td>—</td><td>0.039</td></tr><tr><td>f<sub>NL</sub><sup>lin</sup></td><td>—</td><td>29</td><td>29</td></tr><tr><td>φ<sub>lin</sub><sup>B</sup></td><td>—</td><td>4.6</td><td>4.6</td></tr><tr><td>A<sub>log</sub></td><td>25%</td><td>13%</td><td>11%</td></tr><tr><td>ω<sub>log</sub></td><td>4.7%</td><td>2.2%</td><td>2.0%</td></tr><tr><td>φ<sub>log</sub></td><td>0.052</td><td>—</td><td>0.023</td></tr><tr><td>f<sub>NL</sub><sup>log</sup></td><td>—</td><td>30</td><td>30</td></tr><tr><td>φ<sub>log</sub><sup>B</sup></td><td>—</td><td>4.8</td><td>4.8</td></tr></tbody></table>

**Notes.** We show results for LIN and LOG models in the pessimistic and optimistic settings in the case, using the Euclid GC<sub>sp</sub> bispectrum and the joint power spectrum+bispectrum signal.

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), Miranda et al. (2015), Ballardini et al. (2016), Finelli et al. (2018), Hazra et al. (2018), Braglia et al. (2023).

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 (2022a) for a complete characterisation and quantification of the importance of cross-correlation between CMB and Euclid measurements. Following Euclid Collaboration (2022a), 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)

$$N_{\ell}^{X} = w_{\chi}^{-1} b_{\ell}^{-2}, \quad b_{\ell} = \exp[-\ell(\ell + 1)\theta_{\text{FWHM}}^{2}/(16 \ln 2)], \quad (63)$$

where $\theta_{\text{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<fnref n="8" />; see Okamoto & Hu (2003).

<footnote n="8">⁸ https://github.com/dhanson/quicklens</footnote>

<footnote n="7">⁷ Neglecting the off-diagonal terms in the bispectrum covariance can have a significant impact on the constraints on primordial features, since these contributions are important even on large scales for squeezed triangles (Gualdi & Verde 2020; Biagetti et al. 2022; Flöss et al. 2023). Considering the analytic expression of the full bispectrum covariance is left for future work.</footnote>
matrix for one redshift bin $z_i$ is given by

$$F_{\alpha\beta}^{B}(z_i) = \frac{1}{4\pi} \sum_{k_1,k_2,k_3} \int_{-1}^{1} d\mu_1 \int_{0}^{2\pi} d\phi \frac{1}{\Delta B^2(k_1, k_2, k_3, z_i)} \times \frac{\partial B_{g}^{\text{obs}}(k_1, k_2, k_3, z_i) \partial B_{g}^{\text{obs}}(k_1, k_2, k_3, z_i)}{\partial p_{\alpha} \partial p_{\beta}}, \quad (60)$$

where the sum over triangles has $k_{min} \le k_3 \le k_2 \le k_1 \le k_{max}$, and $k_1, k_2$ and $k_3$ satisfy the triangle inequality. The bin size $\Delta k$ is taken to be the fundamental frequency of the survey, $k_f = 2\pi/L$, where for simplicity we approximate the survey volume as a cube, $L = V_{\text{survey}}^{1/3}$. The $k$ mode is binned with a bin size of $\Delta k$, and here we consider $\Delta 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<fnref n="7" />. The variance for the bispectrum is then (Sefusatti et al. 2007)

$$\Delta B^2(k_1, k_2, k_3, z) = s_{123} \pi k_f(z)^3 \frac{P_g(k_1, z) P_g(k_2, z) P_g(k_3, z)}{k_1 k_2 k_3 \Delta k(z)^3}, \quad (61)$$

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

$$\Theta_{\text{final}}^{B} = \{\Omega_{m,0}, \Omega_{b,0}, h, n_s, \sigma_8, A_X, \omega_X, \phi_X, f_{\text{NL}}^{X}, \phi_{\text{NL}}^{B}\}. \quad (62)$$

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 $\Theta_{\text{final}}^B$. The forecasted marginalised $1\sigma$ 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_{\text{lin}}(k_i)P_{\text{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_{\text{NL}}^X$ and $\phi_X^B$ feature parameters is the primordial component $B_I$ of the galaxy bispectrum, justifying the less rigorous forecasts presented here.

## 7. 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 (LiteBIRD Collaboration 2023), will be able to reduce the uncertainties on the PPS, in particular in the presence of oscillatory

<table><caption>**Table 3.** Fisher-forecast 68.3% CL marginalised uncertainties on primordial feature parameters, relative to their corresponding fiducial values, using Euclid GC<sub>sp</sub> bispectrum.</caption><thead><tr><th></th><th>P</th><th>B</th><th>P+B</th></tr></thead><tbody><tr><td colspan="4">Pessimistic setting (k<sub>max</sub> = 0.25 h Mpc<sup>-1</sup>)</td></tr><tr><td>A<sub>lin</sub></td><td>23%</td><td>16%</td><td>13%</td></tr><tr><td>ω<sub>lin</sub></td><td>4.1%</td><td>2.1%</td><td>1.7%</td></tr><tr><td>φ<sub>lin</sub></td><td>0.11</td><td>—</td><td>0.049</td></tr><tr><td>f<sub>NL</sub><sup>lin</sup></td><td>—</td><td>31</td><td>31</td></tr><tr><td>φ<sub>lin</sub><sup>B</sup></td><td>—</td><td>4.9</td><td>4.9</td></tr><tr><td>A<sub>log</sub></td><td>25%</td><td>17%</td><td>13%</td></tr><tr><td>ω<sub>log</sub></td><td>4.7%</td><td>2.9%</td><td>2.4%</td></tr><tr><td>φ<sub>log</sub></td><td>0.053</td><td>—</td><td>0.026</td></tr><tr><td>f<sub>NL</sub><sup>log</sup></td><td>—</td><td>33</td><td>32</td></tr><tr><td>φ<sub>log</sub><sup>B</sup></td><td>—</td><td>5.2</td><td>5.1</td></tr><tr><td colspan="4">Optimistic setting (k<sub>max</sub> = 0.3 h Mpc<sup>-1</sup>)</td></tr><tr><td>A<sub>lin</sub></td><td>23%</td><td>12%</td><td>11%</td></tr><tr><td>ω<sub>lin</sub></td><td>3.6%</td><td>1.4%</td><td>1.2%</td></tr><tr><td>φ<sub>lin</sub><sup>B</sup></td><td>0.10</td><td>—</td><td>0.039</td></tr><tr><td>f<sub>NL</sub><sup>lin</sup></td><td>—</td><td>29</td><td>29</td></tr><tr><td>φ<sub>lin</sub><sup>B</sup></td><td>—</td><td>4.6</td><td>4.6</td></tr><tr><td>A<sub>log</sub></td><td>25%</td><td>13%</td><td>11%</td></tr><tr><td>ω<sub>log</sub></td><td>4.7%</td><td>2.2%</td><td>2.0%</td></tr><tr><td>φ<sub>log</sub></td><td>0.052</td><td>—</td><td>0.023</td></tr><tr><td>f<sub>NL</sub><sup>log</sup></td><td>—</td><td>30</td><td>30</td></tr><tr><td>φ<sub>log</sub><sup>B</sup></td><td>—</td><td>4.8</td><td>4.8</td></tr></tbody></table>

**Notes.** We show results for LIN and LOG models in the pessimistic and optimistic settings in the case, using the Euclid GC<sub>sp</sub> bispectrum and the joint power spectrum+bispectrum signal.

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), Miranda et al. (2015), Ballardini et al. (2016), Finelli et al. (2018), Hazra et al. (2018), Braglia et al. (2023).

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 (2022a) for a complete characterisation and quantification of the importance of cross-correlation between CMB and Euclid measurements. Following Euclid Collaboration (2022a), 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)

$$N_{\ell}^{X} = w_{\chi}^{-1} b_{\ell}^{-2}, \quad b_{\ell} = \exp[-\ell(\ell + 1)\theta_{\text{FWHM}}^{2}/(16 \ln 2)], \quad (63)$$

where $\theta_{\text{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<fnref n="8" />; see Okamoto & Hu (2003).

<footnote n="8">⁸ https://github.com/dhanson/quicklens</footnote>

<footnote n="7">⁷ Neglecting the off-diagonal terms in the bispectrum covariance can have a significant impact on the constraints on primordial features, since these contributions are important even on large scales for squeezed triangles (Gualdi & Verde 2020; Biagetti et al. 2022; Flöss et al. 2023). Considering the analytic expression of the full bispectrum covariance is left for future work.</footnote>
<table><caption>**Table 4.** Fisher-matrix-based forecasted 68.3% CL marginalised uncertainties using Euclid in combination with the CMB.</caption><thead><tr><th></th><th>Euclid<br>+ Planck-like</th><th>Euclid<br>+ SO-like<br>+ Planck low-ℓ</th><th>Euclid<br>+ CMB-S4<br>+ Planck low-ℓ</th></tr></thead><tbody><tr><td colspan="4"><strong>Pessimistic setting</strong></td></tr><tr><td>A<sub>lin</sub></td><td>18%</td><td>14%</td><td>12%</td></tr><tr><td>ω<sub>lin</sub></td><td>2.4%</td><td>1.2%</td><td>1.0%</td></tr><tr><td>φ<sub>lin</sub></td><td>0.066</td><td>0.046</td><td>0.043</td></tr><tr><td>A<sub>log</sub></td><td>17%</td><td>11%</td><td>7.9%</td></tr><tr><td>ω<sub>log</sub></td><td>1.9%</td><td>1.6%</td><td>1.5%</td></tr><tr><td>φ<sub>log</sub></td><td>0.048</td><td>0.036</td><td>0.033</td></tr><tr><td colspan="4"><strong>Optimistic setting</strong></td></tr><tr><td>A<sub>lin</sub></td><td>16%</td><td>13%</td><td>12%</td></tr><tr><td>ω<sub>lin</sub></td><td>1.0%</td><td>0.92%</td><td>0.85%</td></tr><tr><td>φ<sub>lin</sub></td><td>0.042</td><td>0.038</td><td>0.037</td></tr><tr><td>A<sub>log</sub></td><td>15%</td><td>10%</td><td>7.6%</td></tr><tr><td>ω<sub>log</sub></td><td>1.1%</td><td>1.0%</td><td>0.98%</td></tr><tr><td>φ<sub>log</sub></td><td>0.037</td><td>0.027</td><td>0.024</td></tr></tbody></table>

**Notes.** We show results for LIN and LOG models in the pessimistic and optimistic settings in the case, relative to their corresponding fiducial values ($\mathcal{A}_X = 0.01$, $\omega_X = 10$), using Euclid ($\text{GC}_{\text{sp}}+\text{WL}+\text{GC}_{\text{ph}}+\text{XC}$) in combination with the CMB.

For the Planck-like experiment, we use an effective sensitivity in order to reproduce the Planck 2018 results for the ΛCDM model (Planck Collaboration VI 2020) avoiding the complexity of the real-data likelihood (Planck Collaboration V 2020). We use noise specifications corresponding to in-flight performances of the 143 GHz channel of the high-frequency instrument (HFI; Planck Collaboration I 2020) with a sky fraction $f_{sky} = 0.7$ and a multipole range for temperature and polarisation from $\ell_{min} = 2$ to $\ell_{max} = 1500$. The E-mode polarisation noise is inflated by a factor of 8 to reproduce the uncertainty of the optical depth parameter $\tau$; 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 $\tau$ 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<fnref n="9" />. We use spectra from $\ell_{min} = 40$ to $\ell_{max} = 3000$ for temperature and temperature-polarisation cross-correlation, and $\ell_{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\mu\text{K arcmin}$ in temperature and $\sqrt{2}\mu\text{K arcmin}$ in polarisation, with resolution of $\theta_{\text{FWHM}} = 1$ arcmin. We assume data over the same multipole ranges of SO and with the same sky coverage.

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 Euclid+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 $\omega_X = 10$.

## 8. 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 $\mathcal{A}_X = 0.01$, frequency $\omega_X = 10$, and phase $\phi_X = 0$, where X = {lin, log}.

Following previous studies for ΛCDM and simple extensions in Euclid Collaboration (2020), we have calculated Fisher-matrix-based uncertainties from Euclid primary probes, i.e. spectroscopic galaxy clustering ($\text{GC}_{\text{sp}}$) and the combination of photometric galaxy clustering ($\text{GC}_{\text{ph}}$) and photometric weak lensing (WL); see Sect. 3. We have considered two sets of Euclid specifications: a pessimistic setting with $k_{\text{max}} = 0.25\,\text{h Mpc}^{-1}$ for $\text{GC}_{\text{sp}}$, $\ell_{\text{max}} = 1500$ for WL, $\ell_{\text{max}} = 750$ for $\text{GC}_{\text{ph}}$ and XC, and a cut in redshift of $z < 0.9$ to $\text{GC}_{\text{ph}}$; an optimistic setting with $k_{\text{max}} = 0.30\,\text{h Mpc}^{-1}$ for $\text{GC}_{\text{sp}}$, $\ell_{\text{max}} = 5000$ for WL, and $\ell_{\text{max}} = 3000$ for $\text{GC}_{\text{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 $\text{GC}_{\text{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. $\text{GC}_{\text{sp}}+\text{WL}+\text{GC}_{\text{ph}}+\text{XC}$, we have found that Euclid alone is able to constrain the amplitude of a primordial oscillatory feature with frequency $\omega_X = 10$ to $\mathcal{A}_X = 0.010 \pm 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.003 from 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 $\mathcal{A}_{\text{lin}} = 0.0100 \pm 0.0014$ ($\mathcal{A}_{\text{log}} = 0.0100 \pm 0.0015$) in the pessimistic setting and $\mathcal{A}_{\text{lin}} = 0.0100 \pm 0.0012$ ($\mathcal{A}_{\text{log}} = 0.0100 \pm 0.0012$) in the optimistic setting, both at 68.3% CL, from $\text{GC}_{\text{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 $\mathcal{A}_{\text{lin}} = 0.0100 \pm 0.0009$ ($\mathcal{A}_{\text{log}} = 0.0100 \pm 0.0009$) in the pessimistic setting and $\mathcal{A}_{\text{lin}} = 0.0100 \pm 0.0008$ ($\mathcal{A}_{\text{log}} = 0.0100 \pm 0.0008$) in the optimistic setting, both at 68.3% CL, from $\text{GC}_{\text{sp}}$ (rec)+WL+GC$_{ph}$+XC+BS.

<footnote n="9">⁹ We use version 3.1.0 available at https://github.com/simonsobs/so-noise_models.</footnote>
<table><caption>**Table 4.** Fisher-matrix-based forecasted 68.3% CL marginalised uncertainties using Euclid in combination with the CMB.</caption><thead><tr><th></th><th>Euclid<br>+ Planck-like</th><th>Euclid<br>+ SO-like<br>+ Planck low-ℓ</th><th>Euclid<br>+ CMB-S4<br>+ Planck low-ℓ</th></tr></thead><tbody><tr><td colspan="4"><strong>Pessimistic setting</strong></td></tr><tr><td>A<sub>lin</sub></td><td>18%</td><td>14%</td><td>12%</td></tr><tr><td>ω<sub>lin</sub></td><td>2.4%</td><td>1.2%</td><td>1.0%</td></tr><tr><td>φ<sub>lin</sub></td><td>0.066</td><td>0.046</td><td>0.043</td></tr><tr><td>A<sub>log</sub></td><td>17%</td><td>11%</td><td>7.9%</td></tr><tr><td>ω<sub>log</sub></td><td>1.9%</td><td>1.6%</td><td>1.5%</td></tr><tr><td>φ<sub>log</sub></td><td>0.048</td><td>0.036</td><td>0.033</td></tr><tr><td colspan="4"><strong>Optimistic setting</strong></td></tr><tr><td>A<sub>lin</sub></td><td>16%</td><td>13%</td><td>12%</td></tr><tr><td>ω<sub>lin</sub></td><td>1.0%</td><td>0.92%</td><td>0.85%</td></tr><tr><td>φ<sub>lin</sub></td><td>0.042</td><td>0.038</td><td>0.037</td></tr><tr><td>A<sub>log</sub></td><td>15%</td><td>10%</td><td>7.6%</td></tr><tr><td>ω<sub>log</sub></td><td>1.1%</td><td>1.0%</td><td>0.98%</td></tr><tr><td>φ<sub>log</sub></td><td>0.037</td><td>0.027</td><td>0.024</td></tr></tbody></table>

**Notes.** We show results for LIN and LOG models in the pessimistic and optimistic settings in the case, relative to their corresponding fiducial values ($\mathcal{A}_X = 0.01$, $\omega_X = 10$), using Euclid ($GC_{sp}+WL+GC_{ph}+XC$) in combination with the CMB.

For the Planck-like experiment, we use an effective sensitivity in order to reproduce the Planck 2018 results for the ΛCDM model (Planck Collaboration VI 2020) avoiding the complexity of the real-data likelihood (Planck Collaboration V 2020). We use noise specifications corresponding to in-flight performances of the 143 GHz channel of the high-frequency instrument (HFI; Planck Collaboration I 2020) with a sky fraction $f_{sky} = 0.7$ and a multipole range for temperature and polarisation from $\ell_{min} = 2$ to $\ell_{max} = 1500$. The E-mode polarisation noise is inflated by a factor of 8 to reproduce the uncertainty of the optical depth parameter $\tau$; 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 $\tau$ 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<fnref n="9" />. We use spectra from $\ell_{min} = 40$ to $\ell_{max} = 3000$ for temperature and temperature-polarisation cross-correlation, and $\ell_{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\mu K$ arcmin in temperature and $\sqrt{2}\mu K$ arcmin in polarisation, with resolution of $\theta_{FWHM} = 1$ arcmin. We assume data over the same multipole ranges of SO and with the same sky coverage.

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 Euclid+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 $\omega_X = 10$.

## 8. 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 $\mathcal{A}_X = 0.01$, frequency $\omega_X = 10$, and phase $\phi_X = 0$, where X = {lin, log}.

Following previous studies for ΛCDM and simple extensions in Euclid Collaboration (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}$, $\ell_{max} = 1500$ for WL, $\ell_{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}$, $\ell_{max} = 5000$ for WL, and $\ell_{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 $\omega_X = 10$ to $\mathcal{A}_X = 0.010 \pm 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.003 from 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 $\mathcal{A}_{lin} = 0.0100 \pm 0.0014$ ($\mathcal{A}_{log} = 0.0100 \pm 0.0015$) in the pessimistic setting and $\mathcal{A}_{lin} = 0.0100 \pm 0.0012$ ($\mathcal{A}_{log} = 0.0100 \pm 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 $\mathcal{A}_{lin} = 0.0100 \pm 0.0009$ ($\mathcal{A}_{log} = 0.0100 \pm 0.0009$) in the pessimistic setting and $\mathcal{A}_{lin} = 0.0100 \pm 0.0008$ ($\mathcal{A}_{log} = 0.0100 \pm 0.0008$) in the optimistic setting, both at 68.3% CL, from $GC_{sp}$ (rec)+WL+GC_{ph}+XC+BS.

<footnote n="9">⁹ We use version 3.1.0 available at https://github.com/simonsobs/so.noise.models.</footnote>
<figure><img src="image_7.png" /><figcaption>**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 $\omega_{\text{lin}} = 10$ (top panels) and the LOG model with $\omega_{\text{log}} = 10$ (bottom panels). Left (right) panels correspond to the pessimistic (optimistic) setting for $\text{GC}_{\text{sp}}+\text{WL}+\text{GC}_{\text{ph}}+\text{XC}$ (red), $\text{GC}_{\text{sp}}+\text{WL}+\text{GC}_{\text{ph}}+\text{XC}$ with numerical reconstruction of $\text{GC}_{\text{sp}}$ (green), $\text{GC}_{\text{sp}}(\text{rec})+\text{WL}+\text{GC}_{\text{ph}}+\text{XC}$ in combination with Euclid $\text{GC}_{\text{sp}}$ bispectrum (orange), and their combination plus SO-like CMB (blue).</figcaption></figure>

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 (2022a); see Sect. 7. Combining the Fisher matrix information from a SO-like experiment (complemented with a Planck-like one at low-$\ell$) with the Euclid primary probes, we have found $\mathcal{A}_{\text{lin}} = 0.0100 \pm 0.0014$ ($\mathcal{A}_{\text{log}} = 0.0100 \pm 0.0011$) in the pessimistic setting and $\mathcal{A}_{\text{lin}} = 0.0100 \pm 0.0013$ ($\mathcal{A}_{\text{log}} = 0.0100 \pm 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-$\ell$, correspond to $\mathcal{A}_{\text{lin}} = 0.0100 \pm 0.0008$ at a 68.3% CL and to $\mathcal{A}_{\text{log}} = 0.0100 \pm 0.0008$ for $\text{GC}_{\text{sp}}(\text{PS rec} + \text{BS})+\text{WL}+\text{GC}_{\text{ph}}+\text{XC}+\text{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 $\mathcal{A}_\chi$. 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 VI 2020): $\mathcal{A}_{\text{lin}} = 0.015$ and
<figure><img src="image_7.png" /><figcaption>**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 $\omega_{\text{lin}} = 10$ (top panels) and the LOG model with $\omega_{\text{log}} = 10$ (bottom panels). Left (right) panels correspond to the pessimistic (optimistic) setting for $\text{GC}_{\text{sp}}+\text{WL}+\text{GC}_{\text{ph}}+\text{XC}$ (red), $\text{GC}_{\text{sp}}+\text{WL}+\text{GC}_{\text{ph}}+\text{XC}$ with numerical reconstruction of $\text{GC}_{\text{sp}}$ (green), $\text{GC}_{\text{sp}}(\text{rec})+\text{WL}+\text{GC}_{\text{ph}}+\text{XC}$ in combination with Euclid $\text{GC}_{\text{sp}}$ bispectrum (orange), and their combination plus SO-like CMB (blue).</figcaption></figure>

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 (2022a); see Sect. 7. Combining the Fisher matrix information from a SO-like experiment (complemented with a Planck-like one at low-$\ell$) with the Euclid primary probes, we have found $\mathcal{A}_{\text{lin}} = 0.0100 \pm 0.0014$ ($\mathcal{A}_{\text{log}} = 0.0100 \pm 0.0011$) in the pessimistic setting and $\mathcal{A}_{\text{lin}} = 0.0100 \pm 0.0013$ ($\mathcal{A}_{\text{log}} = 0.0100 \pm 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-$\ell$, correspond to $\mathcal{A}_{\text{lin}} = 0.0100 \pm 0.0008$ at a 68.3% CL and to $\mathcal{A}_{\text{log}} = 0.0100 \pm 0.0008$ for $\text{GC}_{\text{sp}}(\text{PS rec} + \text{BS})+\text{WL}+\text{GC}_{\text{ph}}+\text{XC}+\text{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 $\mathcal{A}_\chi$. 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 VI 2020): $\mathcal{A}_{\text{lin}} = 0.015$ and