# Reconstructing Primordial Curvature Perturbations via Scalar-Induced Gravitational Waves with LISA

Jonas El Gammal, Aya Ghaleb, Gabriele Franciolini, Theodoros

Papanikolaou, Marco Peloso, Gabriele Perna, Mauro Pieroni, Angelo

Ricciardone, Robert Rosati, Gianmassimo Tasinato, et al.

► To cite this version:

Jonas El Gammal, Aya Ghaleb, Gabriele Franciolini, Theodoros Papanikolaou, Marco Peloso, et al.. Reconstructing Primordial Curvature Perturbations via Scalar-Induced Gravitational Waves with LISA. JCAP, 2025, 05, pp.062. 10.1088/1475-7516/2025/05/062. hal-04928410

HAL Id: hal-04928410

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

Submitted on 31 Oct 2025

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
PAPER • OPEN ACCESS

# Reconstructing primordial curvature perturbations via scalar-induced gravitational waves with LISA

To cite this article: Jonas El Gammal et al JCAP05(2025)062

View the article online for updates and enhancements.

You may also like

- Unraveling the early universe's equation of state and primordial black hole production with PTA, BBN, and CMB observations
  Qing-Hua Zhu, , Zhi-Chao Zhao et al.

- Constraints on Primordial Black Holes from $N_{eff}$: Scalar-induced Gravitational Waves as an Extra Radiation Component
  Jing-Zhi Zhou, Yu-Ting Kuang, Zhe Chang et al.

- Gauge-dependence of scalar induced gravitational waves
  Anjali Abirami Kugarajh
RECEIVED: January 28, 2025

REVISED: March 29, 2025

ACCEPTED: April 14, 2025

PUBLISHED: May 23, 2025

# Reconstructing primordial curvature perturbations via scalar-induced gravitational waves with LISA

Jonas El Gammal<sup>1, a, *</sup> Aya Ghaleb<sup>1, b</sup> Gabriele Franciolini<sup>1, c, 1</sup>

Theodoros Papanikolaou<sup>1, d, e, f</sup> Marco Peloso<sup>1, g, h</sup> Gabriele Perna<sup>1, g, h, *</sup>

Mauro Pieroni<sup>1, c</sup> Angelo Ricciardone<sup>1, i, j</sup> Robert Rosati<sup>1, k, 1, 2</sup>

Gianmassimo Tasinato<sup>1, b, l, m</sup> Matteo Braglia<sup>1, n</sup> Jacopo Fumagalli<sup>1, o</sup>

Jun'ya Kume<sup>1, g, h, r</sup> Enrico Morgante<sup>1, p, q</sup> Germano Nardini<sup>1, a</sup> Davide Racco<sup>1, s, t</sup>

Sébastien Renaux-Petel<sup>1, u</sup> Hardi Veermäe<sup>1, v</sup> Denis Werth<sup>1, u</sup> and Ivonne Zavala<sup>1, b</sup>
for the LISA Cosmology Working Group

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

<footnote><sup>a</sup>Department of Mathematics and Physics, University of Stavanger, NO-4036 Stavanger, Norway</footnote>

<footnote><sup>b</sup>Department of Physics, Faculty of Science and Engineering, Swansea University,
Singleton Park, SA2 8PP, Swansea, U.K.</footnote>

<footnote><sup>c</sup>CERN, Theoretical Physics Department, Esplanade des Particules 1, Geneva 1211, Switzerland</footnote>

<footnote><sup>d</sup>Scuola Superiore Meridionale, Largo San Marcellino 10, 80138 Napoli, Italy</footnote>

<footnote><sup>e</sup>Istituto Nazionale di Fisica Nucleare (INFN), Sezione di Napoli,
Via Cinthia 21, 80126 Napoli, Italy</footnote>

<footnote><sup>f</sup>National Observatory of Athens, Lofos Nymfon, 11852 Athens, Greece</footnote>

<footnote><sup>g</sup>Dipartimento di Fisica e Astronomia “G. Galilei”, Università degli Studi di Padova,
via Marzolo 8, I-35131 Padova, Italy</footnote>

<footnote><sup>h</sup>INFN, Sezione di Padova, via Marzolo 8, I-35131 Padova, Italy</footnote>

<footnote><sup>i</sup>Dipartimento di Fisica “Enrico Fermi”, Università di Pisa,
Largo Bruno Pontecorvo 3, Pisa I-56127, Italy</footnote>

<footnote><sup>j</sup>INFN, Sezione di Pisa, Largo Bruno Pontecorvo 3, Pisa I-56127, Italy</footnote>

<footnote><sup>k</sup>NASA Marshall Space Flight Center, Huntsville, AL 35812, U.S.A.</footnote>

<footnote><sup>l</sup>Dipartimento di Fisica e Astronomia, Università di Bologna, Bologna, Italy</footnote>

<footnote><sup>m</sup>INFN, Sezione di Bologna, I.S. FLAG, viale B. Pichat 6/2, 40127 Bologna, Italy</footnote>

<footnote><sup>n</sup>Center for Cosmology and Particle Physics, New York University,
726 Broadway, New York, NY 10003, U.S.A.</footnote>

<footnote><sup>o</sup>Departement de Física Quàntica i Astrofísica and Institut de Ciències del Cosmos (ICC),
Universitat de Barcelona, Martí i Franquès 1, 08028 Barcelona, Spain</footnote>

<footnote><sup>p</sup>Dipartimento di Fisica, Università di Trieste, Strada Costiera 11, I-34151 Trieste, Italy</footnote>

<footnote>*Corresponding author.</footnote>

<footnote n="1">¹Project coordinator.</footnote>

<footnote n="2">²NASA Postdoctoral Program Fellow.</footnote>
qINFN, Sezione di Trieste, Via Valerio 2, 34127 Trieste, Italy

†Research Center for the Early Universe (RESCEU), Graduate School of Science, The University of Tokyo, Hongo 7-3-1 Bunkyo-ku, Tokyo 113-0033, Japan

sInstitut für Theoretische Physik, ETH Zürich, Wolfgang-Pauli-Str. 27, 8093 Zürich, Switzerland

tPhysik-Institut, Universität Zürich, Winterthurerstrasse 190, 8057 Zürich, Switzerland

uInstitut d'Astrophysique de Paris, UMR 7095 du CNRS et de Sorbonne Université, 98 bis bd Arago, 75014 Paris, France

vKeemilise ja bioloogilise füüsika instituut, Rävala pst. 10, 10143 Tallinn, Estonia

E-mail: jonas.elgammal@hotmail.com, gabriele.franciolini@cern.ch, gabriele.perna@phd.unipd.it, robert.j.rosati@nasa.gov

**ABSTRACT:** Many early universe scenarios predict an enhancement of scalar perturbations at scales currently unconstrained by cosmological probes. These perturbations source gravitational waves (GWs) at second order in perturbation theory, leading to a scalar-induced gravitational wave (SIGW) background. The LISA detector, sensitive to mHz GWs, will be able to constrain curvature perturbations in a new window corresponding to scales $k \in [10^{10}, 10^{14}]$ Mpc$^{-1}$, difficult to probe otherwise. In this work, we forecast the capabilities of LISA to constrain the source of SIGWs using different approaches: i) agnostic, where the spectrum of curvature perturbations is binned in frequency space; ii) template-based, modeling the curvature power spectrum based on motivated classes of models; iii) ab initio, starting from first-principles model of inflation featuring an ultra-slow roll phase. We compare the strengths and weaknesses of each approach. We also discuss the impact on the SIGW spectrum of non-standard thermal histories affecting the kernels of SIGW emission and non-Gaussianity in the statistics of the curvature perturbations. Finally, we propose simple tests to assess whether the signal is compatible with the SIGW hypothesis. The pipeline used is built into the **SIGWAY** code.

**KEYWORDS:** cosmological perturbation theory, gravitational waves / sources, primordial black holes, primordial gravitational waves (theory)

ARXIV EPRINT: 2501.11320
# Contents

<table><tr><td>1</td><td>Introduction</td><td>2</td></tr><tr><td>2</td><td>Early universe models leading to enhanced curvature spectra</td><td>6</td></tr><tr><td>2.1</td><td>Single field inflation</td><td>6</td></tr><tr><td>2.2</td><td>Multi-field inflation</td><td>7</td></tr><tr><td>2.3</td><td>Other classes of models</td><td>7</td></tr><tr><td>2.4</td><td>A simple benchmark scenario: single field USR</td><td>8</td></tr><tr><td>3</td><td>Modeling the curvature power spectrum</td><td>10</td></tr><tr><td>3.1</td><td>Model independent parameterization: binned spectrum approach</td><td>11</td></tr><tr><td>3.2</td><td>Analytical templates for curvature spectra</td><td>11</td></tr><tr><td>3.3</td><td>Computation of P<sub>ζ</sub> in single field USR scenarios</td><td>15</td></tr><tr><td>4</td><td>Computation of the scalar-induced GW background</td><td>19</td></tr><tr><td>4.1</td><td>Source of GWs at second order in the scalar perturbations</td><td>20</td></tr><tr><td>4.2</td><td>Radiation-dominated era</td><td>22</td></tr><tr><td>4.3</td><td>Transition from an early matter-dominated to the radiation-dominated era</td><td>22</td></tr><tr><td>4.4</td><td>Computation of SIGW with the binned spectrum approach</td><td>23</td></tr><tr><td>4.5</td><td>Non-Gaussian imprints on the SIGW spectrum</td><td>24</td></tr><tr><td>5</td><td>Mock signal reconstructions with the SGWBinner and SIGWAY codes</td><td>27</td></tr><tr><td>5.1</td><td>Data streams from LISA TDI channels</td><td>28</td></tr><tr><td>5.2</td><td>Analysis of the simulated data</td><td>31</td></tr><tr><td>6</td><td>Results</td><td>33</td></tr><tr><td>6.1</td><td>Binned spectrum method</td><td>34</td></tr><tr><td>6.2</td><td>Template based method</td><td>35</td></tr><tr><td>6.3</td><td>Single field USR inference</td><td>45</td></tr><tr><td>6.4</td><td>Non standard thermal histories</td><td>49</td></tr><tr><td>6.5</td><td>Non-Gaussian effects on SIGWs</td><td>51</td></tr><tr><td>7</td><td>Testing the scalar-induced hypothesis</td><td>55</td></tr><tr><td>8</td><td>Conclusions</td><td>59</td></tr><tr><td>A</td><td>SIGWAY code: technicalities</td><td>62</td></tr><tr><td>A.1</td><td>Perturbations in USR scenarios: code structure</td><td>62</td></tr><tr><td>A.2</td><td>Computation of SIGWs from the spectrum of curvature perturbations</td><td>63</td></tr><tr><td>A.3</td><td>Computation of SIGWs using binned coefficients</td><td>64</td></tr><tr><td>A.4</td><td>Computation of SIGWs including primordial NGs</td><td>65</td></tr><tr><td>A.5</td><td>Inference</td><td>66</td></tr></table>
B Challenges with binned analyses and a large number of bins 67

C Testing the resolvability of non-Gaussian corrections: additional plots 70

D Testing the scalar-induced hypothesis: additional plots 70

# 1 Introduction

The Laser Interferometer Space Antenna (LISA) [1] represents a groundbreaking gravitational wave (GW) observatory aimed to probe and impact our understanding of fundamental physics, astronomy, and cosmology [2–4]. With the first-ever direct probe of the stochastic gravitational-wave background (SGWB) in the millihertz frequency range, LISA provides the opportunity to unveil processes that occurred in the first stages of the Universe, including inflation [5, 6]. Probing a primordial SGWB at millihertz frequencies corresponds to exploring comoving scales that lie between those accessible to ground-based GW interferometers and those probed by pulsar-timing arrays, cosmic microwave background (CMB), or large-scale structure surveys. These scales correspond to comoving wavenumbers in the range $k \sim [10^{10}, 10^{14}] \text{ Mpc}^{-1}$. Several cosmological models (see e.g. refs. [7–10] for reviews) predict a measurable SGWB at these scales, often without any other complementary distinctive signature, placing LISA in a unique position to test these scenarios.

Inflationary models exhibiting amplified scalar fluctuations are one of the candidates for sourcing an SGWB in the LISA frequency band. Enhanced scalar fluctuations generate scalar-induced gravitational waves (SIGWs) at second order in perturbations [11–20], resulting in a potentially large SGWB. Amplified scalar fluctuations naturally arise in single-field inflationary scenarios with features in the potential like those leading to ultra-slow-roll (USR) phases, multi-field setups, and mechanisms such as preheating or early matter-dominated eras. Intriguingly, the same perturbations that seed the SIGWs can also trigger the formation of primordial black holes (PBHs) [21–25]. SIGWs in the millihertz frequency band arise in correspondence with the asteroidal-mass window for PBHs, a viable candidate for addressing the dark matter puzzle [26–30]. By detecting or setting upper bounds on SIGWs, LISA would not only shed light on the inflationary epoch but also on dark matter and non-astrophysical black hole formation channels. This makes SIGWs a high-gain, well-motivated target for LISA.

SIGWs are also powerful tools to investigate non-Gaussianity (NG) in the early universe since their production is highly sensitive to the statistical properties of scalar curvature fluctuations [31–40]. NG is typically characterized by parameters such as $f_{\text{NL}}$, which quantifies NG at the bispectrum level (three-point correlation function), and $\tau_{\text{NL}}$, which appears in the trispectrum (four-point correlation function). Earlier analysis have focused on contributions to the SIGW involving $f_{\text{NL}}$, [41–52]. More recent studies have extended this analysis to higher-order NG terms [47, 53–55]. Recently, a Fisher forecast analysis for LISA about NG has been performed in [55]. The tensor power spectrum of SIGWs is directly related to the four-point correlation function of the curvature fluctuations. Such a correlation function has both connected and disconnected contributions. While the latter contributes only to the Gaussian SIGW power spectrum, the former is directly linked to the trispectrum through the $\tau_{\text{NL}}$ parameter [51], which is the key observable to constrain NG from SIGWs.
The detection and characterization of the primordial SGWB is one of the most challenging objectives of the LISA mission [56]. LISA is a signal-dominated detector, where a multitude of transient or quasi-monochromatic events overlap in time and frequency with the stochastic superposition of all unresolved astrophysical events and, potentially also with a significant primordial SGWB. Additionally, the stationary component of the instrumental noise can mimic a SGWB to some extent. Completing the LISA science program for the SGWB therefore requires:

i) Determining whether a primordial SGWB is present in the data.

ii) Reconstructing the SGWB frequency shape and, if possible, its statistical properties.

iii) Setting upper limits on cosmological sources of SGWB not supported by the data.

iv) Constraining the parameters of the most likely SGWB source candidates.

ESA and NASA plan to address these tasks through the so-called “global fit”, a data analysis procedure where modules fitting each class of sources (galactic binaries, supermassive black hole binaries, SGWB, etc.) iterate until convergence [56]. Recently, successful prototype global fit analyses became available in the literature [57–59], tested on the *Sangria* LISA Data Challenge (LDC) dataset [60], which contains no primordial SGWB. It is still an open question how the global fit should support a primordial SGWB search, and how the SGWB properties should be represented in the detection catalogs that the space agencies will publish. A recent study [61] has attempted to perform an SGWB search directly on the global fit residual.

In this work, we aim to bridge these gaps by providing elements of the global fit SGWB module useful for the tasks i)–iv) in the presence of a SGWB due to SIGWs. To develop and test our rationale, we work in the limit that all resolvable events have been precisely reconstructed,<fnref n="1" /> leaving us with data containing the stationary component of the noise, the SIGW background, and the foregrounds from the unresolved galactic and extragalactic binaries.

As we focus on SIGW sources, we perform an analysis starting from the properties of the source curvature power spectra $P_ζ(k)$ of the source, instead of the GW energy density $\Omega_{GW}(f)$ generated by these power spectra.

Concerning i) and ii), we prototype a model-agnostic method that reconstructs the power spectrum $P_ζ(k)$ by binning it in frequency space. This approach allows for maximal flexibility in capturing unknown SIGW features. It is however not as agnostic as other generic SGWB searches [62–68] since it requires, by construction, a SIGW source, i.e. a SGWB that can be derived as the proper convolution of a generic $P_ζ(f)$. Due to this additional information, the method is expected to be more sensitive to SIGW signals than

<footnote n="1">¹Although this optimistic working hypothesis may seem unrealistic, it is the correct one to use in a global fit module. All current implementations of the global fit are based on a blocked Gibbs sampling scheme, where each source type is sampled independently, assuming a perfect subtraction of the other source types. By periodically alternating which source type is being sampled over, imperfect source subtraction and source type confusion are properly included and fully modeled in the resulting posterior. When working after the global fit, with one sample of the residual as we assume in this work, these possible degeneracies are not fully modeled. Additionally, if the global fit has experienced a convergence failure (the MCMC is still “burning in”), unmodeled source power may still be in the data and lead to false detections of an SGWB. Ref. [61] studies how this convergence failure affects stochastic background recovery in the available prototype global fit residuals.</footnote>
other fully-agnostic approaches. It can be particularly useful for placing upper bounds on the SIGW amplitude if no signal is detected in the LISA data, or act as a key ingredient for SIGW model selection if a signal is present.

On the other hand, if the computational resources available to the mission allow running the global fit for every SGWB template, several modules for the template-based SIGW reconstruction have been conceived since the first iterations of the global fit.<fnref n="2" /> The advantage of such a possibility is clear: the more signal characterization is included in the search, the higher the sensitivity to that signal. This process also reduces the risk of SGWB misreconstructions that the global fit might absorb into the parameter estimation of other sources. To address the points i) and ii) above within this framework, we collect several well-motivated inflationary models with known $P_ζ(k)$ predictions, we design template classes that effectively parameterize these $P_ζ(k)$, and we prototype the SIGW template-based searches for them.

Accurately performing i) and ii) enables LISA to identify the most favored SIGW models and then proceed with tasks iii) and iv). Accordingly, we implement a prototype data analysis pipeline, choosing the USR inflationary setup as a representative example. In particular, we develop a fast numerical algorithm determining $P_ζ(k)$ once the inflationary model parameters are known. Thanks to its speed, the algorithm allows for rapid likelihood evaluations in the fundamental-parameter space, enabling direct inference on the USR model parameters from GW data.

As a proof of concept, we further perform inference on $P_ζ(k)$ in cosmological scenarios where standard assumptions on Gaussianity of curvature perturbations and on the standard thermal history of the Universe are relaxed. We evaluate the SIGW energy density sourced by non-Gaussian contributions parametrized by $\tau_{NL}$, resorting to the local ansatz emerging from a perturbative expansion of scalar fluctuations. This contribution is known to modify its frequency profile compared to the Gaussian counterpart [43, 44, 48, 51, 54, 55]. Since for models of inflation with local type NG, where the curvature perturbation is dominated by one degree of freedom, $\tau_{NL}$ and $f_{NL}$ are related, we take advantage of such a relation to forecast the ability of LISA of probing NG, focusing the analysis on $f_{NL}$ and discussing its implications for $\tau_{NL}$. It is worth mentioning that, as recently argued by [69], the effects of NG on the SIGW background may not always be accurately captured by an expansion around a Gaussian field. Properly accounting for the full impact of intrinsic non-linearities may significantly suppress or enhance the spectrum compared to the predictions based on the local ansatz. Achieving this would require the development of fully non-perturbative approaches to compute the SIGW spectrum, which are beyond the scope of this work. Finally, we perform some diagnostic tests to assess whether the reconstructed signal is consistent with the SIGW hypothesis. Such tests could help to rule out a scalar-induced origin as a viable explanation for some SGWB spectral shapes.

The core of our numerical analysis is implemented in the **SIGWAY** code.<fnref n="3" /> This stand-alone **Python** code addresses tasks i)–iv) for SIGW signals by offering the following functionalities:

* A fast vectorized numerical integrator for computing the SGWB resulting from any spectrum of primordial curvature fluctuations of modes reentering the Hubble radius during radiation domination or a phase of early matter domination.

<footnote n="2">²See ref. [6] for other template-based reconstructions suitable for inflationary models.</footnote>

<footnote n="3">³https://github.com/jonaselgammal/SIGWAY.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Flowchart describing the inputs and modules of the SIGWAY code developed for this project, alongside the SGWBinner adopted to make forecasts for LISA.</figcaption></figure>

• An integration algorithm for computing the SGWB assuming a binned spectrum of $P_ζ(k)$ for agnostic reconstructions of $P_ζ$.

• Solvers for the background- and perturbation equations of motion for the inflaton in a single-field scenario that can be called by the SIGW integrator starting from the inflaton Lagrangian.

• Capabilities for computing the SGWB including non-Gaussian contributions for a lognormal shape of $P_ζ$.

• Functionality for pairing to the SGWBinner pipeline [63, 64] for computing the LISA likelihood, and performing inference on the parameters governing the primordial curvature fluctuations.

In figure 1 we report a flowchart that schematically describes the inputs and outputs of the SIGWAY code developed for this project, alongside the SGWBinner adopted to make forecasts for LISA.

The paper is organized as follows. In section 2 we review some representative models predicting enhanced power spectra of curvature perturbations. In section 3 we identify functional forms that describe the shapes of the aforementioned spectra in terms of effective parameters. In section 4 we describe the analytic and numerical tools that we implement in the SIGWAY code. The functionality of the SGWBinner code that is relevant for performing inference is briefly described in section 5. Section 6 illustrates, for representative benchmark signals, how the elements built in this work help tackle the key tasks *i)-iv*), including testing whether the SIGW hypothesis is compatible with the putative signal in section 7. Finally, section 8 presents our main conclusions, while appendix A and appendix B discuss technicalities and subtleties regarding the proposed SIGW signal reconstruction and interpretation.

**Notation.** We indicate with $k$ the comoving wavenumber while with $f$ the associated frequency $f = k/2\pi$. For scales relevant to LISA, we translate wavenumbers in units of Hz using $c = \text{Hz Mpc}/(1.023 \times 10^{14})$. For presentation purposes, we differentiate between frequencies and momenta arbitrarily denoting them with units of Hz and $s^{-1}$, respectively. As usually done in the literature, we report the GW spectral energy density $\Omega_{\text{GW}}$ multiplied by the rescaled Hubble rate $h = H_0/(100 \text{ km/s/Mpc})$ squared. Finally, we indicate vectors with bold symbols (e.g. $\vec{x} \equiv x$) while their magnitude with the lower-case letter.
## 2 Early universe models leading to enhanced curvature spectra

In this section, we summarise the main classes of models predicting enhanced curvature power spectra at small scales, which can lead to potential GW signatures in LISA.

### 2.1 Single field inflation

In the simplest models of inflation, a single scalar field known as the inflaton moves gradually down its potential under the influence of Hubble friction, resulting in a slow-roll (SR) phase. The generated fluctuations are nearly scale-invariant, Gaussian, and adiabatic: they freeze on super-Hubble scales, producing a universe that is statistically homogeneous and isotropic [70–74].

In certain models, however, the inflaton potential in the Einstein frame can contain features such as a flat region or a mild bump that causes the field velocity to decrease rapidly in a brief USR<fnref n="4" /> phase [75–80] followed by another SR or a constant-roll phase [81–83].<fnref n="5" /> The shape of the enhanced spectral features is determined by the amplification of the curvature perturbations during the USR phase, as well as by the specifics of the transition into the USR era. If the field accelerates significantly during that transt epoch, it can cause sizeable spectral modulation, and even be the dominant source of amplification. The resulting perturbations deviate significantly from scale invariance, exhibiting the strongest amplification for wavelengths exiting the horizon around the USR era.

Generically, models of this kind can be subdivided into four categories: i) Quasi-inflection points and plateaus [46, 80, 86–117]; ii) Upward [118, 119] or downward [119–122] steps; iii) Models in which the inflaton rolls through a global minimum/double-well potentials [123–128]; and iv) Potentials with stacked features/oscillating potentials [129–135]. It was also suggested that models going beyond a non-minimally coupled inflaton, e.g. within modified gravity theories, can introduce features in terms other than the inflaton potential or the non-minimal coupling [84, 136–151]. See e.g. [152] for a model-building review. Even though the literature on these models is quite vast, many scenarios predict a curvature power spectrum that is enhanced at small scales with similar properties. In particular, most models, especially those in category i), produce a single peak in the power spectrum that is approximately captured by a broken power law. However, models with sharp features can produce spectral oscillations at and after the peak in the power spectrum. Such enhancement mechanisms are mainly in the categories ii) and iii). Non-standard potentials in category iv) can also deviate strongly from this picture.

We finish this section with a note on the theoretical consistency of these enhancing mechanisms. There has been significant debate about the potential impact of loop corrections, induced by the enhanced modes during an USR phase, on long-wavelength scales, with some even challenging the validity of perturbative computations in these scenarios. This discussion is also relevant for the possible interpretation of a stochastic signal as originating from SIGW.

<footnote n="4">⁴We consider a slightly broader definition of USR, which is often characterized by $\eta_H \equiv -\ddot{H}/(2H\dot{H}) = 3$ (or $\eta \equiv \dot{\epsilon}/(H\epsilon) = -6$) and can be realized on a flat plateau. However, models considered in the literature often exhibit a small bump instead, which results in $\eta_H > 3$ due to the curvature of the potential. We will include such deviations in our definition of USR.</footnote>

<footnote n="5">⁵See e.g. [84, 85] for scenarios based on the EFT of inflation.</footnote>
The question of whether these corrections can become sufficiently large to undermine the predictive power of inflationary models related to PBH and SIGW has been also explored in Res. [153-163], while other studies have questioned the very existence of these corrections and proposed an argument against their presence [164-168].

## 2.2 Multi-field inflation

**Hybrid/multi-field inflation.** Given the high number of degrees of freedom within multi-field inflationary setups, we can separate and control more efficiently two stages responsible respectively for the generation of nearly scale-invariant primordial curvature perturbations on large CMB scales, and for enhanced curvature perturbations on small scales which lead to the production of SIGWs. Hybrid/multi-field models of inflation tend to generate slowly-growing lognormal-like peaks in the curvature power spectrum [169-177] while strong deviations from a geodesic trajectory in field space may lead to sharp peaks and features such as spectral oscillations in $P_ζ$ [178-185].

**Curvaton models.** Within curvaton scenarios [186], one can realize setups with the curvaton field being characterized by a steep blue spectrum either due to interactions with the inflaton or other degrees of freedom during inflation [187, 188] or due to a non-trivial kinetic term [189-191]. In particular, within axion-like curvaton setups, the curvaton field is identified with the phase of a complex field whose modulus decreases rapidly during inflation [192-195]. We should also highlight that in any curvaton model the curvature perturbations on small scales originate from non-adiabatic curvaton field fluctuations during inflation, leading to a non-Gaussian probability distribution function for the primordial curvature perturbations [196-199] with important consequences at the level of the SIGW signal [44, 55, 200, 201] (but see also [202]).

**Axion-gauge field coupling.** Enhanced scalar [203] and tensor [203, 204] modes can be produced by gauge fields amplified by their pseudo-scalar $\phi F\tilde{F}$ coupling with a rolling inflaton or spectator [205, 206] axion during inflation. The gauge field amplitude is exponentially sensitive to the axion velocity, thus providing naturally blue signals. These enhanced curvature modes can lead to PBH and SIGW [42, 207-209]. The precise shape of these signals is sensitive to the axion evolution, which is significantly impacted by the backreaction of the amplified gauge fields, which is recently being explored via lattice simulations [210-213].

## 2.3 Other classes of models

**Preheating.** During the preheating phase following inflation, as the inflaton field undergoes coherent oscillations around the minimum of its potential, a striking phenomenon emerges: the resonant amplification of quantum inflaton field fluctuations, which drives particle production [214, 215]. These enhanced quantum fluctuations are accompanied by a resonant amplification of the scalar metric fluctuations (usually quoted as metric preheating [216-219]), or, in other terms, with enhanced curvature perturbations, responsible for the generation of SIGWs [220] and potentially for PBH formation [221-223] (see however [224] for an assessment on the role of non-linearities and anharmonicities). Most studies have focused on multi-field inflationary setups since in such scenarios the enhancement of entropic (isocurvature) fluctuations can give rise to the enhancement of the adiabatic/curvature fluctuations in the
broad resonance regime [225–230]. This leads to a notable amplification of the primordial curvature power spectrum, deviating from the standard scale-invariant behavior at small scales [231–233]. Interestingly, recent works also suggest that a parametric amplification of the curvature perturbations can occur even in the narrow regime in the case of single field inflation [221–224, 234]. We should also note that these gravitational waves are typically peaked at MHz or GHz frequencies far above those of LISA, although some scenarios do allow for a peak in LISA’s range [235].

**Matter bouncing scenarios.** In non-singular matter bouncing cosmological models [236], the matter contracting phase inevitably amplifies super-horizon curvature perturbations. This enhancement can lead to an enhanced primordial curvature power spectrum on small scales compared to the ones probed by CMB. As these perturbations cross the cosmological horizon, either during the contracting phase [237–239] or the expanding Hot Big Bang phase [240, 241], they can lead to the abundant production of SIGWs.

**Early PBH domination.** At distances much larger than the mean PBH separation length, a population of PBHs can be viewed as an effective pressureless fluid. One can then treat this PBH fluid within the context of cosmological perturbation theory showing that the PBH energy fluctuations are isocurvature in nature [242, 243] and can convert to adiabatic curvature perturbations in an early matter-dominated era driven by light PBHs ($m_{PBH} < 10^9$g) occurring before BBN. Interestingly enough, these PBH-induced curvature perturbations can source abundant SIGWs detectable by GW observatories [242–244]. Notably, these PBH associated SIGWs [245, 246] can serve as a novel portal to probe primordial non-Gaussianities (NGs) [247, 248] at small scales ($k > \text{Mpc}^{-1}$) as well the underlying gravity theory [249–251] and Hawking evaporation [252–259].

## 2.4 A simple benchmark scenario: single field USR

In this work, we consider one of the simplest realizations of the single-field scenarios discussed in the previous section. This class of models is described by the inflaton potential in the Einstein frame, $V(\phi)$. The corresponding action can be written as

$$S = \int d^4 x \sqrt{-g} \left( \frac{1}{2} M_P^2 R - \frac{1}{2} (\partial^\mu \phi)^2 - V(\phi) \right), \qquad (2.1)$$

where $R$ is the Ricci scalar and $M_P$ is the reduced Planck mass. Assuming a flat FLRW background geometry $ds^2 = -dt^2 + a^2 dx_i^2$, where $a$ is the scale factor, the background evolution is governed by the Friedmann equation (dots indicate time derivatives)

$$3M_P^2 H^2 = \dot{\phi}^2/2 + V(\phi), \qquad (2.2)$$

with $H = \dot{a}/a$, and the Klein-Gordon equation (a prime denotes a derivative with respect to the field)

$$\ddot{\phi} + 3H\dot{\phi} + V'(\phi) = 0. \qquad (2.3)$$

In order to produce an enhancement of perturbations at LISA scales, and at the same time comply with CMB bounds at large scales, the inflationary potential should feature a shallower
region or an inflection point, which breaks the SR evolution exponentially decelerating the field velocity.

In the class of models considered here, the dynamics can be understood in relatively simple terms. In SR, the inflaton evolves with a negligible acceleration, and the SR solution gives $\dot{\phi} = -V'(3H)$. As the inflaton begins to approach the inflection point, the SR conditions are violated primarily due to a rapid change in the second SR parameter. Having almost reached the local maximum $\phi_*$, the inflaton will spend $O(10)$ e-folds crossing it and its evolution is thus dictated by $\ddot{\phi} + 3H\dot{\phi} + \eta_V(\phi_*)H^2(\phi - \phi_*) \simeq 0$, where $\eta_V \equiv M_P^2V''/V$ denotes the second potential SR parameter. The two solutions of this equation describe two phases: first, a USR-like phase in which the inflaton rapidly decelerates, which leads to an amplification of the power spectrum. Second, a subsequent constant roll or a SR phase that is dual to the initial USR-like phase [82, 260].

The linear superposition of these solutions describes a smooth transition between these epochs. The second SR parameter $\eta_V(\phi_*)$ determines the spectral slope after the peak $n_s - 1 = 3(1 - \sqrt{1 - (4/3)\eta_V(\phi_*)})$. Thus, as exact USR ($\eta_V(\phi_*) = 0$) would produce a scale-invariant spectrum at scales above the spectral peak, the violation of scale invariance in the UV is directly related to the deviation from an exact USR.

As an example, we consider the potential given by the rational function proposed in [86] (see also [89, 261])

$$V(\phi) = \frac{\lambda}{12} \phi^2 (v M_P)^2 \left( 6 - 4b_l \frac{\phi}{v M_P} + 3 \frac{\phi^2}{(v M_P)^2} \right) \left( 1 + b \frac{\phi^2}{(v M_P)^2} \right)^{-2} . \quad (2.4)$$

The presence of an inflection point is enforced by setting

$$b = (1 + b_f) \left[ 1 - \frac{b_l^2}{3} + \frac{b_l^2}{3} \left( \frac{9}{2b_l^2} - 1 \right)^{2/3} \right], \quad (2.5)$$

where we included a tuning parameter $b_f$ allowing for deviation from perfect inflection points (with $b_f > 0$ the inflection point becomes a shallow minimum). The field $\phi$ appearing in the action (2.1) is canonically normalized, and minimally coupled to gravity. This is a proxy for more realistic models in which the inflaton field has a quartic potential and couples non-minimally to gravity via a $\xi R\phi^2$ term, see e.g. [87, 88, 103, 262]. After moving to the Einstein frame, the factor in the denominator appears, which flattens the potential at large field values. In this case, one would further need to canonically normalize the field, and possibly add logarithmic corrections to the coefficients of the monomials in eq. (2.4). We do not discuss the origin of such a potential, as our goal is solely to provide a simple representative model to work with.

We define our benchmark potential by choosing the following parameters (close to the ones used in [261])

$$\begin{align}
\lambda &= 1.4731 \times 10^{-6}, & v &= 0.19688, \tag{2.6} \\
b_l &= 0.71223, & b_f &= 1.87 \times 10^{-5}, \nonumber
\end{align}$$

leading to good agreement with CMB (within $3\sigma$ of current Planck 2018 data [74]) and at the same time to a peak of the curvature spectrum: $P_{\zeta}(k_{\text{peak}}) \simeq 10^{-3}$ at LISA scales. In eq. (2.6)
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Example of a single field model (2.4), with our benchmark parameters (2.6), leading to an early SR phase consistent with CMB data as well as an USR phase which leads to an enhancement of perturbations within the LISA frequency range.</figcaption></figure>

we report 5 significant digits because of the required tuning of the USR potential [261]. In section 3.3 we describe how to compute the spectrum of curvature perturbations in detail. The benchmark potential is depicted in figure 2, where we have arbitrarily normalized the axes using the initial values $\phi_0 = 3M_P$ and $V_0 = 2.3 \cdot 10^{-10} M_P^4$, which are set well before the SR phase that governs the CMB scales.

### 3 Modeling the curvature power spectrum

Throughout this paper, we describe the metric as a small perturbation of the FLRW metric in the longitudinal (conformal Newtonian) gauge

$$ds^2 = g_{\mu\nu}dx^\mu dx^\nu = -a^2(1+2\Phi)d\eta^2 + a^2\left[(1-2\Phi)\delta_{ij} + \frac{1}{2}h_{ij}\right]dx^i dx^j, \quad (3.1)$$

where $\eta$ is the conformal time, and we neglect vector perturbations and anisotropic stress (and so we can identify the two scalar Bardeen potentials, $\Phi = \Psi$).

At linear order in perturbation theory, the time and momentum dependence of the scalar potential $\Phi$ can be factored out by introducing the transfer function $T(\eta, k)$: the Bardeen's potential is related to the gauge-invariant comoving curvature perturbation $\zeta(\mathbf{k})$ by

$$\Phi(\eta, \mathbf{k}) = \frac{3 + 3w}{5 + 3w} T(\eta, k) \zeta(\mathbf{k}) \quad (3.2)$$

where $w$ is the equation of state characterizing the fluid dominating the energy density in the universe. In a radiation-dominated universe, $w = 1/3$, and the prefactor becomes $2/3$. As long as a mode $k$ is super-Hubble ($k\eta < 1$), its evolution is frozen and the transfer function tends to 1. For a homogeneous and isotropic universe, the two-point function of curvature perturbations reads

$$\langle \zeta(\mathbf{k}_1) \zeta(\mathbf{k}_2) \rangle \equiv (2\pi)^3 \delta(\mathbf{k}_1 + \mathbf{k}_2) P_\zeta(k_1) \equiv (2\pi)^3 \delta(\mathbf{k}_1 + \mathbf{k}_2) \frac{2\pi^2}{k_1^3} P_\zeta(k_1), \quad (3.3)$$
where we introduced the dimensionful power spectrum $P_ζ(k)$, and the dimensionless one $P_ζ(k)$. We adopt different methodologies to model the primordial curvature power spectrum $P_ζ(k)$ generated in various inflationary models, by:

* Developing a model-independent parametrization of the spectrum, enabling an agnostic reconstruction approach;

* Constructing analytical templates for $P_ζ(k)$. These templates are motivated by specific scenarios and are consequently less flexible than the previous case. However, given the lower number of parameters controlling the templates as well as their simple analytical descriptions, they alleviate the computational challenges posed by models which would require expensive computation of the curvature power spectrum;

* Establishing a robust pipeline for the computation of primordial curvature power spectra $P_ζ(k)$ for the benchmark model described in section 2.

## 3.1 Model independent parameterization: binned spectrum approach

In this section, we discuss a parameterization of the power spectrum $P_ζ$ in terms of a binned spectrum approach, which is useful in developing a template-independent procedure for computing the SIGW. We represent the curvature power spectrum as a sum of discrete components defined within small momentum intervals. Specifically, we decompose $P_ζ$ as a sum of $N-1$ top hat functions as

$$P_ζ(p) = \sum_{i=1}^{N-1} A_i \Theta(p - p_i) \Theta(p_{i+1} - p). \quad (3.4)$$

The coefficients $A_i$ — which can be approximated as constants — represent the amplitude of the spectrum within the $i$-th bin. The latter is defined within the momentum boundaries $p_i$ and $p_{i+1}$, using the Heaviside step functions $\Theta(p)$.<fnref n="6" /> Any curvature spectrum $P_ζ$ can then be associated with a specific vector $A_i$ of the coefficients appearing in eq. (3.4). Therefore, our computation does not require an a priori choice of a particular template for the momentum dependence of $P_ζ(p)$.

The SIGW depends quadratically on $P_ζ$ using convolution integrals. Hence, we can expect that the ansatz (3.4) allows us to express the SIGW in terms of the double sum of (model-dependent) vector components $A_i$, contracted over a general, model-independent matrix kernel. In section 4.4 we elaborate a simple procedure aimed at performing this task.

## 3.2 Analytical templates for curvature spectra

We divide the list of templates into classes depending on their properties. We can define “smooth shapes” and shapes with features.

<footnote n="6">⁶One could extend this approach to include a tilt of the spectrum at each bin and enforcing continuity for adjacent bins, at the cost of doubling the parameters controlling the power in each bin. We leave this, and other possible extensions for future work.</footnote>
### 3.2.1 Smooth templates

We divide this class of templates into two cases:

**Lognormal (LN).** A typical class of spectral peaks can be characterized by a LN-shape

$$ P_{\zeta}^{\text{LN}}(k) = \frac{A_s}{\sqrt{2\pi}\Delta^2} \exp \left[ -\frac{1}{2\Delta^2} \log^2(k/k_*) \right]. \qquad (3.5) $$

Such spectra appear e.g. in a subset of hybrid inflation and curvaton models, as well as from axion-gauge field coupling, see the discussion in section 2. This template can describe scenarios in which the peak is typically narrow and symmetric in log space. Interestingly, this template allows for an analytic derivation of the GW power spectrum in the broad/narrow peak approximations $\Delta \gg 1$ or $\Delta \ll 1$ (see ref. [263] and improvements in [264]). We will consider the following benchmark scenario choosing

$$ \log_{10} A_s = -2.50, \quad \log_{10} \Delta = \log_{10}(0.5) \approx -0.301, \quad \log_{10} (k_*/\text{s}^{-1}) = -2.00. \qquad (3.6) $$

**Broken power law (BPL).** Another broad class of spectra is encountered, for instance, in single field inflation and curvaton models and can be described by a BPL

$$ P_{\zeta}^{\text{BPL}}(k) = A_s \frac{(\alpha + \beta)^{\gamma}}{(\beta (k/k_{*})^{-\alpha/\gamma} + \alpha (k/k_{*})^{\beta/\gamma})^{\gamma}}, \qquad (3.7) $$

where $\alpha, \beta > 0$ describe respectively the growth and decay of the spectrum around the peak, while $\gamma$ is an $O(1)$ model dependent parameter. The normalization is such that $P_{\zeta}^{\text{BPL}}(k) = A_s$ at the peak. This template provides a very close approximation to the shape of $P_{\zeta}$ obtained from single-field USR scenarios. There, one typically finds $\alpha \simeq 5 - |1 - 2\eta_H| \approx 4$ [81, 82, 265], where $\eta_H$ is the second Hubble SR parameter (see footnote 4) before the USR phase, which is typically close to SR, that is, $\eta_H \ll 1$. For this template, it is possible to find an analytic GW spectrum with $\gamma = 1$, see [266]. We will use the following benchmark with parameters

$$ \log_{10} A_s = -2.71, \quad \log_{10}(k_*/\text{s}^{-1}) = -1.58, \qquad (3.8) \\ \alpha = 3.11, \quad \beta = 0.221, \quad \gamma = 1.25. $$

which provides a curvature power spectrum whose peak is within LISA's sensitivity and that fits the spectrum produced in the ab initio USR scenario described below (see section 3.3).

### 3.2.2 Templates with oscillations

Oscillations in the primordial power spectrum have been extensively studied to seek for deviations from the standard SR inflationary paradigm mainly at scales relevant for CMB or large-scale structure (see e.g. [36, 40, 267, 268]). These oscillations, denoted as primordial features, typically arise due to a sudden transition during inflation [269], occurring over a short time-scale of the order of one e-fold. SIGW provides an opportunity to probe such primordial features at small-scales ($\ll$ Mpc) [180, 182, 270–272]. Specifically, if the mechanism responsible for the enhancement in $P_{\zeta}$ is active when modes are sufficiently deep inside the horizon, the resulting spectrum exhibits oscillations of order one. Due to their large
amplitude, these oscillations could potentially leave their imprints and be detected in the SGWB. This phenomenon can occur in both single-field and multi-field inflationary models. As benchmarks for these large amplitude oscillations, we will consider a small-scale feature induced by a genuine multi-field mechanism and, secondly, weaker oscillations arising from a rapid transition between SR and USR phases in single-field inflation.

**Turns in multi-field inflation.** In multi-field inflationary setups, a common phenomenon is the presence of turns in the field space, being equivalent to the inflationary trajectory deviating from a geodesic in field space. This bending is quantified through the parameter $\eta_{\perp}$, measuring the acceleration of the inflationary trajectory perpendicular to its direction [273, 274] or equivalently the deviation of the trajectory from a geodesic in the target field space. One then can show that sharp and strong (large $\eta_{\perp}$) turns can lead to the following curvature power spectrum, modulated by order-one rapid oscillations<fnref n="7" /> [178–180]

$$ P_{\zeta}^{\text{ST}}(\kappa) = P_{\zeta}^{\text{env}}(\kappa) \left[ 1 + (\kappa - 1) \cos \left( 2e^{-\frac{\delta}{2}} \eta_{\perp} \kappa \right) + \sqrt{(2-\kappa)\kappa} \sin \left( 2e^{-\frac{\delta}{2}} \eta_{\perp} \kappa \right) \right] \Theta(2-\kappa), \quad (3.9) $$

and the envelope

$$ P_{\zeta}^{\text{env}}(\kappa) = A_s \exp(-2\eta_{\perp}\delta) \exp \left[ \frac{2\sqrt{(2-\kappa)\kappa \eta_{\perp}\delta}}{4(2-\kappa)\kappa} \right], \quad (3.10) $$

where $A_s$ denotes the amplitude of the power spectrum in the absence of transient instability and $\kappa \equiv k/k_*$ with $k_*$ being associated with the maximally enhanced scale, deep inside the cosmological Hubble sphere at the time of the sharp turn. The parameter $\delta$ is the duration in e-folds of the turn. $\delta \gtrsim \log \eta_{\perp}$ stands for broad turns and $\delta \lesssim \log \eta_{\perp}$ stands for sharp turns. Finally, $\Theta(x)$ denotes the Heaviside theta function.

We define a parameterization in which the oscillations are switched off via the parameter $F \in [0, 1]$, continuously interpolating between (3.9) when $F = 1$ and (3.10) when $F = 0$

$$ P_{\zeta}(k) = F P_{\zeta}^{\text{ST}}(k) + (1-F) P_{\zeta}^{\text{env}}(k). \quad (3.11) $$

We consider the benchmark scenario whose parameters are given by

$$ \begin{array}{l} \log_{10} A_s = -1.5, \quad \log_{10} (k_*/s^{-1}) = -1.5, \\ \delta = 0.5, \quad \eta_{\perp} = 14, \quad F = 1. \end{array} \quad (3.12) $$

**Rapid transitions between SR and USR phases.** In most USR/inflection-point scenarios, the transition from the initial SR phase to the USR-like phase depends on the properties of the potential. In the presence of sufficiently sharp transitions between these phases, the spectrum of curvature perturbations deviates from a BPL profile due to oscillatory features.

<footnote n="7">⁷Equation (3.9) assumes the entropic field to be massless during the turn. Generalized expressions for other mass choices can be found in [180], where the qualitative features remain analogous.</footnote>
A fully analytic spectrum can be obtained when the initial SR to USR transition is instantaneous. In that case, the peak can be described as (eq. (3.8) in [82])

$$ \mathcal{P}_{\zeta}(k) = \frac{\kappa^2}{4\pi} \left| -\frac{\Gamma(1+\nu_{II})}{\zeta_1} (\kappa/2)^{-\nu_{II}+\frac{1}{2}} J_{\nu_{II}}(\kappa) H_{\nu_I}(\kappa) + \frac{\Gamma(\nu_{II})}{\zeta_2} (\kappa/2)^{-\nu_{II}+\frac{3}{2}} [J_{\nu_{II}}(\kappa) H_{\nu_I-1}(\kappa) + J_{\nu_{II}+1}(\kappa) H_{\nu_I}(\kappa)] \right|^2, \quad (3.13) $$

where, $\kappa = k/k_*$, $H_\nu(x)$ and $J_\nu(x)$ denote the Hankel- and Bessel functions of the first kind, respectively, $k_*$ is the scale of the mode that exits the Hubble sphere during the SR to USR transition, and $\nu_I$ and $\nu_{II}$ are related to the spectral slopes in the attractor phases before and after the transition as $n_s = 2(2 - \nu)$. The first line does not contribute to the peak and can be omitted in case there is no sensitivity to spectral features away from the peak. The parameters $\zeta_1, \zeta_2$ control the amplitude at the IR (or CMB) scales and the peak, respectively.

This spectrum resembles a broken power law with a modulation around the peak. This modulation has a period of $2k_*$ and is damped as $1/k$. It is generated in an instantaneous SR to USR transition and is typically suppressed or removed when the transition is non-instantaneous [82, 275–277]. In this way, such oscillations carry information about the evolution of the inflaton during the transition. Moreover, spectral oscillations can be greatly enhanced in some cases. For instance, when the inflaton rolls through a deep minimum before entering the USR phase [123–128].

We can test for the sensitivity of LISA to resolve these oscillations. To this aim, we only consider the second line of the spectrum (3.13), which describes the peak, and consider a generalized form<fnref n="8" />

$$ \mathcal{P}_{\zeta}^{\text{osc}}(k) = F \mathcal{P}_{\zeta}^{\text{osc,B}}(k) + (1-F) \mathcal{P}_{\zeta}^{\text{osc,BPL}}(k), \quad (3.14) $$

where $F \in [0, 1]$ and with

$$ \mathcal{P}_{\zeta}^{\text{osc,B}}(k) = \pi^2 \kappa^{5-2\nu_{II}} A_s \left| J_{\nu_{II}}(2\kappa) H_{\nu_I-1}^{(1)}(2\kappa) + J_{\nu_{II}+1}(2\kappa) H_{\nu_I}^{(1)}(2\kappa) \right|^2, \quad (3.15) $$

$$ \mathcal{P}_{\zeta}^{\text{osc,BPL}}(k) = A_s \left[ \left( \kappa^{7/2-\nu_I} \frac{(\nu_I+\nu_{II})\Gamma(\nu_I-1)}{\Gamma(\nu_{II}+2)} \right)^{-\gamma} + \left( \kappa^{3/2-\nu_{II}} \right)^{-\gamma} \right]^{-2/\gamma}. \quad (3.16) $$

The rewriting in terms of the two contributions in eqs. (3.15) and (3.16) is equivalent to (3.13) when imposing $F = 1$, but it is done to separate the smooth BPL contribution from the one including the oscillations. This way, changing $F$ smoothly transitions between an oscillating ($F = 1$) and a non-oscillating $F = 0$ power spectra.

The shape of the BPL template is recovered using (3.7) with $\nu_I = (7 - \alpha)/2$ and $\nu_{II} = (3 + \beta)/2$ when $\alpha \le 5$. To obtain an exact match with (3.7) the normalization and the momenta must also be rescaled so that $\mathcal{P}_{\zeta}^{\text{BPL}}(k) = n \mathcal{P}_{\zeta}^{\text{osc,BPL}}(bk)$, where

$$ n \equiv b^{\beta}(1 + \beta/\alpha)^{\gamma}, \quad b \equiv \left[ \frac{4(\alpha/\beta)^{\gamma} \Gamma((\beta + 7)/2)^2}{(-\alpha + \beta + 10)^2 \Gamma((5 - \alpha)/2)^2} \right]^{\frac{1}{\alpha + \beta}} \quad (3.17) $$

<footnote n="8">⁸For notational simplicity, the template uses a different normalization and scaling of the argument than (3.13).</footnote>
or, equivalently, as $A_s \to nA_s$ and $k_* \to k_*/b$. The benchmark scenario corresponds to

$$ \begin{gathered} \log_{10} A_s = -2.58, \quad \log_{10}(k_*/\text{s}^{-1}) = -2.02 \\ \nu_I = 1.95, \quad \nu_{II} = 1.61, \quad \gamma = 1.67, \quad F = 1. \end{gathered} \qquad (3.18) $$

We fix the parameters to match the BPL template in the limit $F = 0$.

### 3.3 Computation of $P_ζ$ in single field USR scenarios

The benchmark model we introduced in section 2.4 is based on a single-field model of inflation featuring a phase of USR. We now describe in detail how to compute the spectrum of curvature perturbations using linear perturbation theory.

As a warm-up, let us define the system of equations in terms of dimensionless variables rescaled to the corresponding relevant quantities. This typically improves the numerical stability of a code computing the spectrum of curvature fluctuations. We define

$$ x \equiv \phi/M_P, \quad U(x) \equiv V(\phi)/V_0, \qquad (3.19) $$

where we introduce the suffix 0 to indicate the quantities evaluated at the initial conditions of the background evolution. We can also define the dimensionless time and Hubble rate using

$$ \tau \equiv tV_0^{1/2}/M_P, \quad h \equiv HM_P/V_0^{1/2}. \qquad (3.20) $$

We can now change the evolution variable from time to the number of e-folds $N$, defined as $dN \equiv Hdt = hd\tau$, setting $N = 0$ at $\tau_0$. In this way, the background equations of motion (2.2) and (2.3) become

$$ y' + 3 \left[ 1 - \frac{y^2}{6} \right] y + \frac{U_{,x}}{h^2} = 0, \qquad (3.21a) $$

$$ x' - y = 0, \qquad (3.21b) $$

$$ h' + \frac{(x')^2}{2}h = 0, \qquad (3.21c) $$

where prime denotes derivatives with respect to $N$, while $U_{,x} \equiv dU(x)/dx$. The initial conditions for the inflaton and the Hubble parameters can be found assuming that initially, SR is satisfied, which for a given initial $x_0$ leads to

$$ y_0 = -\frac{3}{U_{0,x}} \left( \sqrt{1 + \frac{2}{3} U_{0,x}^2} - 1 \right), \quad h_0 \equiv \frac{1}{\sqrt{6}} \left( \sqrt{1 + \frac{2}{3} U_{0,x}^2} - 1 \right)^{1/2} . \qquad (3.22) $$

Finally, in the code, we keep track of the equation of state during inflation $w_{\text{inf}} \equiv p_{\phi}/\rho_{\phi}$, which can be written in our notation as $w = [(x')^2 - 3]/3$. As inflation stops when the equation of state becomes larger than $w > -1/3$, we stop the evolution when $(x')^2 = 2$.

The resulting inflationary background can be described by the evolution of the Hubble rate $H$. This is dictated by dynamical equations relating $H$ to the SR parameters, which are defined as

$$ \epsilon_H \equiv -\frac{\dot{H}}{H^2} = \frac{1}{2M_P^2} \left( \frac{d\phi}{dN} \right)^2 \equiv \frac{y^2}{2}, \qquad (3.23) $$

$$ \eta_H \equiv -\frac{\dot{H}}{2H\dot{H}} = \epsilon_H - \frac{1}{2} \frac{d \log \epsilon_H}{dN} = \frac{y^2}{2} - \frac{y'}{y}, \qquad (3.24) $$
where we introduced $y = x'$. Notice that our definition of $\eta_H$ differs from another definition often used in the literature, which is expressed as $\dot{\epsilon}_H/(\epsilon_H H)$. During the USR phase, this definition equals approximately $-2\eta_H$ neglecting $\epsilon_H$ corrections.

As long as the SR approximation is valid, i.e. $\epsilon_H \ll 1$ and $\eta_H \ll 1$, the power spectrum of curvature perturbations (see e.g. [278]) is given by

$$ \mathcal{P}_{\zeta}(k) = \frac{H_k^2}{8\pi^2 M_P^2 \epsilon_{H,k}} = \frac{V_0}{M_P^4} \left( \frac{h_k}{2\pi y_k} \right)^2, \qquad (3.25) $$

where the suffix $k$ indicates that these quantities are evaluated at Hubble crossing $k = aH$. Consequently, one can also show that the scalar spectral index $n_s - 1 \equiv \text{d} \ln \mathcal{P}_{\zeta}/\text{d} \ln k$ and the tensor-to-scalar ratio $r \equiv \mathcal{P}_h/\mathcal{P}_{\zeta}$ are given by the well-known expressions

$$ n_s \approx 1 - 4\epsilon_H + 2\eta_H, \qquad (3.26) $$

$$ r \approx 16\epsilon_H, \qquad (3.27) $$

at leading order in the SR parameters.

The approximations above can not be used in the context of USR scenarios, as the large deceleration of the inflaton velocity causes $\eta_H \sim \mathcal{O}(1)$ ($\eta_H = 3$ in the limit of exponential deceleration). We go beyond the SR approximation by solving the Mukhanov-Sasaki (MS) equation [279, 280]. Introducing momentarily the conformal time $\eta$ such that $d\eta \equiv dt/a$, we can define the MS variable $v \equiv a\delta\phi$ in terms of the inflaton perturbations in the spatially flat gauge, which satisfies the EoM

$$ \frac{\partial^2 v}{\partial \eta^2} + \left( k^2 - \frac{1}{z} \frac{\partial^2 z}{\partial \eta^2} \right) v = 0, \qquad (3.28) $$

where $z \equiv a^2(\partial\phi/\partial\eta)/(\partial a/\partial\eta)$. We can relate the curvature perturbation to $v$ at linear order in perturbation theory using $\zeta = v/z = H\delta\phi/(\partial\phi/\partial t) = \delta x/y$. Assuming initial adiabatic vacuum, each mode $k$ is fixed in the asymptotic past at $\eta_{\text{in}} \ll -1/k$ as

$$ v(\eta) = e^{-ik(\eta-\eta_{\text{in}})} / \sqrt{2k}, \qquad (3.29) $$

which translates into boundary conditions for eq. (3.28) of the form

$$ v_{\text{in}} = \frac{1}{\sqrt{2k}}, \qquad \frac{\partial v_{\text{in}}}{\partial \eta} = -i \sqrt{\frac{k}{2}}. \qquad (3.30) $$

One can improve the stability of the numerical computation by defining a rescaled variable $\tilde{\delta\phi} \equiv a_{\text{in}}e^{ik(\eta-\eta_{\text{in}})}\delta\phi\sqrt{2k}$ and the corresponding dimensionless quantity $\tilde{\delta x} \equiv \tilde{\delta\phi}/M_P$. The prefactor absorbs the sub-Hubble oscillations, simplifying the time evolution of the sub-Hubble phase. Next, we introduce the dimensionless momentum $\kappa \equiv kM_P/V_0^{1/2}$. Finally, moving to the number of e-folds as the time variable, the MS equation becomes

$$ \frac{\mathrm{d}^2\tilde{\delta x}{\mathrm{d}N^2} + \left(3 - \frac{1}{2}y^2 - \frac{2i\kappa}{e^N h}\right) \frac{\mathrm{d}\tilde{\delta x}{\mathrm{d}N} + \left(\frac{U_{,xx} + 2U_{,xy}}{h^2} + 3y^2 - \frac{1}{2}y^4 - \frac{2i\kappa}{e^N h}\right) \tilde{\delta x}{\mathrm{d}N}{\mathrm{d}N} = 0, \qquad (3.31) $$
<figure><img src="image_5.png" /><figcaption>Figure 3. From top to bottom: curvature power spectrum $\mathcal{P}_{\zeta}$; Hubble SR parameters $\epsilon_H$ and $\eta_H$; dimensionless inflaton field $x$; dimensionless inflaton field velocity $y$; rescaled Hubble parameter $h$. In the second panel, $\eta_H$ is initially small and negative and it transitions to positive values close to $\eta_H \sim 3$ during the USR phase. All quantities are plotted as a function of the number of e-folds to the end of inflation. The vertical dashed lines indicate the epoch when the CMB pivot scale $k_{\text{ref}} = 0.05/\text{Mpc}$ crossed the Hubble sphere.</figcaption></figure>

to be solved with initial conditions $\tilde{\delta\phi}_{\text{in}} = 1$, and $\tilde{\delta\phi}'_{\text{in}} = -1$, while the curvature power spectrum is extracted using

$$ \mathcal{P}_{\zeta}(k) = \frac{1}{4\pi^2} \frac{V_0}{M_P^4} \left( \frac{\kappa}{e^{N_{\text{in}}}} \right)^2 \frac{|\tilde{\delta x}|^2}{y^2}. \qquad (3.32) $$

In figure 3 we show the background evolution of $x(N)$, $y(N)$ and $h(N)$ obtained in the benchmark scenario (2.6). For convenience, we define $N_{\text{end}}$ as the number of e-folds at the end of inflation. In the same plot (top panels), we report the evolution of the Hubble SR parameters and the curvature power spectrum. For simplicity, we identify the momentum with the corresponding time of Hubble crossing $k = a(N)H(N)$. On top, we also indicate the momentum in units of $1/s$. The CMB reference scale crosses the Hubble sphere $N - N_{\text{end}} = -58$ e-folds before the end of inflation, which is indicated with a dashed
vertical line. In the top panel, the curvature power spectrum is shown both using the SR approximation (3.25) (blue line), which however fails to reproduce the spectrum around the peak where the USR phase takes place. We show with an orange line the full spectrum computed solving eq. (3.31). At $k = 0.05$/Mpc where the SR expressions (3.25)–(3.27) are valid, we find $P_ζ = 2.12 \times 10^{-9}$, $n_s = 0.952$, and $r = 0.00726$. These numbers are compatible with the latest observational bounds [74, 281]

$$P_ζ = (2.10 \pm 0.03) \times 10^{-9}, \quad n_s = 0.9649 \pm 0.0042, \quad r < 0.036 \quad (3.33)$$

at the scale $k_{\text{ref}} \equiv 0.05$/Mpc, reported here at 68% CL for $A_s$ and $n_s$, and at 95% for $r$. The low value of $n_s$ is a rather common feature of models featuring an USR phase not sufficiently far from the region of the potential controlling the CMB scales [87, 88, 103, 262], as it is the case if one considers enhancements in the LISA band. This benchmark USR scenario leads to a power spectrum within the LISA band which can be fitted with a BPL template with parameters (3.8).

Let us mention here that it was recently suggested that USR dynamics, in the extremal case of large spectral enhancement with $P_ζ \sim O(10^{-2})$ leading to PBH formation, may violate perturbativity and induce loop corrections that could also affect much longer modes associated with the scales observed through the CMB [153]. While the existence and magnitude of this effect for soft modes is still under debate [133, 154–160, 162–168, 282–291], recent analyses show that for realistic transitions into and out of USR perturbativity is retained (see e.g. for analytical [156, 159, 163] and lattice results [292]). In this work, we restrict to adopting linear perturbation theory to compute the $P_ζ$ at the LISA scales and leave further refinements including loop corrections to the spectrum of curvature perturbations for future work.

**Non-Gaussianities.** As a final note, let us comment on the NGs expected in this benchmark scenario. In this class of USR models, the peak in $P_ζ$ is generated during a brief phase of USR inflation, which transitions into constant-roll (i.e. constant $\eta_H$) inflation afterwards [82, 293, 294]. In these scenarios, NGs are controlled by the details of the transition between USR and the subsequent phase, which is related to the UV tilt of the spectrum [97, 295–298]. We can expect the non-linear curvature perturbation to take the form

$$\zeta = -\frac{2}{\beta} \log \left( 1 - \frac{\beta}{2} \zeta_G \right) = \zeta_G + \frac{3}{5} f_{\text{NL}} \zeta_G^2 + \dots, \quad (3.34)$$

with $f_{\text{NL}} = 5\beta/12 \simeq 0.092$ in our benchmark scenario (3.8). Note, however, that this only takes into account the local generation of non-Gaussianity on super-Hubble scales by considering Gaussian inflaton fluctuations, while it has recently been shown using simulations that nonlinear interactions can also generate a large amount of non-Gaussianity intrinsic to the inflaton [292]. Finally, let us mention that such small NG for perturbations of the typical amplitude considered here would lead to negligible contributions to the SIGW spectrum (see more discussion in section 4.5), even beyond standard perturbation theory [69].
## 4 Computation of the scalar-induced GW background

Primordial scalar perturbations are frozen during their super-Hubble evolution; only when re-entering the Hubble radius after inflation ends do they start evolving in time. Moreover, even if scalar, vector and tensor modes are independent at the first order in perturbation theory, they couple when going to higher orders in fluctuations. For example, scalar modes source tensor modes and thus produce GWs.<fnref n="9" />

In the metric defined in eq. (3.1), we ignore tensor perturbations generated at first order $h_{ij}^{(1)}$ and consider scalar perturbations that act as a source of second order tensor modes at $h_{ij}^{(2)}$. The exact evolution of $h_{ij}^{(2)}$ can be obtained from the spatial part of the Einstein equations after applying the projection tensor $\mathcal{T}_{ij}^{lm}$, which selects the transverse-traceless component. In the absence of anisotropic stress, at second order (note that we drop the superscript indicating the order in perturbation theory from now on) one obtains [15, 16, 18, 19]

$$h''_{ij}(\eta, \mathbf{x}) + 2\mathcal{H}h'_{ij}(\eta, \mathbf{x}) - \nabla^2 h_{ij}(\eta, \mathbf{x}) = -4\mathcal{T}_{ij}^{lm}\mathcal{S}_{lm}(\eta, \mathbf{x}), \quad (4.1)$$

where $'$ is the derivative with respect to conformal time $\eta$, $\mathcal{H} = a'/a$ denotes the conformal Hubble parameter, and $\mathcal{S}_{ij}$ is the source term

$$\mathcal{S}_{ij}(\eta, \mathbf{x}) = 4\Phi\partial_i\partial_j\Phi + 2\partial_i\Phi\partial_j\Phi - \frac{4}{3(1+w)}\partial_i\left(\frac{\Phi'}{\mathcal{H}} + \Phi\right)\partial_j\left(\frac{\Phi'}{\mathcal{H}} + \Phi\right). \quad (4.2)$$

In the last equation $w$ is the equation-of-state parameter of the Universe, and the scalar perturbation $\Phi(\eta, \mathbf{x})$ can be related to the gauge-invariant comoving curvature perturbation $\zeta$. Solving eq. (4.1) in Fourier space, as we review in section 4.1, we obtain the (time-averaged) dimensionless power spectrum $\mathcal{P}_h(\eta_f, k)$ of GWs at a time $\eta_f$ after the end of their production. The fractional energy density of GWs per logarithmic interval in frequency is given by

$$\Omega_{\text{GW}}(\eta_f, k) \equiv \frac{\rho_{\text{GW}}(\eta_f, k)}{\rho_c(\eta_f)} = \frac{1}{24} \left(\frac{k}{\mathcal{H}(\eta_f)}\right)^2 \overline{\mathcal{P}_h(\eta_f, k)}, \quad (4.3)$$

with

$$\rho_{\text{GW}}(\eta) = \int d\ln k \rho_{\text{GW}}(\eta, k), \quad (4.4)$$

the bar denoting an oscillation average and $\rho_c$ is the critical energy density. GWs redshift as relativistic species and the current abundance of the GWB can be obtained by accounting for entropy injections in the standard thermal history of the Universe:

$$\Omega_{\text{GW}}(k)h^2 = \Omega_{r,0}h^2 \left(\frac{g_*(\eta_f)}{g_*^0}\right) \left(\frac{g_{*s}^0}{g_{*s}(\eta_f)}\right)^{4/3} \Omega_{\text{GW}}(\eta_f, k). \quad (4.5)$$

Here $g_*(g_{*,s})$ are the relativistic degrees of freedom in energy (entropy), $\Omega_{r,0}h^2 = 4.2 \cdot 10^{-5}$ is the current radiation density if the neutrino were massless [303]. Assuming standard model degrees of freedom, one finds that [304]

$$c_g(\eta_f) \equiv \left(\frac{g_*(\eta_f)}{g_*^0}\right) \left(\frac{g_{*s}^0}{g_{*s}(\eta_f)}\right)^{4/3} \approx 0.39, \quad (4.6)$$

for $\eta_f$ of relevance for LISA.

<footnote n="9">⁹The production of scalar modes from primordial tensor modes is discussed in [299, 300] Moreover, very recently, SIGWs sourced by scalar-tensor perturbations have also been analyzed [301, 302].</footnote>
## 4.1 Source of GWs at second order in the scalar perturbations

The solution to eq. (4.1) can be easily found in Fourier space, where the equation for the GW amplitude $h$ for each polarization state $s$ becomes

$$h_s''(\mathbf{k}, \eta) + 2\mathcal{H}h_s'(\mathbf{k}, \eta) + k^2h_s(\mathbf{k}, \eta) = 4S_s(\mathbf{k}, \eta) \quad (4.7)$$

and where $S_s(\mathbf{k}, \eta)$ encloses the Fourier transform of the (projected) source given by

$$S_s(\mathbf{k}, \eta) = \int \frac{d^3\mathbf{p}}{(2\pi)^3} Q_s(\mathbf{k}, \mathbf{p}) f(|\mathbf{k} - \mathbf{p}|, p, \eta) \zeta(\mathbf{p}) \zeta(\mathbf{k} - \mathbf{p}). \quad (4.8)$$

In the latter equation, we introduced

$$f(|\mathbf{k} - \mathbf{p}|, p, \eta) = \frac{3(1+w)}{(5+3w)^2} \left[ 2(5+3w)T(|\mathbf{k} - \mathbf{p}|\eta)T(p\eta) + \frac{4}{\mathcal{H}^2}T'(|\mathbf{k} - \mathbf{p}|\eta)T'(p\eta) + \frac{4}{\mathcal{H}}(T(|\mathbf{k} - \mathbf{p}|\eta)T'(p\eta) + T'(|\mathbf{k} - \mathbf{p}|\eta)T(p\eta)) \right] \quad (4.9)$$

and we introduced the curvature perturbation transfer function $T(k\eta)$, and the spherical coordinates $(p, \theta, \phi)$ for the internal momentum $\mathbf{p}$ and aligned the axes $(\hat{x}, \hat{y}, \hat{z})$ with $(e_+(\mathbf{k}), e_\times(\mathbf{k}), \mathbf{k})$ (with $e_+$, $e_\times$ being the polarisation tensors for the GW) so that

$$Q_s(\mathbf{k}, \mathbf{p}) = e_s^{ij}(\hat{\mathbf{k}}) p_i p_j = \frac{p^2}{\sqrt{2}} \sin^2 \theta \times \begin{cases} \cos 2\phi, & s = + \\ \sin 2\phi, & s = - \end{cases} \quad (4.10)$$

A solution to the Fourier transform of the inhomogeneous equation of motion for $h_{ij}(\eta, k)$, eq. (4.1), can be obtained using the Green's function $G_k(\eta, \bar{\eta})$, that solves the homogeneous equation

$$G_k''(\eta, \bar{\eta}) + \left(k^2 - \frac{a''}{a}\right) G_k(\eta, \bar{\eta}) = \delta(\eta - \bar{\eta}), \quad (4.11)$$

with the boundary conditions $\lim_{\eta \to \bar{\eta}} G_k(\eta, \bar{\eta}) = 0$ and $\lim_{\eta \to \bar{\eta}} G'_k(\eta, \bar{\eta}) = 1$. The Green's function depends on $k = |\mathbf{k}|$ by isotropy and can be expressed analytically in terms of Bessel functions as

$$k G_k(\eta, \bar{\eta}) = \sqrt{x\bar{x}} [y_\nu(x)j_\nu(\bar{x}) - j_\nu(x)y_\nu(\bar{x})], \quad (4.12)$$

where $x = k\eta$, $\bar{x} = k\bar{\eta}$, and $\nu = 3(1-w)/[2(1+3w)]$ and $j_\nu$ and $y_\nu$ are respectively the spherical Bessel function of the first and second kind. For example, during RD, $kG_k(\eta, \bar{\eta}) = \sin(x - \bar{x})\Theta(\eta - \bar{\eta})$, where $\Theta$ is the Heaviside function.

The amplitude of the tensor modes can then be written as

$$\begin{aligned} h_s(\mathbf{k}, \eta) &= 4 \int_{\eta_i}^{\eta} d\bar{\eta} G_k(\eta, \bar{\eta}) \frac{a(\bar{\eta})}{a(\eta)} S_s(\mathbf{k}, \bar{\eta}) \\ &= 4 \int_{\eta_i}^{\eta} d\bar{\eta} G_k(\eta, \bar{\eta}) \frac{a(\bar{\eta})}{a(\eta)} \int \frac{d^3\mathbf{p}}{(2\pi)^3} Q_s(\mathbf{k}, \mathbf{p}) f(|\mathbf{k} - \mathbf{p}|, p, \bar{\eta}) \zeta(\mathbf{p}) \zeta(\mathbf{k} - \mathbf{p}), \end{aligned} \quad (4.13)$$
with s indicating the polarisation and $\eta_i$ the initial emission time. The GW two-point function, needed to obtain the energy density of GWs as a function of the scalar power spectrum, is defined in terms of the GW power spectrum

$$ \langle h^r(\eta, \mathbf{k}_1) h^s(\eta, \mathbf{k}_2) \rangle = (2\pi)^3 \delta(\mathbf{k}_1 + \mathbf{k}_2) \delta^{rs} \frac{2\pi^2}{k_1^3} \mathcal{P}_h(k_1). \quad (4.14) $$

After substituting eq. (4.12) into eq. (4.13), the final expression for the second-order induced tensor power spectrum reads [305, 306]

$$ \langle h^{s1}(\eta, \mathbf{k}_1) h^{s2}(\eta, \mathbf{k}_2) \rangle = 16 \int \frac{d^3 \mathbf{p}_1}{(2\pi)^3} \frac{d^3 \mathbf{p}_2}{(2\pi)^3} Q_{s1}(\mathbf{k}_1, \mathbf{p}_1) Q_{s2}(\mathbf{k}_2, \mathbf{p}_2) \\ \times I(|\mathbf{k}_1 - \mathbf{p}_1|, p_1, \eta_1) I(|\mathbf{k}_2 - \mathbf{p}_2|, p_2, \eta_2) \langle \zeta_{\mathbf{p}_1} \zeta_{\mathbf{k}_1 - \mathbf{p}_1} \zeta_{\mathbf{p}_2} \zeta_{\mathbf{k}_2 - \mathbf{p}_2} \rangle. \quad (4.15) $$

In the Gaussian case, where only the disconnected part of the trispectrum survives (hence the “d” in the following equation, we discuss in section 4.5 the impact of primordial NG), one obtains

$$ \overline{\mathcal{P}_{h,d}(\eta, k)} = 4 \int_{0}^{\infty} dv \int_{|1-v|}^{1+v} du \left[ \frac{4v^2 - (1+v^2-u^2)^2}{4uv} \right]^2 \overline{I^2(u, v, k, \eta)} \mathcal{P}_{\zeta}(kv) \mathcal{P}_{\zeta}(ku), \quad (4.16) $$

where we introduced the dimensionless variables

$$ v \equiv \frac{p}{k}, \qquad u \equiv \frac{|\mathbf{k} - \mathbf{p}|}{k}. \quad (4.17) $$

The overline stands for an oscillation average, and the kernel function $I(u, v, \eta)$ is defined in terms of Green's function as

$$ I(|\mathbf{k} - \mathbf{p}|, p, \eta) = \int_{\eta_i}^{\eta} d\bar{\eta} G_k(\eta, \bar{\eta}) \frac{a(\bar{\eta})}{a(\eta)} f(|\mathbf{k} - \mathbf{p}|, p, \bar{\eta}). \quad (4.18) $$

Since the integration domain of (4.16) is rectangular, for computational purposes, it is convenient to rotate the coordinates into

$$ t \equiv u + v - 1, \qquad s \equiv u - v \quad (4.19) $$

and the SIGW spectrum $\overline{\mathcal{P}_h(\eta, k)}$ can be then rewritten as

$$ \overline{\mathcal{P}_h(\eta, k)} = 4 \int_0^\infty dt \int_0^1 ds \left[ \frac{t(2+t)(1-s^2)}{(1-s+t)(1+s+t)} \right]^2 \overline{I^2(t, s, k, \eta)} \\ \times \mathcal{P}_\zeta \left( k \frac{t+s+1}{2} \right) \mathcal{P}_\zeta \left( k \frac{t-s+1}{2} \right), \quad (4.20) $$

where the integration in s is restricted to positive values due to the integrand being an even function of s. The integration kernel $I(t, s, \eta)$ contains information about the time evolution of the source during emission, as well as the propagation of the emitted GWs after emission, and thus depends on the thermal history when the relevant modes re-entered the Hubble radius. Let us consider different assumptions on the thermal history in the following sections.
The generation of tensor modes at second-order in perturbation theory raises concerns about the potential gauge dependence of results commonly calculated in the Newtonian gauge. Unlike first-order tensor modes, which are gauge invariant, second-order tensor modes are gauge dependent [307]. During the phase when the source is still active, the result is expected to remain gauge-dependent, as one cannot directly identify the tensor mode with the freely propagating GW. However, when the source becomes inactive after the GWs are produced, it effectively decouples from the GWs. This happens for example during the radiation-dominated era of the Universe and in the other cases considered in this draft. Therefore, in the late-time limit well inside the cosmological horizon, tensor mode behaves as linear metric perturbations and the initial gauge dependence no longer affects the final result [308–311], ensuring that the spectra computed in this work are unaffected by this issue.

## 4.2 Radiation-dominated era

If the emission takes place in a RD universe, the kernel function in the deep sub-horizon regime $k\eta \to \infty$ takes the form [51, 305, 306]

$$ \overline{I_{\text{RD}}^2(t, s, k\eta \to \infty)} = \frac{1}{2(k\eta)^2} I_A^2(u, v) [I_B^2(u, v) + I_C^2(u, v)] \quad (4.21) $$

where

$$ \begin{align} I_A(u, v) &= \frac{3(u^2 + v^2 - 3)}{4u^3v^3}, \nonumber \\ I_B(u, v) &= -4uv + (u^2 + v^2 - 3) \log \left| \frac{3 - (u+v)^2}{3 - (u-v)^2} \right|, \tag{4.22} \\ I_C(u, v) &= \pi(u^2 + v^2 - 3)\Theta(u + v - \sqrt{3}), \nonumber \end{align} $$

$u$ and $v$ have been introduced above and again $\Theta(x)$ is the Heaviside function. Notice that the unphysical divergence obtained in the limit $|k - p| + p = \sqrt{3}k$, which is also retained in the spectrum produced by monochromatic scalar perturbations [18]. The factor of $\sqrt{3}$ originates from the (inverse) sound speed in RD appearing in the transfer function, and in the limit $|k - p| + p = \sqrt{3}k$ the contributions from some of the source terms add up constructively and build up logarithmically over time [305]. Assuming RD up to $\eta \to \infty$ leads to this logarithmic divergence, which is regularized in the integral for $P_h$ if $P_\zeta(k)$ is smooth enough (e.g. has a nonzero width), or if the emission time is finite. We do not introduce this regulator here, as it is numerically irrelevant for spectra with finite width [6].

## 4.3 Transition from an early matter-dominated to the radiation-dominated era

We further consider the alternative thermal history wherein an early matter-dominated (eMD) era may precede the RD era. We follow the prescription of [312], assuming a sudden transition at a conformal time $\eta = \eta_R$, where the subscript $R$ indicates the reheating time [306]. In this scenario, the dominant contribution comes from GWs induced during the RD era by the scalar perturbations that have entered the horizon during an eMD era.

Interestingly, in this case, one observes a resonantly enhanced production of GWs. In particular, during the transition, the time derivative of the Bardeen potential $\Phi$, which is the source of the SIGW signal, goes very quickly from $\Phi' = 0$ (since in an eMD era $\Phi =$
constant) to $\Phi' \neq 0$ (see [243, 309] for more details), leading to an enhanced secondary tensor mode production sourced mainly by the $H^2\Phi'^2$ term in eq. (4.2). In a more physical scenario, the transition happens more gradually [313–315].

For the kernel function $I(t, s, k, \eta)$, one finds two dominant contributions at the onset of the late RD era, i.e. at $\eta = \eta_R$, when most of the GWs are expected to be produced [309, 316]. The first contribution to $I(t, s, k, \eta)$ is given by $k_{\max}/k \sim 1$, at $t_0 = \sqrt{3} - 1$ [312]

$$ \overline{I_{\text{IRD, res}}^2}(t_0, s, k, \eta_R) \approx Y \frac{9(-5+s^2+2t+t^2)^4 x_R^8}{2^{17}5^4(1-s+t)^2(1+s+t)^2} \text{Ci}^2(y), \quad (4.23) $$

where $x_R = k\eta_R$. The variable $y$ is defined as $y \equiv \frac{|t+1-c_s^{-1}|x_R}{2c_s^{-1}} = \frac{|t+1-\sqrt{3}|x_R}{2\sqrt{3}}$, and $Y$ is a fudge factor to absorb the uncertainty in the integration boundary, set here to be 2.3 as in [312]. At $t_0 = \sqrt{3} - 1$, the logarithmic singularity of the function Ci is reached, giving rise to the peak in the spectrum.

The second contribution to $I(t, s, k, \eta)$ comes from the wave-numbers satisfying $k_{\max}/k \gg 1$, hence the integrations are dominated by the large $t$ region $u \sim v \sim t \gg 1$. Therefore, the dependence on $s$ is lost. Setting $s = 0$, the kernel function reads

$$ \overline{I_{\text{IRD,LV}}^2}(t \gg 1, s, k, \eta_R) \approx \frac{9t^4 x_R^8}{2^{17}5^4} \left[ 4\text{Ci}^2(x_R/2) + (\pi - 2\text{Si}(x_R/2))^2 \right]. \quad (4.24) $$

The integration region is

$$ 0 \le s \le 1 \quad \text{and} \quad 0 \le t \le -s + 2 \frac{k_{\max}}{k} - 1 \quad \text{for} \quad k \le k_{\max}, \quad (4.25) $$

and the result obtained from this integration region is then doubled to account for $s \to -s$. The two contributions are computed separately and added to give

$$ \Omega_{\text{GW}}^{\text{eMDRD}} \approx \Omega_{\text{GW}}^{(\text{LV})} + \Omega_{\text{GW}}^{(\text{res})}. \quad (4.26) $$

## 4.4 Computation of SIGW with the binned spectrum approach

In this section, we discuss in more detail the procedure sketched in section 3.1 for computing the SIGW in terms of a template-free approach to the curvature power spectrum. We express the power spectrum as a sum over momentum bins, as in eq. (3.4). The specific profile for the power spectrum $P_\zeta$ is then associated with a vector of coefficients $A_i$. By plugging eq. (3.4) into eqs. (4.3) and (4.20), we recast the SIGW density into the sum

$$ \Omega_{\text{GW}}(k) = \sum_{i,j}^{N-1} \Omega_{\text{GW}}^{(i,j)}(k) A_i A_j \quad (4.27) $$

performed over the momentum bins. The kernel for this sum is the matrix

$$ \Omega_{\text{GW}}^{(i,j)}(k) = \frac{1}{12} \left(\frac{k}{aH}\right)^2 \int_0^\infty dt \int_0^1 ds \left[ \frac{t(2+t)(s^2-1)}{(1-s+t)(1+s+t)} \right]^2 \overline{I^2(t,s,k,\eta)} \\ \times \Theta(kv(s,t)-p_i) \Theta(p_{i+1}-kv(s,t)) \Theta(ku(s,t)-p_j) \Theta(p_{j+1}-ku(s,t)), \quad (4.28) $$

where $p_i, p_j$ are the boundaries of the momentum bins entering in eq. (3.4).
Importantly, we stress that the matrix $\Omega_{GW}^{(i,j)}$ of eq. (4.27) is *independent* of the specific scalar spectrum considered, and the information on $\mathcal{P}_{\zeta}$ is stored only in the coefficients $A_i$ appearing quadratically in eq. (4.27). This implies that the computation of eq. (4.28) depends only on the kernel function, and its entries can be computed once for all for any given cosmology: for example, we can use one of the kernels discussed in section 4.2 or section 4.3. The nested integrals appearing in eq. (4.28) should then be performed a single time for each kernel. Once $\Omega_{GW}^{(i,j)}$ is determined, it can be used to swiftly compute the resulting $\Omega_{GW}$ for any $\mathcal{P}_{\zeta}$, by means of the contractions in eq. (4.27). In fact, with this method, we reduce the problem of computing $\Omega_{GW}$ to perform the simple sum of eq. (4.27).

This approach is useful in scenarios where the underlying shape of the curvature spectrum is not accurately known, for example, due to the presence of peaks or breaks, whose position depends on the underlying physics we wish to probe. In fact, the method allows us to scan over different sets of $A_i$ components, swiftly computing the SIGW frequency profile, which can then be compared with data. Other approaches for reconstructing the properties of the underlying $\mathcal{P}_{\zeta}$ from SIGW data can be found for example in [317–319].

We tabulate the matrix (4.27) assuming a varying number of bins $N$ in the range of relevance for LISA, which is $k \in [1.26 \times 10^{-4}, 6.28]/\text{s}$ both for the internal $(i, j)$ indices, as well as external momentum $k$. This range is chosen to match the one used by the SGWBinner code adopted to perform the LISA forecasts. See the discussion in section 5.

## 4.5 Non-Gaussian imprints on the SIGW spectrum

From the solution of the SIGW, eq. (4.13), one can relate the tensor power spectrum to the four-point correlation function of the curvature perturbation, see eq. (4.15). As anticipated above, the latter can be decomposed into disconnected and connected contributions, where the connected part vanishes when primordial fluctuations are drawn from a Gaussian distribution. The disconnected contribution gives rise to eq. (4.16) that can be solely expressed in terms of the scalar power spectrum $\mathcal{P}_{\zeta}(k)$. However, the connected contribution depends on the primordial trispectrum $\langle \zeta_{k_1} \zeta_{k_2} \zeta_{k_3} \zeta_{k_4} \rangle'_c = T_{\zeta}(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3, \mathbf{k}_4)$, whose corresponding tensor power spectrum reads [51]$^{10}$

$$ \overline{\mathcal{P}_{h,c}} = \frac{1}{4\pi} \int_0^\infty dv_1 \int_{|1-v_1|}^{1+v_1} du_1 \int_0^\infty dv_2 \int_{|1-v_2|}^{1+v_2} du_2 \int_0^{2\pi} d\psi \\ \times \frac{\cos(2\psi)}{(u_1 v_1 u_2 v_2)^{5/4}} [4v_1^2 - (1+v_1^2-u_1^2)^2] [4v_2^2 - (1+v_2^2-u_2^2)^2] \\ \times I(u_1, v_1, k, \eta) I(u_2, v_2, k, \eta) \mathcal{T}_{\zeta}(u_1, v_1, u_2, v_2, \psi), \quad (4.29) $$

and we use the variables $u_i$ and $v_i$ defined above. The dimensionless trispectrum function $\mathcal{T}_{\zeta}$ is defined as

$$ \mathcal{T}_{\zeta}(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3, \mathbf{k}_4) = \frac{(k_1 k_2 k_3 k_4)^{9/4}}{(2\pi)^6} T_{\zeta}(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3, \mathbf{k}_4), \quad (4.30) $$

<footnote>$^{10}$In [51], which uses different conventions, their eq. (2.6) should be divided by 4, as well as subsequent results. This has been corrected in [52] in a study of SIGWs including parity violation.</footnote>
and is evaluated at $\mathbf{k}_1 = \mathbf{p}_1, \mathbf{k}_2 = \mathbf{k} - \mathbf{p}_1, \mathbf{k}_3 = -\mathbf{p}_2, \mathbf{k}_4 = -\mathbf{k} + \mathbf{p}_2$, with $\psi = \phi_1 - \phi_2$ the difference between the azimuthal angles of $\mathbf{p}_1$ and $\mathbf{p}_2$ with respect to $\mathbf{k}$. The integration kernel for emission during radiation domination is given by

$$ \overline{I(u_1, v_1, k, \eta)I(u_2, v_2, k, \eta)} = \frac{1}{2(k\eta)^2} I_A(u_1, v_1)I_A(u_2, v_2) \times [I_B(u_1, v_1)I_B(u_2, v_2) + I_C(u_1, v_1)I_C(u_2, v_2)], \quad (4.31) $$

in terms of $I_{A,B,C}$ defined in eq. (4.22) and with $x = k\eta$. As for the disconnected contribution, it is numerically convenient to change the integration variables from $(u_i, v_i)$ to $(t_i, s_i)$, with

$$ \int_0^\infty dv_i \int_{|1-v_i|}^{1+v_i} du_i(\dots) = \frac{1}{2} \int_0^\infty dt_i \int_{-1}^1 ds_i(\dots), \quad (4.32) $$

where we did not assume symmetry between positive and negative $s$, to retain full generality in this case.

On general grounds, the properties of the trispectrum and the symmetries of the kernel eq. (4.31) enable one to split the connected contribution (4.29) into three inequivalent channels [51]. Hence, any trispectrum function of the unordered set $\{\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3, \mathbf{k}_4\}$, can be written as

$$ T_{\zeta}[\{\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3, \mathbf{k}_4\}] = \tilde{T}_{\zeta}[(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3, \mathbf{k}_4)] + 23 \text{ perm.}, \quad (4.33) $$

where the individual contributions $\tilde{T}_{\zeta}$ are not in general invariant under permutations of their arguments. Moreover, it can be conveniently written as

$$ T_{\zeta} = (T_{\mathfrak{s}} + T_{\mathfrak{t}} + T_{\mathfrak{u}}) + 7 \text{ perm.}, \quad \text{with} \quad \left\{ \begin{array}{l} T_{\mathfrak{s}} = \tilde{T}_{\zeta}[(\mathbf{k}_1, \mathbf{k}_2; \mathbf{k}_3, \mathbf{k}_4)] \\ T_{\mathfrak{t}} = \tilde{T}_{\zeta}[(\mathbf{k}_1, \mathbf{k}_3; \mathbf{k}_2, \mathbf{k}_4)], \\ T_{\mathfrak{u}} = \tilde{T}_{\zeta}[(\mathbf{k}_1, \mathbf{k}_4; \mathbf{k}_2, \mathbf{k}_3)] \end{array} \right. , \quad (4.34) $$

where the channels $\mathfrak{s}, \mathfrak{t}$ and $\mathfrak{u}$ correspond to the three unordered pairs $\{\{\mathbf{k}_1, \mathbf{k}_2\}, \{\mathbf{k}_3, \mathbf{k}_4\}\}$, $\{\{\mathbf{k}_1, \mathbf{k}_3\}, \{\mathbf{k}_2, \mathbf{k}_4\}\}$ and $\{\{\mathbf{k}_1, \mathbf{k}_4\}, \{\mathbf{k}_2, \mathbf{k}_3\}\}$ together with their “exchanged momenta” $\mathfrak{s} = |\mathbf{k}_1 + \mathbf{k}_2|, \mathfrak{t} = |\mathbf{k}_1 + \mathbf{k}_3|$ and $\mathfrak{u} = |\mathbf{k}_1 + \mathbf{k}_4|$ respectively, and where the seven permutations in (4.34) preserves the exchanged momentum of each channel. Furthermore one can show that the trispectrum-induced GW spectrum (4.29) can be written in terms of only three fundamental contributions corresponding to the seeds $T_{\mathfrak{s}}, T_{\mathfrak{t}}$ and $T_{\mathfrak{u}}$:

$$ \overline{\mathcal{P}_{h,c}} = 8(\overline{\mathcal{P}_{h,c}^{\mathfrak{s}}} + \overline{\mathcal{P}_{h,c}^{\mathfrak{t}}}) + \overline{\mathcal{P}_{h,c}^{\mathfrak{u}}}, \quad (4.35) $$

with $\overline{\mathcal{P}_{h,c}^{\mathfrak{s}}}$ simply corresponding to $\overline{\mathcal{P}_{h,c}}$ with $T_{\zeta}$ replaced by $T_{\mathfrak{s}}$, etc.

Computing from first principles the trispectrum generated in models relevant for GW astronomy is a difficult task, as these scenarios often involve a strong breaking of scale invariance as well as enhanced fluctuations that can jeopardize perturbative computations, see e.g. [119, 133, 135, 180, 292, 320–323]. In the following, we assume that the trispectrum is of the local $\tau_{\text{NL}}$ type:

$$ T_{\zeta}^{\text{loc}} = \tau_{\text{NL}} [(P_{\zeta}(\mathfrak{s})P_{\zeta}(k_1)P_{\zeta}(k_3) + 3 \text{ perm.}) + (\mathfrak{s} \leftrightarrow \mathfrak{t}) + (\mathfrak{s} \leftrightarrow \mathfrak{u})] . \quad (4.36) $$
We give more details on the computation of the spectrum in the presence of local NGs in appendix A.4.

Local NGs, typical of multi-field models, generically arise from the non-linear evolution of cosmological fluctuations on super-Hubble scales (see e.g. [38] for a review). Besides the $\tau_{\text{NL}}$ type, which emerges microscopically from the exchange of scalar particles through cubic interactions, the local trispectrum also acquires in general a $g_{\text{NL}}$ component, coming from contact quartic interactions. However, its momentum dependence is such that it does not contribute to the GW spectrum, and hence we can disregard it for our purpose. We stress that our choice, eq. (4.36), is a first methodological step motivated by simplicity, in particular because the momentum dependence of the trispectrum is fully characterized by the one of the power spectrum $P_{\zeta}(k)$, which is not the case in general (see [51] for a study of the impact of various trispectrum shapes on the GW spectrum). Let us also highlight a conceptual aspect. Several works in the literature consider a local ansatz in which the real-space curvature perturbation $\zeta(\boldsymbol{x})$ is expanded in powers of a Gaussian variable $\zeta_G(\boldsymbol{x})$ as

$$ \zeta = \zeta_G + \frac{3}{5} f_{\text{NL}} \zeta_G^2 + \dots, \qquad (4.37) $$

and compare the corresponding GW spectrum with the one obtained by keeping only the first term, fully characterized by the power spectrum $P_{\zeta_G}$, see e.g. [43, 44, 47, 48, 50, 54, 55, 286, 319, 324-327]. At leading order, such an expansion does lead to the trispectrum (4.36) with $P_{\zeta}$ replaced by $P_{\zeta_G}$, and $\tau_{\text{NL}} = (6f_{\text{NL}}/5)^2$. However, the nonlinear terms in eq. (4.37) also imply that the curvature power spectrum does not coincide with $P_{\zeta_G}$. Instead, keeping only the quadratic term shown in eq. (4.37) for definiteness, one finds the power spectrum

$$ P_{\zeta}(k) = P_{\zeta_G}(k) + \frac{1}{2} \left(\frac{6}{5} f_{\text{NL}}\right)^2 \int \frac{d^3\boldsymbol{p}}{(2\pi)^3} P_{\zeta_G}(p) P_{\zeta_G}(|\boldsymbol{k} - \boldsymbol{p}|). \qquad (4.38) $$

Hence, as described in [51], in this approach, one considers on similar grounds the impact of primordial NG on the SIGW, through the trispectrum, and the difference between a putative $P_{\zeta_G}$ to which we have no access, and the power spectrum of $\zeta$ which is anyway the only observable quantity. Again, we emphasize that the effects of non-linearities on the SIGW spectrum may not always be fully captured by the local ansatz, eq. (4.37). As a result, the predicted amplitude of the resulting SGWB could differ significantly [69], potentially being suppressed or enhanced by several orders of magnitude.

In our analysis, whose results are shown in section 6.5, we find it conceptually clearer to take as a benchmark the disconnected prediction from the purely Gaussian theory (4.16) with a given power spectrum $P_{\zeta}(k)$, which will take to be of the log-normal form (3.5), and to compare it with the addition of the non-Gaussian, connected, contribution, eq. (4.29) with trispectrum (4.36). Note that for the latter, the s-channel contribution vanishes as the corresponding $T_{\zeta}$ in eq. (4.29) does not depend on the azimuthal angle $\psi$. We are thus left with the two t and u contributions in eq. (4.35). Overall, we stress that the parameter to be constrained from observations is $\tau_{\text{NL}}$, which measures the non-Gaussian contribution to the SIGW spectrum coming from the trispectrum of curvature perturbations. On scales relevant to LISA, there is a priori no constraint on $\tau_{\text{NL}}$ except that it is positive in known concrete realizations of inflation. Its size is also a priori arbitrary, although from a theoretical perspective, perturbative control during inflation typically implies $\tau_{\text{NL}}P_{\zeta} < 1$.
## 5 Mock signal reconstructions with the SGWBinner and SIGWAY codes

This section outlines the analysis method employed in this work. Before presenting the LISA data model adopted in our analysis (section 5.1) and functionalities of the code (section 5.2), let us briefly illustrate the measurement of GWs with LISA.

The observatory will consist of three satellites ($\alpha = 1, 2, 3$) that orbit at the vertices of an approximately equilateral triangle with sides about 2.5 million kilometers long. Each satellite contains two Test Masses (TMs), whose positions are constantly monitored, and two lasers emitting toward the other satellites. By monitoring the fractional Doppler frequency shifts of photons traveling along the arms between satellites, LISA measures the relative displacements of the TMs. The path connecting two satellites is typically dubbed “link” and the single link measurement can be denoted as $\eta_{\alpha\beta}(t)$, where the laser emitted from the satellite $\beta$ at time $t - L_{\alpha\beta}/c$ is recorded at time $t$ in the satellite $\alpha$. These measurements are, however, dominated by laser frequency noise, which is expected to be several orders of magnitude greater than the required sensitivity [1]. To suppress this noise contribution, LISA will employ a post-processing technique called Time-Delay Interferometry (TDI) [328–336]. In practice, TDI can be understood as the operation of 3 × 6 matrix on the six single link measurements $\eta_{\alpha\beta}(t)$ [66, 337] that returns the three TDI channels where the laser frequency noise is strongly suppressed.

As in the previous studies using the SGWBinner code [6, 63, 64, 338, 339], in this work we assume for simplicity i) equal and static arm lengths and ii) equality of noise at each link. While, in reality, these hypotheses will not be perfectly satisfied,<fnref n="11" /> it has been shown that the signal reconstruction is almost unaffected by unequal (but static) arm length and unequal noise amplitudes [337, 342]. Under the equal and static arm length assumption, the so-called first-generation TDI variables suffice to achieve laser noise cancellation.<fnref n="12" /> In the {X, Y, Z} basis, they are expressed as

$$X = (1 - D_{13}D_{31})(\eta_{12} + D_{12}\eta_{21}) + (D_{12}D_{21} - 1)(\eta_{13} + D_{13}\eta_{31}), \quad (5.1)$$

with Y and Z being cyclic permutations of X. Here $D_{\alpha\beta}$ is the delay operator acting on any time-dependent function $x(t)$ as $D_{\alpha\beta}x(t) = x(t - L_{\alpha\beta})$ and we take $L_{\alpha\beta} = L = 2.5 \times 10^9$ m. For SGWB signal searches, it is convenient to combine the XYZ variables to obtain the so-called AET basis [344, 345], defined as

$$A = \frac{Z - X}{\sqrt{2}}, \quad E = \frac{X - 2Y + Z}{\sqrt{6}}, \quad T = \frac{X + Y + Z}{\sqrt{3}}, \quad (5.2)$$

which, in the limit of equal arms and equal noises, can be shown to have vanishing cross-correlations and simplify the likelihood computation. Moreover, due to its symmetric structure, the T channel strongly suppresses GW signals at small frequencies compared to instrumental noise. For this reason, the T channel can be treated as a quasi-null channel that is mostly sensitive to instrumental noise.<fnref n="13" />

<footnote n="11">¹¹With realistic orbits, LISA will not be perfectly equilateral and arm-lengths vary at the percent level [340] (see also appendix A of [341]).</footnote>

<footnote n="12">¹²To account for non-static arm lengths and the associated Doppler shifts, the second-generation TDI variables [333, 334, 336, 343] would be required.</footnote>

<footnote n="13">¹³The T channel does not remain a null channel in general with unequal and flexing arms and differing noise levels in the different spacecraft, although other quasi-null channels are available [346].</footnote>
## 5.1 Data streams from LISA TDI channels

We represent the three time-domain data streams as $d_i(t)$, where $i$ runs over the channels of the TDI basis. These quantities are real-valued functions defined on the interval $[-\tau/2, \tau/2]$ with $\tau$ being the duration of a data segment. The Fourier transforms of these data streams are then given by

$$ \tilde{d}_i(f) = \int_{-\tau/2}^{\tau/2} dt \, e^{2\pi i f t} d_i(t). \qquad (5.3) $$

Our central assumption is that all transients including loud deterministic signals and glitches in the noise are subtracted from the time stream through some appropriate methods within the LISA global fit scheme [57–59, 347–349].<fnref n="14" /> That is, as adopted in previous studies [6, 63, 64, 338, 339, 342], the data considered in our analysis only contain the stochastic contributions to the noise, $\tilde{n}_i^\nu$, and the stochastic signal $\tilde{s}_i^\sigma$ due to the unresolved binary signals and, possibly, the SGWB:

$$ \tilde{d}_i(f) = \sum_{\nu} \tilde{n}_i^{\nu}(f) + \sum_{\sigma} \tilde{s}_i^{\sigma}(f), \qquad (5.4) $$

where $\nu$ and $\sigma$ run over different noise and signal components, respectively. In the following, we will assume that all these components are stationary and obey Gaussian statistics with zero mean and variance given by

$$ \langle \tilde{n}_i^{\nu}(f) \tilde{n}_j^{\nu*}(f') \rangle = \frac{1}{2} \delta(f - f') P_{N,ij}^{\nu}(f), \quad \langle \tilde{s}_i^{\sigma}(f) \tilde{s}_j^{\sigma*}(f') \rangle = \frac{1}{2} \delta(f - f') P_{S,ij}^{\sigma}(f), \quad (5.5) $$

where we define the one-side power-spectral density (PSD) (for $i = j$) and cross-spectral density (CSD) (for $i \neq j$) of noise and signal components as $P_{N,ij}^{\nu}(f)$ and $P_{S,ij}^{\sigma}(f)$, respectively. Assuming all these components to be uncorrelated with one another, we obtain

$$ \begin{aligned} \langle \tilde{d}_i(f) \tilde{d}_j^*(f') \rangle &= \frac{1}{2} \delta(f - f') \left[ \sum_{\nu} P_{N,ij}^{\nu}(f) + \sum_{\sigma} P_{S,ij}^{\sigma}(f) \right] \\ &\equiv \frac{1}{2} \delta(f - f') [P_{N,ij}(f) + P_{S,ij}(f)], \end{aligned} \qquad (5.6) $$

where $P_{N,ij}(f)$ and $P_{S,ij}(f)$ are the total noise and signal PSDs and CSDs. By denoting the response functions for isotropic SGWB signals as $R_{ij}(f)$ (see refs. [64, 351] for expressions for this quantity), the SGWB (in either strain $S_h^\sigma(f)$ or abundance $\Omega_{GW}^\sigma(f)$) projected onto the data PSDs and CSDs can be expressed as

$$ P_{S,ij}(f) = R_{ij}(f) \sum_{\sigma} S_h^{\sigma}(f) = R_{ij}(f) \frac{3H_0^2}{4\pi^2 f^3} \sum_{\sigma} \Omega_{GW}^{\sigma}(f), \qquad (5.7) $$

where $H_0$ is the present Hubble constant and $h$ is the normalized one as $H_0/h \simeq 3.24 \times 10^{-18}$ 1/s. Note again that under the assumptions stated above in this section, one finds that $R_{ij}(f)$ is diagonal in the AET basis.

<footnote n="14">¹⁴See ref. [350] for the application of simulation-based inference to the SGWB search performed by LISA in the presence of transient signals, which goes beyond the framework of purely stochastic analysis.</footnote>
It is common practice to quantify the predicted primordial SGWB signal in terms of $h^2\Omega_{GW}(f)$; therefore, for later convenience, we define

$$P_{N,ij}^{\Omega}(f) = \frac{4\pi^2 f^3}{3H_0^2} P_{N,ij}(f). \qquad (5.8)$$

In the following, we provide more detailed descriptions of the noise PSDs in the AET basis and of the astrophysical foregrounds which are included in $h^2\Omega_{GW}^{\sigma}(f)$.

### 5.1.1 Instrumental noise

Our current knowledge of the LISA noise is based on the LISA Pathfinder [352] and laboratory tests. As a first approximation, the stochastic component of the noise in each TDI channel can be grouped into two effective components: “Optical Metrology System” (OMS) noise and TM noise. The former accounts for noise in the readout frequency, such as laser shot noise, while the latter models the noise sources causing accelerations of the TMs, e.g., by environmental disturbances. Introducing the transfer functions for these two noise sources $T_{ij,\alpha\beta}^{\nu}(f)$ (for details, see e.g. refs. [64, 337, 343, 353]), which project those contributions onto the TDI channels, the total noise PSDs and CSDs can be expressed as

$$P_{N,ij}(f) = \sum_{\nu} P_{N,ij}^{\nu}(f) = \sum_{\alpha\beta} \left[ T_{ij,\alpha\beta}^{\text{TM}}(f) S_{\alpha\beta}^{\text{TM}}(f) + T_{ij,\alpha\beta}^{\text{OMS}}(f) S_{\alpha\beta}^{\text{OMS}}(f) \right]. \qquad (5.9)$$

As customary in the literature, we assume stationary, Gaussian, and uncorrelated noises at each link with identical spectral shapes given by

$$S_{\alpha\beta}^{\text{TM}}(f) = 7.7 \times 10^{-46} \times A_{\alpha\beta}^2 \left(\frac{f_c}{f}\right)^2 \left[1 + \left(\frac{0.4\text{mHz}}{f}\right)^2\right] \left[1 + \left(\frac{f}{8\text{mHz}}\right)^4\right] \times \text{s}, \qquad (5.10)$$

$$S_{\alpha\beta}^{\text{OMS}}(f) = 1.6 \times 10^{-43} \times P_{\alpha\beta}^2 \left(\frac{f}{f_c}\right)^2 \left[1 + \left(\frac{2\text{mHz}}{f}\right)^4\right] \times \text{s}, \qquad (5.11)$$

where $A_{\alpha\beta}$ and $P_{\alpha\beta}$ represent the amplitudes of the TM and OMS noises in the different links. Moreover, we have introduced $f_c \equiv (2\pi L/c)^{-1} \simeq 19\text{mHz}$ representing the characteristic frequency of the detector. As mentioned above, we assume the noise amplitudes for all links to be identical, i.e. $A_{\alpha\beta} = A_{\text{noise}}$ and $P_{\alpha\beta} = P_{\text{noise}}$, and, following the ESA mission specifications [56], with fiducial values $A_{\text{noise}} = 3$ and $P_{\text{noise}} = 15$. In this case, the noise spectra reduce to $S_{\alpha\beta}^{\text{TM}}(f) = S^{\text{TM}}(f, A)$ and $S_{\alpha\beta}^{\text{OMS}}(f) = S^{\text{OMS}}(f, P)$. With $T_{ij,\alpha\beta}^{\nu}(f)$ in the equal arm length limit, the PSDs in the AET basis read

$$\begin{aligned} P_{N,AA}(f) &= P_{N,EE}(f) \\ &= 8\sin^2 x \left\{ 4 \left[ 1 + \cos x + \cos^2 x \right] S^{\text{TM}}(f, A_{\text{noise}}) + \left[ 2 + \cos x \right] S^{\text{OMS}}(f, P_{\text{noise}}) \right\}, \end{aligned} \qquad (5.12)$$

and

$$P_{N,TT}(f) = 16 \sin^2 x \left\{ 2 \left[ 1 - \cos x \right]^2 S^{\text{TM}}(f, A_{\text{noise}}) + \left[ 1 - \cos x \right] S^{\text{OMS}}(f, P_{\text{noise}}) \right\}, \qquad (5.13)$$

where we have defined $x \equiv f/f_c$. Here the CSDs vanish, i.e. $P_{N,ij}(f) = 0$ ($i \neq j$), so that the noise covariance matrix is diagonal.
### 5.1.2 Astrophysical foregrounds

Numerous weak and unresolvable signals from astrophysical sources will superimpose incoherently generating astrophysical SGWB [4, 354–360]. There are at least two guaranteed components in the LISA band. Below a few millihertz, the dominant contribution will come from Compact Galactic Binaries (CGBs) mostly composed of Double White Dwarfs (DWDs) [361, 362]. At higher frequencies, another contribution is expected from all the extragalactic compact objects including Stellar Origin Binary Black Holes (SOBBHs) and binary neutron stars (BNS) [363]. In the remainder of this section, we provide the templates for these foreground components implemented in the SGWBinner code that was recently used in refs. [6, 338, 339, 342].

**Galactic foreground.** This component represents the contribution from the unresolved sub-threshold mergers of CGBs that remain after the removal of loud signals from the population of CGBs in the galactic disk [364]. Due to the angular dependence of the response functions and LISA yearly orbit, this component exhibits an annual modulation. While, in principle, this characteristic can help distinguish the galactic component from other stationary contributions, e.g., by accounting for variations in each segment [341, 365, 366], we average over anisotropies, which leads to suboptimal (but conservative<fnref n="15" />) foreground extraction. Similarly, because this foreground is formed by the superposition of many unresolvable sources, it is expected to have Gaussian statistics. Recently, refs. [61, 371, 372] have called into question whether the populations entering the foreground are sufficient for the central limit theorem to apply at all frequencies, and imply a Gaussian description of the foreground may be biased. The non-stationarity of the foreground also in principle induces some non-Gaussianity when time-averaging.

Nevertheless, we use the empirical model from ref. [373], which describes the sky-averaged and Gaussian contribution by

$$h^2 \Omega_{\text{GW}}^{\text{Gal}}(f) = \frac{1}{2} \left( \frac{f}{1 \text{ Hz}} \right)^{2/3} e^{-(f/f_1)^\alpha} \left[ 1 + \tanh \frac{f_{\text{knee}} - f}{f_2} \right] h^2 \Omega_{\text{Gal}}, \quad (5.14)$$

where the value of $f_1$ and $f_{\text{knee}}$ depends on the total observation time $T_{\text{obs}}$ as

$$\begin{aligned} \log_{10}(f_1/\text{Hz}) &= a_1 \log_{10}(T_{\text{obs}}/\text{year}) + b_1, \\ \log_{10}(f_{\text{knee}}/\text{Hz}) &= a_k \log_{10}(T_{\text{obs}}/\text{year}) + b_k. \end{aligned} \quad (5.15)$$

The exponential factor $e^{-(f/f_1)^\alpha}$ accounts for the loss of stochasticity at higher frequency [373], while the last tanh term models the expected complete subtraction of CGBs signal at frequencies $f > f_{\text{knee}}$. In order to keep the notation compact, we define $\log_{10}(h^2\Omega_{\text{Gal}}) \equiv \alpha_{\text{Gal}}$. From ref. [373], we set the fiducial values $a_1 = -0.15$, $b_1 = -2.72$, $a_k = -0.37$, $b_k = -2.49$, $\alpha = 1.56$, $f_2 = 6.7 \times 10^{-4}\text{Hz}$ and $\alpha_{\text{Gal}} = -7.84$.

<footnote n="15">¹⁵Keeping track of anisotropic nature of the signal, requires the analysis to be time-frequency. To consistently work this out, one would also need to keep track of the non-stationary nature of the noise (see e.g. [61, 367–370]) and the presence of gaps in the data.</footnote>
**Extragalactic foreground.** The extragalactic foreground, arising from the incoherent superposition of all extragalactic compact object mergers, includes potential contributions from SOBBHs, BNSs, EMRIs, and DWDs in their inspiral phase. In this work, we focus on only the SOBBH+BNS contribution, leaving any potential EMRI and DWD contribution to future work. Recent studies suggest that extreme mass-ratio inspirals can largely contribute to the foreground but only in somewhat extreme population synthesis scenarios [357]. Extragalactic DWDs may also be more abundant than previously estimated, with a relevant impact on the extragalactic foreground [359, 374] which ongoing analyses are verifying [375]. In the lack of a firmer understanding, we assume these contributions to be below the foregrounds of galactic binaries and extragalactic SOBBHs and BNSs.

We now focus on what we will assume to be the dominant contribution, the SOBBH and BNS foreground. The vast majority of these signals cannot be individually resolved by LISA [376–378] and, for the most part individual detections are possible for the few multi-band sources [379] (see [380], for a more accurate study of such sources). The best estimates for the populations of these objects are based on observations from ground-based detectors [381, 382]. Due to the relatively uniform distribution of the sources and the limited angular resolution of LISA, this component can be well modeled as an isotropic SGWB signal with the power-law shape

$$h^2 \Omega_{\text{GW}}^{\text{Ext}}(f) = h^2 \Omega_{\text{Ext}} \left( \frac{f}{1\text{mHz}} \right)^{2/3}, \qquad (5.16)$$

where $h^2\Omega_{\text{Ext}}$ is the amplitude at 1 mHz. Recent observations by LIGO-Virgo-KAGRA collaboration estimate the magnitude of SGWB signal from SOBBHs and BNS as [381]

$$\Omega_{\text{Ext}} = 7.2_{-2.3}^{+3.3} \times 10^{-10} \text{ at } f = 25 \text{ Hz.} \qquad (5.17)$$

In order to keep the notation compact, we define $\log_{10}(h^2\Omega_{\text{Ext}}) \equiv \alpha_{\text{Ext}}$. Extrapolating this amplitude to the LISA band [358], yields the fiducial value $\alpha_{\text{Ext}} = -12.38$.

## 5.2 Analysis of the simulated data

In this section, we summarise the data analysis scheme implemented in the SGWBinner code (see refs. [63, 64] for more details). Let us start with the generation of simulated data. Given the effective observation time $T_{\text{obs}}$ and the number of segments $N_d$ (which define the duration of each segment $\tau = T_{\text{obs}}/N_d$), the code generates the data $\tilde{d}_i^s(f_k)$ ($s = 1, \dots, N_d$) segment-by-segment in the frequency-domain. For each frequency bin $f_k$ (spanning $[3 \times 10^{-5}, 0.5]$ Hz with spacing $\Delta f = 1/\tau$), $N_d$ Gaussian realizations of the signal, noise, and foregrounds are generated with zero mean and variances defined by their respective PSDs. These data are then averaged over segments to define $\bar{D}_{ij}^k \equiv \sum_{s=1}^{N_d} d_i^s(f_k) d_j^{s*}(f_k)/N_d$, which gives an estimate of the total power at all frequencies. The next step consists of coarse-graining the data using inverse variance weighting. This results in a coarser set of frequency bins $f_{ij}^k$ and a data set $D_{ij}^k$ with weights $n_{ij}^k$, retaining similar statistical properties of the original dataset. Similarly to refs. [6, 64, 338, 339, 342], we set $\tau = 11.4$ days ($\Delta f = 10^{-6}$ Hz), $N_d = 126$, and $T_{\text{obs}} = 4$ years in our analysis.
The likelihood employed in the code reads [64]

$$ \ln \mathcal{L}(D|\boldsymbol{\theta}) = \frac{1}{3} \ln \mathcal{L}_G(D|\boldsymbol{\theta}) + \frac{2}{3} \ln \mathcal{L}_{\text{LN}}(D|\boldsymbol{\theta}), \qquad (5.18) $$

with

$$ \ln \mathcal{L}_G(D|\boldsymbol{\theta}) = -\frac{N_d}{2} \sum_{i \in \{\text{AET}\}} \sum_k n_{ii}^k \left[ \frac{\mathcal{D}_{ii}^{th}(f_{ii}^k, \boldsymbol{\theta}) - \mathcal{D}_{ii}^k}{\mathcal{D}_{ii}^{th}(f_{ii}^k, \boldsymbol{\theta})} \right]^2, \qquad (5.19) $$

$$ \ln \mathcal{L}_{\text{LN}}(D|\boldsymbol{\theta}) = -\frac{N_d}{2} \sum_{i \in \{\text{AET}\}} \sum_k n_{ii}^k \ln^2 \left[ \frac{\mathcal{D}_{ii}^{th}(f_{ii}^k, \boldsymbol{\theta})}{\mathcal{D}_{ii}^k} \right], \qquad (5.20) $$

where the index k runs over the coarse-grained data points and $\mathcal{D}_{ii}^{th}(f, \boldsymbol{\theta})$ denotes the theoretical predictions for the data, depending on some parameters $\boldsymbol{\theta}$. The model can be further expressed as $\mathcal{D}_{ii}^{th}(f, \boldsymbol{\theta}) \equiv \mathcal{R}_{ii} h^2 \Omega_{\text{GW}}(f, \boldsymbol{\theta}_{\text{cosmo}}, \boldsymbol{\theta}_{\text{fg}}) + P_{N,ii}^{\Omega}(f, \boldsymbol{\theta}_{\text{n}})$, with $\boldsymbol{\theta}_{\text{cosmo}}, \boldsymbol{\theta}_{\text{fg}}$ and $\boldsymbol{\theta}_{\text{n}}$ denoting the signal, foreground, and noise parameters, respectively. Notice that the diagonality of the AET basis has been exploited, and no cross terms appear in the likelihood. Given some priors $\pi(\boldsymbol{\theta})$ for the parameters, the posterior distribution reads

$$ p(\boldsymbol{\theta}|D) \equiv \frac{\pi(\boldsymbol{\theta})\mathcal{L}(D|\boldsymbol{\theta})}{Z(D)}, \qquad (5.21) $$

where $Z(D)$ is the model evidence defined as

$$ Z(D) \equiv \int d\boldsymbol{\theta} \, \pi(\boldsymbol{\theta}) \mathcal{L}(D|\boldsymbol{\theta}). \qquad (5.22) $$

To compare the validity of two different models $M_i(\boldsymbol{\theta}_i)$, each characterized by a set of parameters $\boldsymbol{\theta}_i$, we can use the evidence $Z_i$ as a measure of the quality of the models given the data. The *Bayes factor* for two models $i, j$ is defined as $B_{ij} \equiv Z_i/Z_j$. This Bayes factor can then be compared to the Jeffreys' scale [383] to determine which model is favored by the data.

As key functionalities, the SGWBinner code offers i) model-agnostic signal reconstruction and ii) template-based signal reconstruction. The former fits the signal in each frequency bin using a power-law template, i.e. the signal parameters are

$$ \boldsymbol{\theta}_{\text{cosmo}} = \{\alpha_1, n_{\text{T},1}, \dots, \alpha_n, n_{\text{T},n}\}, \qquad (5.23) $$

with n denoting the number of bins. The number and width of the bins are dynamically adjusted as described in refs. [63, 64]. In practice, this method enables a preliminary identification of the spectral shape of the signal, which can guide the choice of the template for the template-based analysis. For the latter, the vector of parameters of the cosmological component $\boldsymbol{\theta}_{\text{cosmo}}$ corresponds to the template parameters.

In this work, we assume the fiducial noise and foreground parameters to be

$$ \boldsymbol{\theta}_{\text{n}} = \{A_{\text{noise}}, P_{\text{noise}}\}, \quad \boldsymbol{\theta}_{\text{fg}} = \{\alpha_{\text{Gal}}, \alpha_{\text{Ext}}\}, \qquad (5.24) $$

while we assume that the other foreground parameters are known,<fnref n="16" /> and $\boldsymbol{\theta}_{\text{cosmo}}$ is model dependent. Moreover, when fitting the simulated data, we use the same noise model applied

<footnote n="16">¹⁶The effect of loosening this assumption on the signal reconstruction has been discussed in ref. [342].</footnote>
to generate the data.<fnref n="17" /> Both for the noise and foreground amplitudes, we assume Gaussian priors centered on their fiducial values. For the former, we set the standard deviation to be 20% of the fiducial mean value. For $\alpha_{Gal}$ and $\alpha_{Ext}$, we set the standard deviation to be 0.21 and 0.17, respectively. To sample the parameter space the code relies on the Cobaya [384] inference framework. To facilitate template-based analysis specifically for SIGW signals, we develop the dedicated SIGWAY code. As detailed in appendix A, the SIGWAY code implements the parameterization of curvature perturbations discussed in section 3 and performs the numerical computation of SIGW signals.

Finally, the code also supports Fisher analysis. In practice, the Fisher Information Matrix (FIM) can be computed by the continuous integral over the frequency range, expressed as

$$F_{ab} \equiv T_{\text{obs}} \sum_{i \in \{\text{AET}\}} \int_{f_{\text{min}}}^{f_{\text{max}}} \mathrm{d}f \left. \frac{\partial \ln \mathcal{D}_{ii}^{th}}{\partial \theta^a} \frac{\partial \ln \mathcal{D}_{ii}^{th}}{\partial \theta^b} \right|_{\theta=\theta_{\text{fid}}}, \quad (5.25)$$

where $f_{min}$ and $f_{max}$ represent the detector's minimal and maximum measured frequencies, assumed to be $f_{min} = 3 \times 10^{-5}$ Hz and $f_{max} = 0.5$ Hz [56]. If non-trivial (log-)priors are included in the analysis, the code consistently adds their derivatives to eq. (5.25) to obtain the full FIM. The relative uncertainty on the reconstruction parameters can then be estimated from the covariance $cov_{ab} = F_{ab}^{-1}$. Given its computational efficiency, we also employed the FIM approach to assess the prospect of signal reconstruction with some level of accuracy. Note that in the case of SIGW signals, the FIM can be efficiently computed using the automatic differentiation feature of the JAX library [385] by applying the $\theta$ derivative in eq. (5.25) directly to $P_\zeta$, before it is integrated to yield $P_h$. We stress that the FIM formalism only works under the assumption that the likelihood is well approximated by a Gaussian distribution in the model parameters around the best fit (and that, when dealing with real data, the true values of the parameters lie within the region where the FIM is evaluated).

Finally, to complement the visualization of relative uncertainties in the parameter space, we will plot the signal-to-noise ratio (SNR) defined as

$$\text{SNR} \equiv \sqrt{T_{\text{obs}} \sum_{i \in \{\text{AET}\}} \int_{f_{\text{min}}}^{f_{\text{max}}} \left( \frac{P_{S,ii}^{\sigma}}{P_{N,ii}} \right)^2 \mathrm{d}f}, \quad (5.26)$$

which scales linearly with the signal amplitude.

# 6 Results

In this section, we summarise our main results presenting different analyses based on the SIGWAY code outlined in section 5.2 (see appendix A for more details). We adopt the three approaches discussed in section 3, namely i) binned spectrum agnostic approach; ii) template-based approach; iii) first principle USR model of inflation — all limited to the leading order SIGW and assuming an RD universe. We then consider specific examples including non-standard early universe evolution and non-Gaussianities.

<footnote n="17">¹⁷We note that any differences between the instrumental noise and the model could introduce bias. This issue will have to be closely monitored in future upgrades of the code.</footnote>
We report results for both $\Omega_{GW}$ and $P_{\zeta}$. In the former case, we include the noise curves as well as the foregrounds as discussed in section 5. Since the SIGW backgrounds we consider are emitted at very high redshift, when scales currently associated with mHz re-enter the Hubble sphere, they contribute to the energy budget in the early Universe and can affect cosmological observables as any other relativistic free-streaming component beyond the standard model. In particular, the SIGW contributes to the effective number of neutrino species as $N_{\text{eff}} \equiv 3.044 + \Delta N_{\text{eff}}^{\text{GW}}$, with $\Delta N_{\text{eff}}^{\text{GW}} = \rho_{\text{GW}}/\rho_{\nu,1}$ and $\rho_{\nu,1}$ is the energy density of a single neutrino species. Specifically, the total (integrated) GW abundance is $\Omega_{\text{GW}}h^2 \simeq 1.6 \cdot 10^{-6} (\Delta N_{\text{eff}}^{\text{GW}}/0.28)$ [10]. Measurements of the CMB [303] and Baryon Acoustic Oscillations (BAO) constrain $\Delta N_{\text{eff}} \le 0.28$ at 95% C.L. We report this bound for reference as shaded gray regions in the $\Omega_{\text{GW}}h^2$ plots.

Strong primordial density perturbations can lead to the copious formation of PBHs with masses of the order of the horizon mass $M_H = 1.3 \times 10^{-15} M_{\odot} [(k/\kappa_{\text{rm}})/\text{s}^{-1}]^{-2}$, where we kept track of the additional prefactor $\kappa_{\text{rm}} \equiv kr_m \sim \mathcal{O}(3)$ that relates the perturbation scale to the characteristic perturbation size $r_m$ at Hubble crossing [386–389]. Thus, the overproduction of dark matter in the form of PBHs in the asteroid mass range implies a bound $P_{\zeta} \le \mathcal{O}(10^{-2})$ [265, 390–393] on the scalar curvature perturbations and thus also on the strength of the SIGW in the mHz frequency band. The abundance and the mass distribution of PBHs depend on the shape of the curvature power spectrum, non-Gaussianities, and on the equation of state of the universe during their formation [29], so does the implied upper bound on SIGWs.

## 6.1 Binned spectrum method

In figure 4 we report the constraints obtained with the binned method (see section 3.1 and 4.4) when injecting no SIGW signal. This analysis forecasts the model-independent upper bounds on both the SIGW energy density spectrum (left panel) and the primordial curvature power spectrum (right panel) in case of no SGWB detection at LISA: this only relies on observational data, without assuming a specific signal model. For this analysis, we assume that the spectrum is divided into $N = 15$ bins. The free parameters in this model are

$$ \theta_{\text{cosmo}} = \{A_1, \dots, A_{15}\}. \qquad (6.1) $$

In the left panel of figure 4 we indicate the posterior predictive distribution for $\Omega_{\text{GW}}$ with the shaded light blue region, denoting the 95% credible interval (CI). The upper bound effectively reflects the LISA sensitivity, which falls around two orders of magnitude below the noise components in the AA channel because of the long observation time $T_{\text{obs}} = 4\text{yr}$ (see section 5.2). In the right panel, we then show how the LISA sensitivity translates into the $P_{\zeta}$ parameter space. The figure displays the upper bounds on $P_{\zeta}$ across the range of momenta $k$ considered, which is $k \in [1.26 \times 10^{-4}, 6.28]/\text{s}$. The posterior saturates at the lower edge of the prior on the amplitude parameters $A_i > 10^{-8}$ reflecting the absence of detectable power beyond the noise level. The posterior predictive bands illustrate that the method can constrain $P_{\zeta}$ across several orders of magnitude in $k$. In the most sensitive range, this bound reaches $P_{\zeta} \lesssim 2 \times 10^{-4}$.

This sensitivity is sufficient for probing a wide range of viable scenarios for asteroid mass PBHs. In particular, a non-detection of a SIGW by LISA would close the asteroid mass
<figure><img src="image_6.png" /></figure>

**Figure 4.** Posterior predictive distribution for both $\Omega_{\text{GW}}h^2$ (left panel) and $P_\zeta$ (right panel). We represent the binned reconstruction of section 6.1 in the case without an injected signal. The posterior saturates the lower bound of the prior for the amplitudes $A_i$, due to the absence of a resolvable signal. We therefore can only set upper bounds on both $\Omega_{\text{GW}}h^2$ and $P_\zeta$. The light blue line shows the 95% credible intervals, while the pale black lines individual realisations of a signal sampled from the parameters' posterior distribution. The upper bound from $\Delta N_{\text{eff}}$ is shown with a gray shading.

window for PBH dark matter formed from the collapse of moderately non-Gaussian curvature fluctuations, including models of PBHs from first-order phase transitions [394–396]. However, as with $\mu$-distortion bounds on heavy PBHs [397–400], extremely strong non-Gaussianities could enhance PBH production and potentially allow evading these constraints. Such extreme scenarios and their theoretical consistency should be studied case by case.

In figure 5 we report the constraints obtained with the binned method when injecting the benchmark SIGW signal derived from the single field USR model of section 2.4. The curvature power spectrum has a BPL shape (see eq. (3.7)), with a peak at around $P_\zeta \sim 2 \cdot 10^{-3}$. Again, we use a template with $N = 15$ bins. In the left panel of figure 5, we show the injected SIGW signal (blue dashed line), along with the posterior predictive distribution (light blue band). Although the low number of bins reduces the frequency resolution of our model compared to the one achieved by LISA, the SIGW spectrum is well reconstructed, reaching a precision of the order of a few percent around the peak. At the edges of the observable range of frequencies, the blue bands widen up indicating a poor constraining power on the tail regions. The right panel of figure 5 indicates the SIGW bounds translate into four bins being well constrained in the range $k \sim [10^{-2}, 10^{-1}]$/s with around $\mathcal{O}(10)\%$ precision, while the other ones being subject to an upper bound of similar amplitude as in figure 4.

One could in principle enhance the frequency resolution by using a template with a larger number of bins, at the cost of drastically increasing the computational cost of the Bayesian MCMC inference. In appendix B we discuss these issues in more detail.

## 6.2 Template based method

In this section, we present a forecast on reconstructing the SIGW signal using a template-based method, addressing different scenarios discussed in section 3.2.
<figure><img src="image_7.png" /><figcaption>Figure 5. Same as figure 4, but simulating the observation of a signal obtained in the benchmark USR model scenario. The quantity $P_ζ$ is reconstructed with the model-independent binning method with 15 bins. The blue band in the upper panels shows the 90% (symmetric) credible interval, while the blue band in the bottom shows the 95% upper bound on the residuals.</figcaption></figure>

### 6.2.1 Smooth spectra

**Lognormal scalar spectrum.** The first signal injection we consider is a log-normal shape of $P_ζ$ as defined in eq. (3.5) with the benchmark values defined in eq. (3.6), which we report here $log_{10} A_s = -2.5$, $log_{10} Δ = log_{10}(0.5)$, $log_{10}(k_*/s^{-1}) = -2$. This produces a loud signal with the typical double peak being fully within the LISA band. This choice gives us a concrete measure of the precision that is achievable in measuring the SIGW background in this scenario for a reasonably loud signal. The free parameters for this approach are

$$ \theta_{\text{cosmo}} = \{\log_{10} A_s, \log_{10} Δ, \log_{10} (k_*/s^{-1})\}. \qquad (6.2) $$

Figure 6 shows the posterior distribution for each parameter of the LN template, alongside the ones describing the galactic/extragalactic foregrounds. We omit in these plots the posterior distributions for the noise parameters, as they are weakly correlated with the others in all cases. The injected values are indicated with a dashed gray line. As we can see, due to the relatively high SNR of the injected signal, the parameters of this template are very accurately reconstructed, with $P_ζ$ being reconstructed with a relative error of a few percent at its peak. The correlation between $A_s$ and $Δ$ originates from the definition of $P_ζ$. As customary in the literature [263], the amplitude at the peak is $A_s/Δ$, while $A_s$ is the integrated power spectrum $\int_{-\infty}^{\infty} P_ζ d \log k = A_s$. Therefore, $A_s$ enters the power spectrum only through the ratio $A_s/Δ$, thus the positive correlation between the two parameters. Defining $A_s$ to be the peak amplitude would avoid this degeneracy. The correlation of $A_s$, $k_*$, and $α_{\text{Gal}}$ is instead specific of our choice of fiducial parameters. As can be seen from the posterior predictive distribution in the left panel in figure 7, the injected signal is close enough to the galactic
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Corner plot with the posterior distribution for an injected LN spectrum $P_ζ$, as defined in eqs. (3.5) and (3.6). $k_*$ is expressed in 1/s units.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Same as figure 4, for an injected log-normal spectrum, as defined in eqs. (3.5) and (3.6). In the left panel, we arbitrarily cut the posterior predictive where the signal falls below $\Omega_{GW} \lesssim 10^{-15}$.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 8.** From left to right: relative uncertainties of each of the parameters of the LN template, computed with an FIM forecast injecting an LN $\mathcal{P}_\zeta$ with fixed amplitude $\log_{10} A_s = -2.5$ and varying $k_*$ and $\Delta$. Black (purple) contours show uncertainties without (with) astrophysical foregrounds. The white line indicates SNR values.</figcaption></figure>

foreground that a slight increase in $A_s$ with a decrease in $k_*$ can be compensated by a small decrease in the background amplitude $\alpha_{\text{Gal}}$. We expect that these correlations fade away with a larger injected $k_*$, when the signal and the galactic background are more distinct.

Figure 7 shows the posterior predictive distribution for the SIGW (left panel) and curvature power spectrum (right panel). We see that the signal reconstruction achieves better than percent uncertainties on $\Omega_{\text{GW}}$ and $\mathcal{P}_\zeta$ around the peak. The uncertainty on the low-frequency tail of $\Omega_{\text{GW}}$ saturates at around a few percent, due to the universal behavior of the causality tail [401] sufficiently deep in the IR. See e.g. appendix B of [6] for a discussion of logarithmic corrections to the IR tail of the SIGW spectrum. The uncertainty on the tails of $\mathcal{P}_\zeta$ remains low even at very small values, because of the rigid assumption about the LN template.

In figure 8, we scan the parameter space in $k_*$ and $\Delta$ estimating the relative uncertainties on all LN parameters using the FIM method. We fix the amplitude of $\mathcal{P}_\zeta$ to $\log_{10} A_s = -2.5$. The SIGW amplitude scales like $\Omega_{\text{GW}} \sim A_s^2$, and in the high SNR limit we expect the uncertainties on the parameter to scale inversely $\sim 1/A_s^2$. The results highlight the great sensitivity that is achievable on a SIGW background if the peak lies around the peak sensitivity of LISA, $k_* \sim 10^{-3}-10^{-1} \text{s}^{-1}$. In that range, the width $\Delta$ for an LN scalar spectrum can be measured with an accuracy of order 10% or better if $\Delta \lesssim \mathcal{O}(1)$. Notice from figure 8 that the
<figure><img src="image_11.png" /><figcaption>Figure 9. Same as figure 6, but for a recovered BPL curvature power spectrum assuming an injected signal motivated by the USR benchmark model. $k_*$ is expressed in 1/s units.</figcaption></figure>

purple contours, marking the sensitivity on the primordial SIGW background accounting for astrophysical foregrounds, degrade when $k_*$ coincides with the expected peak of the white-dwarfs (WD) galactic foreground, and the two GW backgrounds are less distinguishable [402].

**Broken power law.** The second injected signal is a BPL that is derived from the USR model discussed in section 2.4 and 3.3, with input parameters as in eqs. (3.7) and (3.8): $\log_{10} A_s = -2.71$, $\log_{10}(k_*/s^{-1}) = -1.58$, $\alpha = 3.11$, $\beta = 0.221$, $\gamma = 1.25$. We reconstruct the signal using a BPL template. The results of the reconstruction of the signal using the USR model will be discussed below. The free parameters for this approach are

$$ \theta_{\text{cosmo}} = \{\log_{10} A_s, \log_{10} (k_*/s^{-1}), \alpha, \beta, \gamma\}. \qquad (6.3) $$

In figure 9 we show a corner plot of the reconstructed parameters of the broken power-law template, while figure 10 displays the posterior predictive distribution for the SIGW (left panel) and curvature power spectrum (right panel). While the amplitude of the $P_\zeta$ peak $A_s$ and the peak position $k_*$ are well reconstructed, the infrared (IR) spectral index $\alpha$ and the smoothing coefficient $\gamma$ are poorly constrained, and $\alpha$ and $\gamma$ appear to be very degenerate as the IR tail of the signal is hidden by the galactic foreground. This can be seen in the corner plot of figure 9 as well as from the right panel of figure 10, where the slope for $k < k_*$ has a large uncertainty. In order to reconstruct $\alpha$ and $\beta$ with some precision, one would need to
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Same as figure 4, but for an injected signal from a BPL $P_ζ$, recovered using the BPL template. The injected parameters are motivated by the USR benchmark model.</figcaption></figure>

be sensitive to the tails of the signal outside the peak region. This would be only possible for a much larger signal, which nevertheless would have to compete with stringent bounds from PBH overproduction. With the relatively low SNR injected here, changes in the tilt can be traded for a smoother turnover around the peak, and vice versa. Furthermore, there is a residual correlation between $β$ and $γ$, although it is less pronounced. Biases in $A_s$, $β$, and $α_{Ext}$ are also evident due to their degeneracy in the high-frequency tail. Specifically, there is a tendency to reconstruct higher foreground values compared to the UV part of the SIGW spectrum. Importantly, we have verified that this bias is not an artifact introduced by the additional degeneracies induced by $γ$, which is correlated with both tilts and amplitude. This conclusion is supported by tests we ran with $γ$ fixed to its injected value.

Notably, this bias does not appear when the same signal is reconstructed using the USR model. The USR model is inherently less flexible, and its UV tilt is better constrained, thereby mitigating the impact of degeneracies.

In figure 11 we show FIM estimates of uncertainties on the tilt parameters depending on the injected BPL shape and $k_*$.$k_*$ can be independently resolved only if one observes the tail of the signal with a sufficiently high SNR. If one only observed the SIGW peak, a shallower (steeper) tilt can be traded off for a smoother (faster) transition. We see this happening in figures 9 and 10, where large degeneracies exist between $γ$ and the other parameters. As a consequence, these posterior can hardly be approximated with the FIM. Furthermore, in part of the parameter space explored in figure 11, we would obtain ill-conditioned FIM because of this degeneracy. In order to avoid this, only in this case we remove $γ$ from the parameters of the model and fix it to the injected value. This should be considered a slightly less flexible model for $P_ζ$, and an optimal choice should be made when performing a real GW data inference.
<figure><img src="image_13.png" /><figcaption>Figure 11. Uncertainty on α (left panel) and β (right panel) computed with the FIM approach for the broken power law $P_{\zeta}$. The remaining parameters injected are those of the benchmark scenario in eq. (3.8). We fix γ to the injection to remove the degeneracies with the two tilts in the small SNR regions of the parameter space.</figcaption></figure>

The left panel of figure 11 shows the absolute uncertainty achievable on α in the $(k_*, \alpha)$ plane, and the right panel shows the same for β. We checked that these uncertainties do not depend on the injected value for the other tilt parameter. The slope of the IR tail of the GW spectrum is only mildly dependent on α when $\alpha \gtrsim 1.5$, as discussed before, so most of the sensitivity comes from the signal in the frequency range around the peak. For this reason, the uncertainty on α reduces to 0.1 or better only if the SIGW is well within the LISA range, and on the right of the galactic WD foreground ($10^{-2} < k_*/s^{-1} < 10^{-1}$). Still, it is very interesting to notice that steep values of the tilt, higher than $\alpha \simeq 1.5$, which are fully covered by the causality tail $\Omega_{GW} \sim f^3$ [401, 403–405] in the SIGW spectrum (up to log-corrections), can be well constrained, as information on the tilt is still retained in the shape of the double peak feature of the signal close to the dominant peak, provided no large degeneracies with other parameters controlling the shape of the peak are present. The sensitivity to β (right panel of figure 11) is instead much better, as it determines the UV slope $\Omega_{GW}(f) \sim f^{-2\beta}$. Therefore, β cannot be measured with an uncertainty smaller than 0.1 only if the SIGW lies outside LISA’s peak range ($k_* > 10^{-1} s^{-1}$) or if $\beta \gtrsim -2$, where the SIGW background falls too quickly in the UV.

### 6.2.2 Spectra with oscillations

**Turns in multi-field inflation.** As a benchmark example of a primordial feature in the power spectrum, we analyze a signal arising from turns in multi-field space as introduced in section 3.2.2. The free parameters for this analysis are

$$ \theta_{\text{cosmo}} = \{\log_{10} A_s, \log_{10}(k_*/s^{-1}), \delta, \eta_{\perp}, F\}. \qquad (6.4) $$

In figure 12 we show the posterior distributions for key parameters governing the sharp-turn scenario in multi-field inflation, along with the foreground parameters. The injection assumes the benchmark scenario where the parameters controlling the template (3.9) are fixed as in eq. (3.12). We see that in this case, the parameter F is constrained to be close
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Same as figure 6, but for a simulated signal motivated by the multi-field scenario with sharp turns from (3.9). $k_*$ is expressed in $1/s$ units.</figcaption></figure>

to unity with better than percent precision, showing the high sensitivity to the template oscillations. In this case, the signal amplitude and central scale $k_*$ are weakly correlated, while the former is still positively correlated to both $\delta$ and $\eta_{\perp}$ which control the enhancement factor. Due to the ideal location of the SIGW peak, we also observe weak correlations between the foreground parameters and the signal parameters.

Figure 13 shows the posterior predictive distribution for $\Omega_{\text{GW}}$ and $P_{\zeta}$. In this example, the main peak of the SGWB lies within the LISA sensitivity band and above both astrophysical foregrounds. As a result, both the shape and amplitude of the peak in $\Omega_{\text{GW}}$, along with the $O(20\%)$ modulations, are reconstructed at the percent level. Since the frequency of these modulations is linked to the oscillations in $P_{\zeta}$ through the assumed thermal history at horizon re-entry, the oscillations are also reconstructed with high accuracy — see the right bottom plot. In this fortunate case, it would be possible to pinpoint the duration and strength of the field-space turn, as well as the inflationary time scale of the phenomenon. The latter is related to the oscillation frequency, as it is customary from the sharp feature phenomenon, while the former two can be disentangled by combining the peak amplitude, its location, and the frequency of the modulations.

Finally, the results of a Fisher analysis, highlighting the uncertainty in the parameter $F$ associated with the oscillatory behavior, are presented in figure 14. There, we vary the power spectrum amplitude $A_s$ and the position of the main peak, while keeping the other parameters fixed to the benchmark values discussed just above. This simplification is useful for illustrative purposes, as the parameters in the current model are not independent. Notably,
<figure><img src="image_15.png" /><figcaption>Figure 13. Same as figure 4, for the signal generated through a multi-field scenario with sharp turns from eq. (3.9).</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 14. Fisher analysis for the oscillation template from multi-field inflation with turns. We vary $k_*$ and $A_s$ while keeping the remaining parameters fixed to the benchmark values (3.12).</figcaption></figure>

when the signal is centered near the LISA sweet spot at $\log_{10}(k_*/s^{-1}) \simeq -2$, the oscillations can be accurately detected even with a moderate enhancement such as $\log_{10} A_s \simeq -3$.

**Rapid transitions between SR and USR phases.** The other injected spectrum with oscillatory features is characteristic of single field models with fast transitions from an SR to a USR phase described by eqs. (3.13), (3.14). The free parameters for this analysis are

$$ \theta_{\text{cosmo}} = \{\log_{10} A_s, \log_{10}(k_*/s^{-1}), \nu_1, \nu_{II}, F\}. \qquad (6.5) $$
<figure><img src="image_17.png" /><figcaption>**Figure 15.** Same as figure 6, but for a simulated signal from eq. (3.14). $k_*$ is expressed in $1/s$ units. Note the injected value $F = 1$ is not visible sitting at the edge of the plot.</figcaption></figure>

We consider the benchmark scenario with the input parameters listed eq. (3.18): $\log_{10} A_s = -2.58$, $\log_{10}(k_*/\text{s}^{-1}) = -2.02$, $\nu_I = 1.95$, $\nu_{II} = 1.61$, $\gamma = 1.67$, $F = 1$. We perform the MCMC Bayesian inference modelling of the signal using the template (3.14), which allows us to turn on the oscillations smoothly by varying the parameter $F$ from 0 to 1. The value $F = 0$ corresponds to a featureless BPL similar to the one considered above. Furthermore, we fix $\gamma$ as the strong degeneracy between $\alpha$ and $\gamma$ (see figure 9) makes sampling challenging and since our main concern lies in constraining $F$.

In figure 15 we show the posterior distribution for each parameter of the signal and foregrounds. First of all, we see the parameter $F$ controlling the relevance of the oscillations over the smooth BPL is very tightly constrained around unity. This tells us the presence of oscillations can be resolved with high accuracy for such a high-SNR signal. The sensitivity to oscillations is mainly driven by the dominant peak, as we will discuss in the following. We also find tight correlations between the parameters, which are non-trivially connected in the signal template (3.13). In particular, we observe a strong correlation between the BPL tilts and the amplitude, due to the large impact of the former on the overall amplitude of the dominant peak. The negative correlation in the $(A_s, k_*)$ plane is probably induced by the way the dominant peak, contributing to most of the SNR, can be adjusted, as one could lower the characteristic scale by enhancing the amplitude.

In figure 16 we show the posterior predictive distribution for both $\Omega_{\text{GW}}$ and $P_\zeta$. The presence of a dominant peak at scales around $k_*$ in the right plot leads to a distinctive large enhancement of the SIGW signal around peak frequencies seen in the left plot. Additional oscillations in the SIGW spectrum can be observed at larger frequencies, although the
<figure><img src="image_18.png" /><figcaption>Figure 16. Same as figure 4, for the signal generated through a fast transition from SR to USR from (3.14). Both $P_ζ$ and $Ω_{GW}$ are reconstructed very well.</figcaption></figure>

second-order emission soon washes out further oscillations in the UV tail. The residuals of $Ω_{GW}$ show that the IR tail of the signal is reconstructed at around the percent level, with a flat behavior due to the causality tail dominating the IR. On the other hand, the relative deviation grows larger than $O(10)\%$ percent in the UV part of the plot, due to the finite precision at reconstructing $ν_{II}$. The best accuracy is obtained around the peak, as expected. Correspondingly, in the right plot, $P_ζ$ is reconstructed with better than percent accuracy around the peak, while the reconstruction degrades in both tails. The oscillations are reconstructed for a few cycles in $k$, while they are lost in the UV as seen in the bottom panel showing the relative deviation from the injected signal. The envelope of the out-of-phase oscillations behaves following the underlying power-law tail. Note that the tails are reconstructed better than in the pure BPL scenario, as in this template (3.14), the shape of the dominant peak also brings information on the parameter $ν_{II}$ controlling the tails.

We perform a Fisher analysis focusing on the uncertainty on $F$, by varying $log_{10} A_s$ and $log_{10} k_*$, see figure 17. The oscillations are very well recovered, to $O(10^{-2})$ uncertainty, if the peak of the signal falls within the LISA band. This is because the oscillations in $P_ζ$ translate into oscillations mainly around the peak in the SGWB, as visible in figure 16.

## 6.3 Single field USR inference

In figure 18 we show the reconstruction capability of the USR model parameters of section 2.4, obtained by running a MCMC Bayesian inference using the USR inflationary model with free parameters

$$ \theta_{\text{cosmo}} = \{\lambda, v, b_l, b_f\}, \qquad (6.6) $$
<figure><img src="image_19.png" /><figcaption>**Figure 17.** Absolute uncertainty on $F$ estimated using the FIM for the oscillation template from a sharp transition between SR and USR. We vary $\log_{10} k_*$ and $\log_{10} A_s$ while keeping the remaining parameters fixed to the benchmark values (3.18).</figcaption></figure>

controlling the inflaton potential. We assume that the CMB scale crosses the Hubble sphere $N = 58$ e-folds before the end of inflation. We therefore avoid modelling the reheating era, and postopone its inclusion for future work (see e.g. [406]). The input values determining the injected signal were introduced in eqs. (2.4) and (2.6), but we report them here for convenience: $\lambda = 1.47312 \times 10^{-6}$, $v = 0.19688$, $b_l = 0.71223$, $b_f = 1.87 \times 10^{-5}$.

We understand the results as follows. The height of the peak in $P_\zeta$ is proportional to $\lambda$ and it is also sensitive to the tuning of $b_f$.<fnref n="18" /> This results in the negative correlation between $\lambda$ and $b_f$. The self-coupling $\lambda$ can be constrained, even though with only $O(1)$ precision, because it controls the slope of the potential and therefore the SR parameter $\eta_H$ before the inflection point, which determines the growth of $P_\zeta$ before the peak as discussed below eq. (3.7). The parameters $b_l$ and $v$ appear to be strongly correlated, meaning that the linear term in $b_l/v$ in eq. (2.4) gives the dominant dependence on $b_l$ in the potential. The galactic background is well reconstructed due to its large magnitude, while the extragalactic one is completely hidden by the USR signal.

Figure 19 shows the posterior predictive means in $\Omega_{\text{GW}}h^2$ and in $P_\zeta$. It is interesting to compare the right panel of this figure with that of figure 10, which is obtained with the same injected signal but a different template for the reconstruction. In the present case, the spectrum of scalar perturbations $P_\zeta$ is reconstructed with excellent precision, even if LISA is sensitive only to the peak. This comes from the fact that the spectrum for the USR model has a universal slope $\sim k^4$ in the IR, whereas the IR slope is a free parameter for the BPL model.

The relatively large uncertainty on the overall potential amplitude $V(\phi)$ in figure 20 is due to the degeneracies between the overall scale $V_0 \sim \lambda v^4$ and the parameter $b_f$ controlling the enhancement. As we are only constraining the enhanced part of the spectrum, there is a tight correlation between $\lambda$ and $b_f$. Adding information from CMB data in the inference would reduce this uncertainty by adding an independent constraint on $V_0$.

<footnote n="18">¹⁸With other potential parameters fixed, we did find the approximate behavior $P_\zeta \propto (1 - b_f/b_{f,*})^{-n}$ in the parameter region supporting peaked $P_\zeta$. Here, $b_{f,*}$ and $n > 2$ are parameters that depend on the remaining parameters of the potential. Such scaling is observed in other models [128].</footnote>
<figure><img src="image_20.png" /><figcaption>Figure 18. Same as figure 6, but for the USR reconstruction of the benchmark scenario.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Figure 19. Same as figure 4, but for the USR reconstruction of the benchmark scenario.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>**Figure 20.** Posterior predictive distribution of the potential $V(\phi)$ for the USR reconstruction of the benchmark scenario.</figcaption></figure>

**Comparison between different methods.** We can compare the performance of different methods when fitting the same injected signal, which is taken to be the USR benchmark scenario. In figure 21 we show the upper bound at 95% C.L. on the relative difference between the posterior predictive distribution and the injected signal for the binned, template-based, and ab initio USR approaches.

We observe that the binned method provides a competitive constraint on $\Omega_{GWh}^2$ in the central frequencies close to the peak (barring oscillations induced by the poor resolution associated with choosing 15 bins). However, the constraint quickly degrades at both ends, due to the unconstrained curvature spectral amplitude there. The template-based method (assuming a priori a SIGW from a BPL scalar power spectrum) improves the reconstruction of the tails, but results in an overall loss of precision of a factor $O(6)$ with respect to the posterior predictive derived assuming the USR scenario. This is most probably due to the larger number of parameters in the BPL template compared to the USR model, and the known degeneracy between $\gamma$ and the two tilts around the peak.

Also, $P_\zeta$, shown in the right panel, is most tightly constrained when using the USR model, showing the effectiveness and robustness of our analysis pipeline in reconstructing the injected high SNR signal. The BPL template gives an intermediate result on the IR tail, which, however, degrades much faster than USR, due to the limited information on the tail of the SIGW, which is mostly controlled by the causality tail. Finally, the binned method gives a worse reconstruction of both the IR and UV tail, due to the small information within the LISA band about these regimes, and the independence assumed in this method between the central bins (best constrained) and the ones on the sides.

Overall, the binned method proves to be a powerful approach to explore the interpretation of a primordial background at LISA within a more agnostic approach. The comparison with
<figure><img src="image_23.png" /><figcaption>Figure 21. Comparison of the residuals between the three different methods for recovering the signal injected assuming the benchmark USR model.</figcaption></figure>

specific SIGW templates does not significantly outperform the binned method for the range of frequencies around the peak, which are the best constrained by LISA. However, consistently with expectations, adopting the correct USR model provides greater accuracy in capturing the features of the signal, leading to a more precise reconstruction. These findings demonstrate the power of also adopting inference analyses based on explicit ab initio models (of which USR is just an example) that could outperform traditional template-based approaches. This, of course, assumes one can identify the best early universe model through model comparison. We will come back to discussing how to compare different scenarios in section 7.

## 6.4 Non standard thermal histories

Using information on the SIGW spectrum, LISA would be able to challenge the vanilla assumption that the SIGW was emitted during a RD era. As discussed in section 4, the kernels entering the computation of the SIGW spectrum bring information about the equation of state around the epoch of SIGW emission.

**A sudden transition from eMD era to the RD era.** We exemplify this case by showing how LISA can constrain the SIGWs emitted within an alternative thermal history by considering an early period of matter domination (eMD). We further assume sudden reheating, as introduced in section 4.3. As we discussed, during this eMD epoch $\Phi$ does not decay, leading to an enhancement of the SIGW spectra around the scale $k \gtrsim 1/\eta_R$. We take as a benchmark a nearly scale-invariant spectrum $P_\zeta$, with a cutoff at placed at $k_{max}$. We simplistically describe the spectrum as

$$P_{\zeta}(k) = A_s \Theta(k_{\max} - k), \qquad (6.7)$$

with benchmark parameters

$$A_s = 2.1 \times 10^{-9}, \quad \eta_R = 2000 \text{ s}, \quad k_{\max} = 0.06 \text{ s}^{-1}. \qquad (6.8)$$

Therefore the free parameters to examine in this case are

$$\theta_{\text{cosmo}} = \{A_s, k_{\max}, \eta_R\}. \qquad (6.9)$$
<figure><img src="image_24.png" /><figcaption>**Figure 22.** Left panel: same as figure 6, but for the case of a nearly scale-invariant power spectrum and a sudden eMD to RD transition. $k_{max}$ and $\eta_R$ are expressed in units of s$^{-1}$ and s, respectively. Right panel: corresponding posterior predictive distribution for $\Omega_{GW}h^2$. The notation used matches the one in figure 4.</figcaption></figure>

Our phenomenological parametrization should be regarded as a toy model, with the UV cut-off scale $k_{max}$ introduced to ensure perturbativity, as assumed when computing the SIGW. For this reason, given the fact that the energy density contrast grows linearly with the scale factor during a MD era, i.e. $\delta\rho/\rho \propto a$, one can associate $k_{max}$ as the scale at which the power spectrum of density contrast becomes unity, i.e. $P_\delta(k_{max}) = 1$ [312, 313, 407], although our actual choice is slightly more restrictive. One then can easily understand why $k_{max}$ depends on the scale we are probing, as the source is largest for modes that spent the most time within the horizon during the eMD era. While we do not model the non-linear part of the spectrum, it may lead to further observational signatures [220, 408–410]. Finally, the template (6.7) can be made more realistic by introducing a smooth cut-off.

The parameters of this template are accurately reconstructed as shown in figure 22 (left panel). We notice that $\eta_R$ and $A_s$ are strongly correlated, as the duration of the MD signal directly controls the growth of perturbations emitting SIGWs, which is therefore degenerate with the primordial amplitude. The cut-off scale is also strongly constrained, with a marginalised posterior distribution which is flat within a narrow range of scales corresponding to the resolution adopted in this forecast. It should be kept in mind, however, that more realistic spectra would feature a smoother drop-off, alongside a contribution from non-linear scales not included here, thus jeopardizing the relevance of the constraining power on $k_{max}$. In the right panel of figure 22, we show the posterior predictive distribution for the SIGW. We find the reconstruction to be accurate up to the cutoff scale (better than a few %). The SIGW spectrum is reconstructed well in the large-scale approximation, with the resonant amplification improving accuracy by an order of magnitude. The resonant peak
is reconstructed with a larger accuracy due to its milder model dependence and due to its tilt being controlled by the resonant conditions (see discussion around eq. (4.23)). Moreover, that part of the signal appears with a larger SNR in the LISA detector. The associated $P_ζ$ is accurately reconstructed as a flat spectrum with a maximum relative error of order 10%.

Our numerical pipeline can also be applied to other scenarios for the thermal history of the early Universe. Of particular interest would be the study of time-dependent EoS parameters on the SIGWs, like in the case of smooth-crossovers [411, 412], analogous to the QCD phase transition.

## 6.5 Non-Gaussian effects on SIGWs

As discussed in section 4.5, the tensor power spectrum of SIGWs receives contributions from the four-point correlation function of curvature perturbations. This contribution can be split into disconnected and connected terms. While the disconnected one depends only on the scalar power spectrum, the connected part arises from the primordial trispectrum, i.e. it is sensitive to primordial NG. And, as stressed in section 4.5, for local NG, $τ_{NL}$ would be the key observable to extract a constraint on NG from SIGWs. However, when the curvature perturbation originates from a single fluctuating degree of freedom beyond the inflaton, the parameters $τ_{NL}$ and $f_{NL}$ are connected as measures of higher-order correlations in the curvature perturbation, satisfying the relation $τ_{NL} = (\frac{6}{5}f_{NL})^2$, which saturates the Suyama-Yamaguchi inequality [413]. This relation has relevant implications for SIGWs, generated by models characterized by local-type NG, since observational constraints on one parameter can indirectly provide bounds on the other, assuming a given model. In the following analysis, we adopt the strategy of performing the analysis considering $f_{NL}$ as a parameter of the model and assume the shape (4.38) of the full power spectrum. We then discuss the implications for $τ_{NL}$ that derive from the constraints on $f_{NL}$. In this case, we assume the curvature power spectrum to have a LN profile, (see eq. (3.5))<fnref n="19" /> and as free parameters we use

$$ \theta_{\text{cosmo}} = \{\log_{10} A_s, \log_{10} \Delta, \log_{10} (k_*/s^{-1}), f_{\text{NL}} \equiv 5/6 \sqrt{\tau_{\text{NL}}}\}. \quad (6.10) $$

The left panel of figure 23 shows the absolute uncertainty associated with $f_{NL}$ varying the parameter $f_{NL}$ against $\log_{10}(k_*)$ computed using the FIM method. We fix the amplitude of $P_ζ$ to $\log_{10} A_s = -2$ and the width to $\log_{10} \Delta = -0.75$. The dashed black (purple) contour lines represent the relative percentage error associated with $f_{NL}$ when astrophysical foregrounds are not included (or are included). The white vertical line indicates the SNR. Notice that for small and large values of $k_*$, $f_{NL}$ exhibits higher uncertainties, whereas the intermediate range of $k_* \sim 10^{-2}-10^{-1} s^{-1}$ shows the minimal uncertainties for $f_{NL}$. This suggests that tight constraints on $f_{NL}$ can only be achieved within this specific range of $k_*$. Outside of this range, the errors increase notably, indicating less reliable measurements for $f_{NL}$. Even in the optimal case, the reconstruction of $f_{NL}$ only reaches the percent level for large $f_{NL} \gtrsim 12.5$. Notice

<footnote n="19">¹⁹Note that our choice of the LN $P_ζ$ assumes that the dimensionless primordial curvature fluctuations — including the higher order term coming from the trispectrum (the left-hand side of eq. (4.38) multiplied by $k^3/2\pi^2$) — describe a lognormal. This practically isolates the effect of computing $P_ζ$ including non-Gaussian contributions *ab-initio* from the effect caused by a NG contribution on the computation of $\Omega_{GW}$ given $P_ζ$. We will also compare our results to the LN case where only the Gaussian contribution is considered.</footnote>
<figure><img src="image_25.png" /><figcaption>**Figure 23.** Left panel: absolute uncertainties of each of the parameters in the case of an injected signal which includes NG and assuming the LN template with parameters fixed to $log_{10} A_s = -2$ and $log_{10} \Delta = -0.75$. Right panel: same as the left panel, but varying $\Delta$ instead of $k_*$, for an injected signal which includes NG and assuming the LN template with parameters fixed to $log_{10} A_s = -2$ and $log_{10}(k_*/s^{-1}) = -1.4$.</figcaption></figure>

that, in the presence of foregrounds, the accuracy on $f_{\text{NL}}$ slightly degrades, in particular when $k_*$ coincides with the expected peak of the galactic foreground, i.e. when $k_* \sim 10^{-2.2} s^{-1}$.

Similarly, the right panel of figure 23 shows the FIM absolute uncertainty associated with $f_{\text{NL}}$ varying the parameter $f_{\text{NL}}$ against $log_{10} \Delta$. We fix the amplitude of $P_{\zeta}$ to $log_{10} A_s = -2$ and the peak scale to the optimal location $log_{10}(k_*/s^{-1}) = -1.4$. In this case, the uncertainty when estimating $f_{\text{NL}}$ is lower ($\lesssim 30\%$) in the region of $log_{10} \Delta$ below -0.3, while it significantly degrades for larger widths. This suggests that the most stringent measurements of $f_{\text{NL}}$ will be obtained for relatively narrow curvature power spectra.

Given the relation between $f_{\text{NL}}$ and $\tau_{\text{NL}}$, improvements in the precision of the former directly translate into tighter constraints on the latter. The FIM analysis shows that percent-level accuracy on $f_{\text{NL}}$ is achievable only for large $f_{\text{NL}}$, which in turn would correspond to a percent-level constraint on $\tau_{\text{NL}}$. In favorable scenarios — where the peak scale $k_*$ lies within $10^{-2}-10^{-1} s^{-1}$ and the spectral width $log_{10} \Delta$ is relatively narrow — uncertainties in $f_{\text{NL}}$ are minimal, restricting the allowed range of $\tau_{\text{NL}}$. Conversely, when $f_{\text{NL}}$ is less precisely determined, $\tau_{\text{NL}}$ remains poorly constrained.

The corner plot in figure 24 illustrates the posterior distributions of the SIGW signal parameters {$log_{10} A_s$, $log_{10}(k_*/s^{-1})$, $log_{10} \Delta$}, including the primordial NG parameter $f_{\text{NL}}$, alongside the extragalactic and galactic background amplitude energy densities $log_{10}(h^2\Omega_{\text{Ext}}) \equiv \alpha_{\text{Ext}}$, $log_{10}(h^2\Omega_{\text{Gal}}) \equiv \alpha_{\text{Gal}}$. As in the other cases, we omit the LISA noise parameters $A_{\text{noise}}$, $P_{\text{noise}}$ from the corner plot as they are tightly constrained and weakly correlated with the rest. As first benchmark, we injected a template with $log_{10} A_s = -2$, $log_{10} \Delta = -0.75$, $log_{10}(k_*/s^{-1}) = -2.3$ and $f_{\text{NL}} = 1$. On the right panel of figure 24, we report the corresponding reconstructed $\Omega_{\text{GW}}h^2$. For comparison, we also plot the GW energy density in the Gaussian case. Due to the low value of $f_{\text{NL}}$ chosen, the reconstructed and Gaussian curves are almost superimposed, showing that the effects of NG are quite mild in this case.
<figure><img src="image_26.png" /><figcaption>**Figure 24.** Posterior distribution for an injected signal which includes NG, $f_{\text{NL}} = 1$ and assuming the LN with $\log_{10} A_s = -2$, $\log_{10} \Delta = -0.75$, $\log_{10}(k_*/s^{-1}) = -2.3$.</figcaption></figure>

The marginalized posterior for $f_{\text{NL}}$ exhibits a broad posterior distribution, ranging from -5 to 5. This distribution indicates a large uncertainty in $f_{\text{NL}}$ spanning roughly ±3 at the 68% C.L. and is compatible with $f_{\text{NL}} = 0$. This suggests a limited constraining power of LISA on $f_{\text{NL}}$. For such a signal, the posterior shows a bimodal structure, that arises because $f_{\text{NL}}$ enters quadratically in the GW spectral energy density through the trispectrum. As the observed value suggests compatibility with a Gaussian primordial distribution, the broad posterior distribution also implies that ruling out moderate NG will be challenging, emphasizing the need for improved precision or additional data to refine these estimates. The other parameters, such as the log-amplitude of the seed power spectrum ($\log_{10} A_s$), show a remarkable reconstruction, in line with the results of section 6.2.1. Note however that we are injecting different values of $\{\log_{10} A_s, \log_{10}(k_*/s^{-1}), \log_{10} \Delta\}$. For a comparison to the fully Gaussian case with the same injection in $P_\zeta$ see figure 31. The effect of adding the small NG correction $f_{\text{NL}} = 1$ on the recoverability of power spectral parameters $\{\log_{10} A_s, \log_{10}(k_*/s^{-1}), \log_{10} \Delta\}$ and foreground parameters $\{\alpha_{\text{Gal}}, \alpha_{\text{Ext}}\}$ is small in this case, indicating that a small NG contribution does not spoil the reconstruction of curvature power spectra parameters. This also indirectly supports our choice of not including NG corrections in the benchmark USR scenario discussed in section 3.3, which is characterized by $f_{\text{NL}} \simeq 0.09$. However, notice the visible non-zero correlation between $f_{\text{NL}}$ and both $\{\log_{10} A_s, \alpha_{\text{Gal}}\}$. The joint posterior $f_{\text{NL}}-\log_{10} A_s$ reflects the multimodality induced by the double peak structure of $f_{\text{NL}}$. Nevertheless, LISA can still strongly constrain the amplitude of the non-linear power spectrum. Hence the detection of the GWB is not strongly influenced by the primordial NG, which is beneficial for simplifying the analysis when considering models predicting small NGs.
<figure><img src="image_27.png" /><figcaption>**Figure 25.** Same as figure 24 with $f_{\text{NL}} = 10$. Note that there is a second, identical mode at $f_{\text{NL}} = -10$ as only the square of $f_{\text{NL}}$ enters $\Omega_{\text{GW}}$. We omitted this mode.</figcaption></figure>

The posterior distributions for astrophysical foreground amplitudes are broader and show limited correlation with signal parameters. In contrast, instrumental noise parameters are tightly constrained and largely independent of signal estimation. This decoupling ensures robust estimation of the signal amplitude, enhancing the reliability of SGWB detection despite noise uncertainties.

Overall, the joint distributions show only weak correlations, indicating that each of these parameters can be inferred with a good degree of independence from the others.

The situation changes when considering a larger value of $f_{\text{NL}}$, as shown in figure 25. Specifically, we run the MCMC with the following injected template with $\log_{10} A_s = -2$, $\log_{10} \Delta = -0.75$, $\log_{10}(k_*/s^{-1}) = -2.3$ and $f_{\text{NL}} = 10$, differing from the previous benchmark only in the choice of $f_{\text{NL}}$. On the right panel, we report the reconstructed $\Omega_{\text{GW}}h^2$ as well as the Gaussian counterpart for comparison. Given the higher value of $f_{\text{NL}}$, the differences are now more evident, resulting in an enhancement of the UV tail, but also in a slightly higher peak and more smoothed minimum (not visible in figure 25). Now, the marginalized posterior distribution for $f_{\text{NL}}$ shows only a narrow spread indicating that the estimate of $f_{\text{NL}}$ is precise. The absolute value of $f_{\text{NL}}$ has a mean reconstructed value that varies about $\pm 0.8$ at the 68% C.L. which makes it incompatible with $f_{\text{NL}} = 0$.

As with the previous case, the non-linear power spectral parameters are tightly constrained and no significant bias is observed with respect to the injected values. Both $k_*$ and $\Delta$ are very weakly correlated with $f_{\text{NL}}$, suggesting that their reconstruction is not heavily influenced when $f_{\text{NL}}$ is large. However, $\log_{10} A_s$ shows a larger anti-correlation with the $f_{\text{NL}}$. An anti-correlation also appears between $f_{\text{NL}}$ and $\alpha_{\text{Gal}}$, probably induced by the similar IR shape behavior. The tighter constraints associated with such a signal imply that higher levels of NG are easier to constrain, yielding clearer and more reliable effects.
The independence between $f_{\text{NL}}$ and other parameters (except for $\log_{10} A_s$), indicated by the weak correlations, suggests robustness in parameter estimation. This means that uncertainties in $f_{\text{NL}}$ do not drastically affect the inference of other parameters, resulting in more precise parameter constraints compared to the $f_{\text{NL}} = 1$ case. When $f_{\text{NL}}$ is larger, the signal is more distinct, allowing for setting more stringent constraints on primordial NG. Note that the prior in figure 25 is restricted to positive values of $f_{\text{NL}}$. Similarly to figure 24 the posterior distribution has a second mode at $f_{\text{NL}} = -10$ since it only enters quadratically in the signal.

For the amplitude $A_s$ considered in this case, the imprints due to possible inaccuracies in accounting for the full non-Gaussian behavior for some models of inflation are expected to be negligible when $f_{\text{NL}} = 1$, but could be substantial when $f_{\text{NL}} = 10$, as recently argued by [69]. For this analysis, we neglected those refinements. We further stress that for the enhanced amplitude of the power spectrum considered here, $f_{\text{NL}} = 10$ represents much larger deviations from Gaussianity than on CMB scales, because the expansion parameter determining the relative size of the trispectrum versus power spectrum is $\tau_{\text{NL}} \cdot P_\zeta$. It is of order one in the current context, while less than $10^{-4}$ on CMB scales.

Concerning the implications for $\tau_{\text{NL}}$, large uncertainties in $f_{\text{NL}}$ directly translate into poor constraints on $\tau_{\text{NL}}$. For example, for $f_{\text{NL}} = 1$ with a quite broad uncertainty of $\pm 3$, $\tau_{\text{NL}}$ could span from values close to zero (if $f_{\text{NL}} \approx 0$) up to $\approx 23$ (if $f_{\text{NL}} \approx 4$), making it challenging to clearly identify a primordial NG signal. In this range, even moderate NGs become difficult to distinguish from a Gaussian spectrum. Without improved precision on $f_{\text{NL}}$, the corresponding $\tau_{\text{NL}}$ will remain poorly determined, limiting our ability to discriminate between different levels of primordial NG. When $f_{\text{NL}} = 10$, providing a more pronounced non-Gaussian signal, the corresponding $\tau_{\text{NL}} = (\frac{6}{5} \cdot 10)^2 = 144$ is now much more tightly constrained. Since the uncertainty in $f_{\text{NL}}$ is roughly $\pm 0.8$, $\tau_{\text{NL}}$ varies up to a $\pm 15\%$ range. This tighter range is obviously better than the scenario with small $f_{\text{NL}}$. Hence, larger $f_{\text{NL}}$ values significantly improve our ability to determine $\tau_{\text{NL}}$, allowing LISA to better distinguish between different levels of primordial NG in the SGWB.

Finally, it is important to highlight that while Planck provides constraints that are very close to zero [74], indicating no evidence for primordial NG at large scales, the analysis we are performing for LISA focuses on NG at much smaller scales. LISA's ability to provide tight constraints on $f_{\text{NL}}$ suggests that GW detection could play a crucial role in refining our understanding of primordial NG, particularly in scenarios where the signal is expected to be strong. In addition, the sensitivity of LISA to different scales compared to Planck provides an important cross-check, helping to verify any scale dependence for $f_{\text{NL}}$ [414, 415]. Overall, while Planck remains a benchmark for CMB-based constraints on $f_{\text{NL}}$ at large scales, LISA shows the potential of GW detectors to significantly advance the search for, and the characterization of primordial NG.

# 7 Testing the scalar-induced hypothesis

In this section, we outline a procedure to test the compatibility of the SIGW hypothesis with a possible SGWB detection. So far, our analyses have assumed that the cosmological contribution of the SGWB originates from SIGWs. There are, however, many alternative
possible physical mechanisms for sourcing SGWBs in the early universe. Our goal is to
offer a practical approach for assessing the validity of the hypothesis explored in this work
— namely, whether or not a hypothetically detected signal originates from enhanced scalar
fluctuations of inflationary origin. To this end, we focus on two illustrative scenarios that are
distinct in nature, leaving a detailed comparison of various early-universe signals — which
is beyond the scope of this paper — for future work. We use the evidence (5.22) as an
estimator for model selection. Specifically, given an injected signal, we consider different
reconstruction techniques, for which we can compute the (log) evidence using the nested
sampler PolyChord [416, 417]. Then, we compare different hypotheses by computing the
Bayes factors to test which one is favoured given the observed signal. It is important to note
that the Bayes factor is a global estimator: it not only assesses the goodness of fit of a model
to the data but also incorporates information about the prior volume and its compression
as the prior transitions to the posterior. Therefore, given a similar fit to data, it naturally
favours simpler models — those with fewer parameters — over more complex ones.

In the following analysis, we compare two different hypothesis:

**Hypothesis 1.** The SGWB can be modelled directly in $\Omega_{\text{GW}}(f)$, independently of whether or not it is scalar-induced. Thanks to the SGWBinner code [63, 64], we can use both *i*) a template-based [6], as well as *ii*) a model-agnostic approach to model $\Omega_{\text{GW}}$ when reconstructing the signal. Unlike the model-agnostic methods presented in this work, SGWBinner's methods do not assume the underlying physics. The templates we use to fit the signal are of the same functional form as the injection template. Of course, with real LISA data, we will not necessarily know the correct template, and must either rely on model selection between several physically-motivated templates or, alternatively, model-agnostic methods. The binned approach of SGWBinner divides the frequency space into bins by SNR and then fits a power law within each bin. This results in an agnostic reconstruction of $\Omega_{\text{GW}}$.

**Hypothesis 2.** The SGWB is scalar-induced, and can be modelled at the level of $P_\zeta(k)$. By contrast, the various techniques presented in the previous sections assume that the signal that we are considering is due to SIGWs. We can reconstruct the SGWB exploiting the SIGWAY code in two different ways. One option is to use *i*) a template-based approach (see section 3.2), in which a template for $P_\zeta$ is specified. A second possibility is to use *ii*) the $P_\zeta$-agnostic (still assuming SIGW to be the source of the SGWB) binned approach, as described in section 3.1.<fnref n="20" /> We choose the number of bins $N_{\text{bins}} = 40$. We assume the SIGW to be produced during the radiation-dominated era.

<footnote n="20">²⁰Let us note that there is a difference in the implementation of the binned $P_\zeta$ approach, compared to the binned $\Omega_{\text{GW}}$. In the latter approach, the SGWBinner code dynamically selects the optimal number of bins before the nested sampling, based on the Akaike information criterion [418] (we refer the interested reader to the discussion around eq. (3.6) of [63] for more details). On the other hand, such a feature has not been implemented in the SIGWAY code, as there are fundamental difficulties to attempting a similar approach in $P_\zeta$-space (see appendix B for a detailed discussion) so the number of bins for $P_\zeta(k)$ has to be chosen by hand.</footnote>
<figure><img src="image_28.png" /><figcaption>**Figure 26.** Reconstruction of $\Omega_{\text{GW}}$ for case 1 (left) and case 2 (right). We only show the reconstruction of the injected cosmological contribution to the SGWB. It is clearly visible that the two models that assume SIGWs cannot reconstruct the narrow log-normal peak in case 1 and therefore lead to much lower evidence.</figcaption></figure>

We simulate the following two qualitatively different signals:

**Case 1.** The injection is not SIGWs. The first signal is a narrow lognormal in $\Omega_{\text{GW}}(f)$. This injection serves as an example of a signal that cannot be produced by SIGWs, assuming the modes reenter during RD. In this case, regardless of how narrow the peak in $P_\zeta$ is, the generated SIGW will always exhibit the so-called “causality tail” proportional to $f^3$. As a benchmark for this injection, we chose to reproduce the main peak of the double-peak background shown in the top-left panel of figure 11 in [6], with a slightly lowered amplitude. This amounts to choosing the following parameters: $\log_{10}(h^2\Omega_*) = -9.5$, $\log_{10}(f_*/\text{Hz}) = -2.21$, $\log_{10}(\rho) = -1.10$ in eq. (2.8) of [6]. We will henceforth call this signal $\Omega_{\text{GW}}^{\text{LN}}$.

**Case 2.** The injection is SIGWs. The second injection is instead derived from a SIGW scenario. We inject a lognormal power spectrum of curvature perturbations, see eq. (3.5) and compute the resulting SGWB numerically. The power spectrum parameters used were: $\log_{10} A_s = -2.3$, $\log_{10} \Delta = -0.70$, $\log_{10}(k_*/\text{s}^{-1}) = -1.5$. The resulting shape in the SGWB can be described through the double-peak template in eq. (2.10) of [6], with parameters $\{\log_{10}(h^2\Omega_*), \log_{10}(f_*/\text{Hz}), \beta, \kappa_1, \kappa_2, \rho, \gamma\} = \{-9.5, -5, 0.242, 0.456, 1.234, 0.08, 6.91\}$. With this choice, the main peak of the SIGW coincides with the injection of Case 1.

We fit each of the two injections using the four models specified above. The results of our analysis are summarized in figure 26 (see also figure 32 in appendix D showing the reconstruction including foregrounds), and table 1, which we now describe in order.
In Case 1, the models that perform the worst in terms of model selection are the two based on the SIGW hypothesis (marked by the row SIGWAY). As shown in the left panels of figure 26, both the lognormal and binned $P_ζ(k)$ models attempt to fit the lognormal SGWB using the primary peak of the SIGW. However, the secondary peak at lower frequencies severely undermines the fit to the data, resulting in a very poor evidence value. According to Jeffreys' scale (see section 5.2), these models are decisively ruled out when compared to the two alternative hypotheses. As expected, the binned $P_ζ(k)$ performs worst, due to its significantly larger number of parameters.

We also fit the same injection using a lognormal template for $\Omega_{GW}$ (which corresponds to the true injection) and a binned $\Omega_{GW}$ reconstruction with SGWBinner. As illustrated in figure 26, both models reconstruct the signal very well, closely matching the injection. Furthermore, they yield very similar best-fit likelihood values with a preference for the lognormal template (matching the injected signal) due to its smaller number of parameters.

The main takeaway from Case 1 is that if a similar signal was detected, we could conclude with very high statistical significance that the signal does not have a scalar-induced origin.

In Case 2, all the models considered successfully recover the injected signal. Unlike in the previous case, the Bayes factors are closer together, while still being orders of magnitude apart (we would like to stress that the differences quoted in table 1 are $\ln(Z)/10^4$). Again, the model-agnostic reconstructions — whether in $\Omega_{GW}$ or $P_ζ(k)$ — perform the worst. However, their flexibility and agnostic nature still makes them useful in the real data analysis to pre-constrain the shape of $\Omega_{GW}$ and $P_ζ$. The Bayes factors are shown in table 1.

The accuracy of the recovery of the injected signal is sufficient to show a strong preference for the injected SIGWAY template in $P_ζ$ when compared to the empirical template for $\Omega_{GW}$ developed in [6, eq. (2.10)]. This outcome arises because the template in [6] is described by seven free parameters compared to three for the lognormal $P_ζ(k)$, and because there is a slight mismatch between the numerical computation and the phenomenological fit. This highlights the need for accurate models for the SGWB. This behaviour can be seen in the right column of figure 26.

On the other hand, when comparing the agnostic approaches, we see that the former gives slightly better $\Omega_{GW}$ reconstructions around the peak of the signal, with an oscillatory behaviour of the residuals due to the finite resolution of the binned approach (40 bins). This result stresses the importance of adopting optimal modelling of the eventual cosmological scenario when reconstructing the signal (see also appendix D).

All in all, the results of this section, although based on two illustrative examples, demonstrate that LISA has the potential to confirm or rule-out the scalar-induced nature of the SGWB with high statistical significance. These results confirm that the true models achieve the best Bayes factors when appropriately matched. Alternative and binned models consistently show inferior fits, highlighting the distinctiveness of the injected signals.

However, is important to stress that not all SGWB that may appear in LISA lead to such a clear difference between signals that can or cannot be generated by SIGW. The characteristic double-peak structure that we observe with the injected $P_ζ^{LN}$ signal is only measurable by LISA if (a) the peak in $P_ζ$ is sufficiently narrow, and (b) both peaks in $\Omega_{GW}$ happen to fall within the sensitivity of LISA. On the other hand, there are many potential shapes for $P_ζ$ where these conditions are not met. In these cases the SIGW signal can easily
<table><caption>**Table 1.** In Bayes factors (normalized to a reference value $10^4$) comparing the SGWB reconstruction (either using the injected model as a template, or a model-agnostic binned method) with the SIGW reconstruction for two signals: i) a log-normal power spectrum in $\Omega_{GW}$ which cannot be generated by SIGW (within the assumptions we are working in), ii) a log-normal power spectrum in $P_\zeta$ which generates detectable SIGW. In bold, we show the Bayes factors for the recovery which assumes the injected template. As expected, they are the best reconstructions for each injection.</caption><thead><tr><th>Package</th><th>Method</th><th>Case 1: not SIGW signal<br>ln [Z(Ω<sub>GW</sub><sup>LN</sup>)]/10<sup>4</sup></th><th>Case 2: SIGW signal<br>ln [Z(P<sub>ζ</sub><sup>LN</sup>)]/10<sup>4</sup></th></tr></thead><tbody><tr><td rowspan="2">SGWBinner</td><td>Template</td><td><b>0.5444</b></td><td>-0.5278</td></tr><tr><td>Binned</td><td>-0.5625</td><td>-0.5479</td></tr><tr><td rowspan="2">SIGWAY</td><td>Template</td><td>-25.6078</td><td><b>-0.5203</b></td></tr><tr><td>Binned</td><td>-25.7934</td><td>-2.712</td></tr></tbody></table>

mimic one expected from other cosmological sources, potentially making it much harder to rule out models. We will leave a more detailed discussion of this for future work.

## 8 Conclusions

In this work, we investigated the potential of the LISA detector for reconstructing the SGWB sourced by second-order scalar perturbations. Three approaches were explored: a binned spectrum reconstruction, template-based methods, and a direct modeling approach rooted in first-principles scenarios (taken to be the single-field USR inflationary model for presentation purposes).

Our results demonstrate that the direct modeling approach yields the tightest constraints on the primordial curvature power spectrum $P_\zeta$, particularly capturing both the IR and UV tails of the signal with better precision than alternative methods, due to the stronger prior information inevitably included in the fit. This highlights the power of incorporating *ab initio* physics into signal reconstruction pipelines to leverage the constraining power of LISA observations at their best.

The binned spectrum reconstruction approach is complementary and proved effective in providing model-independent upper bounds on $P_\zeta$ when the cosmological contribution to the SGWB is below the sensitivity. Capturing the overall shape of the SIGW spectrum with this approach proved to be difficult, due to a combination of missing SNR towards the edges of the LISA window and strong degeneracies between the bins when choosing a large number of them. Despite these shortcomings, it is possible — if a cosmological contribution to the SGWB is detected — to tell apart signals that can be SIGW from those that cannot, by Bayesian model selection.

In comparison, the template-based methods provided a more consistent reconstruction across frequencies, though they inherently rely on prior assumptions about the shape of the spectrum. The complementary strengths of these approaches suggest that an optimal reconstruction strategy would involve their combined use, with model-dependent templates guiding reconstructions and binned methods offering flexibility in capturing unanticipated spectral features or in setting bounds that are agnostic on the spectral shape.
We also examined the impact of going beyond the simplest vanilla cosmological scenarios on the SIGW reconstruction investigating the sensitivity of SIGW signals to early-universe physics. In particular, we included in our analysis the study of the effect of the transition from early matter-dominated to radiation-dominated eras, as well as the role of non-Gaussianity in the SIGW spectrum. Future research will include the study of the effect of a time-dependent equation-of-state parameter on the SIGW spectrum, as would be generated in the case of a smooth crossover [411, 412], such as in the QCD phase transition. Overall our analysis demonstrated how SIGW searches in LISA will provide constraints that vastly outperform those deduced from the effective number of relativistic species $\Delta N_{\text{eff}}$ and PBH overproduction bounds. In this regard, SIGW searches will also be an invaluable tool for probing the asteroid mass window of PBH dark matter.

Looking forward, several key avenues remain open for future work. On the phenomenological side, it remains an open question how well LISA will be able to constrain primordial non-Gaussianity or non-standard thermal histories while allowing for a fully non-parametric curvature power spectrum. In this work we have only quoted template-based constraints on these effects, e.g. in figures 22 and 24, which do not consider possible degeneracies with the shape of the spectrum.

Concerning the binned approach to reconstructing $P_ζ$, a more mature method that incorporates assumptions about the smoothness of the spectrum of scalar perturbations and addresses the computational difficulties with the binning will be needed in the future. It is likely that — even allowing for non-Gaussianities and non-standard thermal histories — some shapes of the SGWB cannot be scalar-induced and can be confirmed as signatures of directly sourced tensor perturbations, even without identifying a specific early-universe source.

On the theory side, future work should consider expanding the scope of single-field scenarios by exploring more general actions in the Jordan frame, incorporating non-minimal coupling and non-canonical kinetic terms (e.g. eq. 2.1 of [87]). Also, going beyond single-field USR models, first-principle multi-field inflationary scenarios merit investigation, as discussed in section 2.2. In some cases, multi-field models can effectively be reduced to single field descriptions, making some of the techniques developed here already applicable, and enabling simpler parameter space scans as done in [419, 420]. A reverse engineering approach could be particularly valuable, where inflationary dynamics are modeled based on a minimal set of parameters, and the corresponding inflationary potential is reconstructed within single- or multi-field frameworks, as demonstrated in refs. [276, 421]. Expanding the framework for computing the non-Gaussian signatures predicted in most SIGW models beyond the lognormal template in section 4.5 could serve as a diagnostic tool for breaking degeneracies in cases where the SGWB spectrum alone is insufficient. Also, implementing the binned approach in scenarios with NGs could also allow us to reduce the computational costs of these analyses. While in this work we only considered the monopole signal, as non-Gaussianities may impact the large scale SIGW anisotropies, it would be interesting to include information from higher order in the multiple expansion of power in the sky [326, 422–424]. Finally, integrating these advanced modeling and reconstruction techniques into the global fit pipeline of LISA, as well as incorporating measurements from other experiments, will be essential for unlocking the mission’s full potential in probing the early universe’s cosmological landscape.
Acknowledgments

We thank Nicola Bartolo, Gaetano Luciano, Marco Merchand, Sabino Matarrese, and Toni Riotto for discussions and interactions in an early stage of this project. We acknowledge the LISA Cosmology Working Group members for seminal discussions. We especially thank the authors of refs. [6, 338, 339, 342] for the collaborative developments of the SGWBinner code upon which we built to produce many of the forecasts we presented in this work. The research of JF is supported by the grant PID2022-136224NB-C22, funded by MCIN/AEI/ 10.13039/501100011033/FEDER, UE, and by the grant/2021-SGR00872. RR was supported by an appointment to the NASA Postdoctoral Program at the NASA Marshall Space Flight Center, administered by Oak Ridge Associated Universities under contract with NASA. The work of EM is supported by the Italian Ministry of University and Research grant Rita Levi-Montalcini "New directions in axion cosmology". EM acknowledges the support of Istituto Nazionale di Fisica Nucleare (INFN) through the iniziativa specifica TAsP. JE and GN acknowledge support by the ROMFORSK grant project no. 302640. TP acknowledges the contribution of the COST Action CA21136 "Addressing observational tensions in cosmology with systematics and fundamental physics (CosmoVerse)" and of the INFN Sezione di Napoli iniziativa specifica QGSKY. He acknowledges as well financial support from the Foundation for Education and European Culture in Greece. DR is supported by the UZH Postdoc Grant 2023 Nr. FK-23-130. During the course of this work, S.RP and DW were supported by the European Research Council under the European Union's Horizon 2020 research and innovation programme (grant agreement No 758792, Starting Grant project GEODESI). The work of AG is supported by the UKRI AIMLAC CDT, funded by grant EP/S023992/1. The work of AG, GT, and IZ is partially funded by the STFC grants ST/T000813/1 and ST/X000648/1. The work of HV was supported by the Estonian Research Council grants PSG869 and RVTT7 and the Center of Excellence program TK202. M.Pe acknowledges support from Istituto Nazionale di Fisica Nucleare (INFN) through the Theoretical Astroparticle Physics (TAsP) project and from the MIUR Progetti di Ricerca di Rilevante Interesse Nazionale (PRIN) Bando 2022 - grant 20228RMX4A, funded by the European Union - Next generation EU, Mission 4, Component 1, CUP C53D23000940006. JK acknowledges support from the JSPS Overseas Research Fellowships and the INFN TAsP project. GP acknowledges partial financial support by ASI Grant No. 2016-24-H.0. AR acknowledges support from Istituto Nazionale di Fisica Nucleare (INFN) through the iniziativa specifica TEONGRAV and by the project BIGA - "Boosting Inference for Gravitational-wave Astrophysics" funded by the MUR Progetti di Ricerca di Rilevante Interesse Nazionale (PRIN) Bando 2022 - grant 20228TLHPE - CUP I53D23000630006.

**Authors' contribution.** JE: Main developer of SIGWAY. Implementation and optimization of techniques discussed in appendix A, coding the interface to SGWBinner. Co-coding different thermal histories. Running analyses and producing figures for all results except sections 6.5 and 7. Writing sections 6 and appendices. Reviewing all draft. AG: Co-coding different thermal histories. Writing sections 4.3 and 6.4. Reviewing section 3. GF: Proposing and coordinating the project with RR, including defining goals and managing tasks. Co-coding parts of SIGWAY, mainly on SIGW computations and the USR module. Writing and Reviewing
all sections of the draft. TP: Co-coding kernel functions and different thermal histories. Writing section 6.4 and parts of sections 2.1, 2.2, 2.3, 3.2 and 8. Reviewing sections 6.1, 6.2 and 6.3. MPe: Devising methods for computations in section 3.3 and developing the algorithm for the binned spectrum approach, contributing to the code implementation. Writing sections 3.1 and 4.4. Reviewing sections 1, 2, 5, and 8. GP: Co-coding SIGW computations with NGs in SIGWAY. Running and producing figures for section 6.5. Writing appendix A.4 and section 4. Reviewing section 4. MPi: Co-coding SIGW computation in the SIGWAY and interfacing with SGWBinner. Writing section 5. AR: Co-coding SIGW computations with NGs in SIGWAY. Running and producing figures for section 6.5. Writing section 6.5. Reviewing section 1. RR: Proposing and coordinating the project with GF, including defining goals and managing tasks. Co-coding the USR module. Writing and Reviewing all sections of the draft. GT: Developing the algorithm for the binned spectrum approach, contributing to the code implementation. Writing sections 3.1 and 4.4. Reviewing the draft. MB: Running the analysis and writing of section 7. JF: Running the analysis and writing of section 7. Reviewing sections 7, 3.2.2, 2.2 and section 2. JK: Coding SGWBinner. Writing section 5. Reviewing section 3. EM: Contributing to an early version of the USR code. Writing sections 2.1, 2.2, 2.4, 3.2, 3.3, 6.2, and 6.3. Reviewing all draft. GN: Writing section 1. Reviewing all draft. DR: Writing sections 3 (intro), 4.1, 4.2, and 6.2.1. Reviewing sections 4 and 7. SRP: Devising and interpreting section 4.5 and 6.5. Writing section 4.5 and related elements in sections 4 and 6.5. Reviewing all draft. HV: Devising sections 2 and 3. Writing sections 2, 3, and 6.2. Reviewing sections 1, 6.4, 6.5, and 8. DW: Devising and interpreting section 4.5. Reviewing sections 4.5, 1, and 8. IZ: Writing parts of section 4 and reviewing sections 2, 3, and 4.

# A SIGWAY code: technicalities

In this appendix, we describe some technical aspects of the SIGWAY code developed for the analysis performed in this work.

## A.1 Perturbations in USR scenarios: code structure

Given a potential $V(\phi)$ inducing a USR phase, the curvature power spectrum is computed in three steps that are described below. $P_{\zeta}(k)$ is then interpolated and $\Omega_{GW}h^2$ is computed as described in A.2.

Using the notation described in section 3.3, the inputs in the code are:

* the inflaton potential $V(\phi)$;

* the number of $e$-folds from when CMB modes exit the Hubble horizon to the end of inflation $N_{\text{CMB}\to\text{end}}$.

* The initial conditions $\phi_0$ and $\pi_0 = \phi'_0$.

The code then automatically defines the dimensionless variables in eq. (3.19). For definiteness, in this paper, we have fixed $N_{\text{CMB}\to\text{end}} = 58$. Fixing the number of $e$-folds from the CMB scale to the end of inflation effectively allows us to set the correspondence between $N$ and
wavenumbers $k$. Fixing $N_{\text{CMB}}$ also implicitly fixes the thermal history of reheating and subsequent phases. We do not model these eras for simplicity, but they would be fixed in a complete USR+reheating model. Also, we checked that the initial SR attractor would quickly pull the field space trajectory to the background evolution, and thus one could also assume negligible initial velocity for simplicity.

The computation then proceeds as follows:

1. Solve for the background evolution using eqs. (3.21). We evolve the dynamics until the first SR condition is violated ($\epsilon_H = 1$), collecting $x(N)$, $y(N)$ and $h(N)$. We denote the number of e-folds at the end of inflation $N_{\text{end}}$.

2. Compute the relation between the wave number $k$ and the number of e-folds $N$ at Hubble crossing using

$$k = k_{\text{CMB}} \frac{h(N)}{h_{\text{CMB}}} \exp(N - N_{\text{CMB}}), \qquad (A.1)$$

where $k_{\text{CMB}} = 0.05/\text{Mpc}$, $N_{\text{CMB}} = N_{\text{end}} - N_{\text{CMB}\to\text{end}}$ and $h_{\text{CMB}} = h(N_{\text{CMB}})$. With this, we can compute $P_{\zeta}(k)$ in the SR approximation according to eq. (3.25).

3. We are only interested in computing the spectrum of curvature perturbations beyond the SR approximation for modes of relevance for LISA. Once we have computed $k(N)$, a set of modes covering the LISA frequency band $\sim (10^{-5}-10^{-1})\text{Hz}$ is selected and those modes are evolved according to eq. (3.31) from $N_{\text{in}}$ to $N_{\text{out}}$, where $N_{\text{in}}$ and $N_{\text{out}}$ can be set by the user. We found $N_{\text{in}} = N-3$ and $N_{\text{out}} = N+7$, where $N$ is the horizon crossing time of the mode, to be long enough for the modes to freeze out as there is no super horizon evolution in these models once USR has ended. Notice that the USR can only last about $\mathcal{O}(3)$ e-folds without making perturbations grow beyond the validity of perturbation theory. Additionally, for simple shapes of the potential, we find that evolving $\sim 100$ modes in $k$ and interpolating between them yields sufficient precision.

4. Lastly, $P_{\zeta}(k)$ is computed without resorting to the SR approximation using eq. (3.32). This quantity is then passed to the algorithm computing $\Omega_{\text{GW}}h^2$ which can then be passed to the LISA likelihood.

Figure 27 shows the evolution of a number of modes in $P_{\zeta}$ as a function of $N$ across the USR phase. As this algorithm is called many times when sampling the LISA likelihood, we are using the package `diffrax` [425] to solve the inflationary perturbation equation of motion.

## A.2 Computation of SIGWs from the spectrum of curvature perturbations

For a given shape of $P_{\zeta}(k)$, computing $\Omega_{\text{GW}}h^2$ is relatively straightforward by evaluating the double integral in eq. (4.20). The main concern here is making this computation as fast as possible, as this dominates the time it takes for each likelihood evaluation. To take full advantage of vectorization, we compute $s, t, k$ on a grid, where $s$ is linearly spaced, and $t, k$ are logarithmically spaced. The integration in $t$ is convergent for any realistic shape of $P_{\zeta}$ as both the value of $\overline{I^2(k, s, t)}$ tends towards 0 for large $t$ and there needs to be some cutoff controlling the amplitude of $P_{\zeta}$ at large momenta not to violate BBN bounds. If
<figure><img src="image_29.png" /><figcaption>**Figure 27.** Power spectrum as a function of $N$ for different modes $k$ indicated by the color palette, close to the onset of the USR phase $N - N_{\text{end}} \simeq -28$. This plot shows the evolution and freezing out of modes for different values of $k$ for the potential from eq. (2.4). The evolution of each mode is traced for around $\Delta N \simeq 10$ e-folds from within the Hubble sphere to sufficiently after Hubble crossing and freezing. The time of Hubble crossing is marked with gray dots.</figcaption></figure>

$P_{\zeta}(k)$ features a scale after which it drops rapidly, it can be advantageous to define a custom $t(k)$ and compute the grid as $s, t(k), k$.

After evaluating the integrand in eq. (4.20) on the grid, the integral is computed with Simpson-integration. Depending on the shape of $P_{\zeta}$, we found a number of points in the grid around $N_s \sim 10-100$, $N_t \sim 300-1000$ and $N_k \sim 100$ would give sufficient precision.

To take full advantage of threading, we use JAX and just-in-time compilation [385] for computing the integrand and performing the integration itself. Altogether, this results in a significant speedup over other publicly available codes (e.g. [426]). We record wall-clock times of $\sim 10^{-2}$s for calculating 100 values of $\Omega_{\text{GW}}$ for $P_{\zeta}$ containing only JAX-native functions. In the case discussed in section 6.2.2 this time rises to about 1s. Crucially, this speedup allows us to sample the posterior distribution efficiently and obtain good MCMC convergence with a laptop on timescales of $\mathcal{O}(\text{hours})$.

## A.3 Computation of SIGWs using binned coefficients

In the case where we bin $P_{\zeta}(p)$ in momentum space, $\Omega_{\text{GW}}(k)$ can be computed in a straightforward manner through eq. (4.27). The most computationally expensive part of this is the computation of the coefficients $\Omega_{\text{GW}}^{(i,j)}(k)$. Luckily we can precompute these coefficients on a grid and save them as a $N_k \times N_p \times N_p$ tensor $K_{ij}^k \equiv \Omega_{\text{GW}}^{(i,j)}(k)$. At runtime, we then compute a $N_p \times N_p$ tensor of $B = A \otimes A$, with which we can conveniently re-write eq. (4.27) as

$$ \Omega_{\text{GW}}^k = \Omega_{ij}^k B^{ij}, \qquad (A.2) $$

where we used Einstein sum convention. The reason for doing this arguably very simple conversion is that the matrix equation is fully vectorizable with JAX.

Using this trick we are able to compute $\Omega_{\text{GW}}$ from $P_{\zeta}$ in $\lesssim 10^{-3}$s for $N_k = N_p = 50$, thus making inference possible despite a large number of parameters to sample.