aainstitutetext: Center for Cosmology and Particle Physics, New York
University, 726 Broadway, New York, NY 10003, USAbbinstitutetext: Instituto de
Estructura de la Materia, CSIC, Serrano 121, 28006 Madrid,
Spainccinstitutetext: Theoretical Physics Department, CERN, 1211 Geneva 23,
Switzerlandddinstitutetext: Departament de Física Quàntica i Astrofísica i
Institut de Ciències del Cosmos, Universitat de Barcelona, Martí i Franquès 1,
08028 Barcelona, Spainhhinstitutetext: Department of Mathematics and Physics,
University of Stavanger, NO-4036 Stavanger, Norwayffinstitutetext:
Dipartimento di Fisica e Astronomia “Galileo Galilei”, Università di Padova,
35131 Padova, Italygginstitutetext: INFN, Sezione di Padova, 35131 Padova,
Italyhhinstitutetext: Institut d’Astrophysique de Paris, UMR 7095 du CNRS et
de Sorbonne Université, 98 bis bd Arago, 75014 Paris, Franceiiinstitutetext:
Dipartimento di Fisica “Enrico Fermi”, Università di Pisa, Largo Bruno
Pontecorvo 3, Pisa I-56127, Italyjjinstitutetext: INFN, Sezione di Pisa, Largo
Bruno Pontecorvo 3, Pisa I-56127, Italyllinstitutetext: Physics Department,
Swansea University, SA28PP, United Kingdommminstitutetext: Dipartimento di
Fisica e Astronomia, Università di Bologna, and
INFN, Sezione di Bologna, I.S. FLAG, viale B. Pichat 6/2, 40127 Bologna,
Italynninstitutetext: Keemilise ja bioloogilise füüsika instituut, Rävala pst.
10, 10143 Tallinn, Estonia

# Gravitational waves from inflation in LISA: reconstruction pipeline and
physics interpretation

Matteo Braglia111Corresponding author<EMAIL_ADDRESS>b , Gianluca Calcagni
c , Gabriele Franciolini 2, d , Jacopo Fumagalli222Corresponding author:
<EMAIL_ADDRESS>3,e , Germano Nardini333Project coordinator:
<EMAIL_ADDRESS>4,f,g , Marco Peloso444Project coordinator:
<EMAIL_ADDRESS>5, c , Mauro Pieroni555Corresponding author:
<EMAIL_ADDRESS>h , Sébastien Renaux-Petel i,j,f , Angelo
Ricciardone l,m , Gianmassimo Tasinato f,g,n , Ville Vaskonen
(For the LISA Cosmology Working Group)

###### Abstract

Various scenarios of cosmic inflation enhance the amplitude of the stochastic
gravitational wave background (SGWB) at frequencies detectable by the LISA
detector. We develop tools for a template-based analysis of the SGWB and
introduce a template databank to describe well-motivated signals from
inflation, prototype their template-based searches, and forecast their
reconstruction with LISA. Specifically, we classify seven templates based on
their signal frequency shape, and we identify representative fundamental
physics models leading to them. By running a template-based analysis, we
forecast the accuracy with which LISA can reconstruct the template parameters
of representative benchmark signals, with and without galactic and
extragalactic foregrounds. We identify the parameter regions that can be
probed by LISA within each template. Finally, we investigate how our signal
reconstructions shed light on fundamental physics models of inflation: we
discuss their impact for measurements of _e.g.,_ the couplings of inflationary
axions to gauge fields; the graviton mass during inflation; the fluctuation
seeds of primordial black holes; the consequences of excited states during
inflation, and the presence of small-scale spectral features.

††preprint: CERN-TH-2024-072

##  1 Introduction

The Laser Interferometer Space Antenna (LISA) LISA:2017pwj  is a planned
space-based gravitational wave (GW) detector designed for the mHz frequency
band. This mission, led by the European Space Agency (ESA), involves ESA
member countries and substantial contributions from NASA and other
international space agencies. The mission was adopted in 2024. During this
phase, established Figures of Merit bind the LISA sensitivity and ensure its
main mission objectives Colpi:2024xhw . With fewer variables affecting mission
performance, the scientific community can design and adjust the data analysis
pipelines of signal searches and science interpretations toward the
established experimental configuration. The core of the adjustments should be
finalized over the next six years or so. This in fact provides ample time for
ESA to validate and properly code the analyses, before the planned launch in
the mid 2030s.

The search for the cosmological Stochastic Gravitational Wave Background
(SGWB) poses a major challenge for LISA due to the similarities between its
statistical properties and those of the instrumental noise and of unresolved
astrophysical events. Therefore, it will be challenging to extract and
characterise the contribution of these different components from the data
stream. This problem is further complicated by the fact that there is no
unique expectation for the cosmological SGWB, neither concerning its amplitude
nor its spectral (frequency) shape, as there are several, non mutually
exclusive, potential phenomena in the early universe that could generate it.
However, there is a vast amount of literature on possible cosmological GW
signals, and one can then envisage classifying and compiling them into a list
of signal predictions to search for.

If the primordial SGWB signal is not sufficiently strong to dominate over most
of the transient events (see _e.g.,_ ref. Braglia:2024siw  for a study of the
strong background case), its reconstruction at LISA is planned to be achieved
via a “global fit”, where all sources are simultaneously reconstructed in an
iterative manner Cornish:2005qw ; Vallisneri:2008ye ;
MockLISADataChallengeTaskForce:2009wir ; Littenberg:2023xpl . Based on current
estimates, the global fit is too computationally expensive to be repeated for
the whole list of primordial SGWB templates Littenberg:2023xpl . A more
feasible approach is to consider a global fit where the primordial SGWB is
firstly isolated with an approximate template-free approach (see _e.g.,_ refs.
Caprini:2019pxz ; Flauger:2020qyi ; Karnesis:2019mph ), and then use this
reconstruction to shortlist the theoretically-motivated templates (_i.e.,_
parameterisation of the expected SGWB power spectra) that best suit the
reconstructed signal. When several theoretically well-motivated templates
exhibit minor differences from each other, it is more efficient to group them
into approximate templates, and to introduce an intermediate step between the
template-free analysis and the refined theory-based template analysis. At this
point, pursuing the global fit with a few branches of the cosmological SGWB
templates becomes affordable. In light of the above reasoning, a comprehensive
list of cosmological SGWB templates is required.

The present work is part of a series of three papers that have the purpose to
start the collection of this ‘template SGWB databank’. While in the works
Caprini:2024hue ; Blanco-Pillado:2024aca , templates motivated by,
respectively, phase transitions and topological defects are considered, here
we propose and study templates that are motivated by broad classes of SGWB
production mechanisms associated with primordial inflation. Inflation
Guth:1980zm ; Starobinsky:1980te ; Linde:1981mu ; Albrecht:1982wi  solves
several shortcomings in cosmology, and can naturally produce primordial
fluctuations with properties in agreement with observations Mukhanov:1981xt ;
Hawking:1982cz ; Starobinsky:1982ee ; Guth:1982ec ; Bardeen:1983qw . However
the precise inflationary model, and even the energy scale at which inflation
took place, are still unknown. The most minimal models of inflation predict a
slightly red-tilted SGWB spectrum, and that, once the cosmic microwave
background (CMB) bounds are taken into account, the amplitude of this signal
is too low to be observed at LISA (as well as all current and next generation
GW detectors). However, several well motivated inflationary models have been
proposed that can lead to an observable signal; see _e.g.,_ refs.
Bartolo:2016ami ; LISACosmologyWorkingGroup:2022jok  for an extensive list of
references. In this work we present a set of templates that can mimic a large
number of these proposed signals, motivated by theoretical considerations on
existing inflationary scenarios. For each of these templates, we build a
dedicated LISA data analysis pipeline, reconstructing the cosmological SGWB
signal as well as the expected astrophysical foregrounds and instrumental
noise. With this pipeline, we forecast the accuracy at which LISA will
reconstruct the parameters characterising each template, if the signal is
drawn from it. For some illustrative cases, we show how the reconstruction can
shed light on the fundamental parameters of such inflationary setups.

Due to the nature of the global fit, definitive conclusions on the LISA
capabilities for the detection and reconstruction of the cosmological SGWB
require data analysis pipelines for all possible LISA sources. In particular,
quantifying such capabilities requires a global-fit pipeline exploiting all
possible means to isolate the cosmological SGWB from the other sources,
_e.g.,_ tools based on its statistical properties (as for instance Gaussianity
and stationarity), anisotropic features LISACosmologyWorkingGroup:2022kbp , or
polarisation Domcke:2019zls . There is no doubt that the pipeline presented in
this paper is just one of the first steps towards the final result. It is
however sufficient for our main purposes which are:

  1. 1.

To initiate a LISA template bank of well-motivated SGWB models justified by
inflationary mechanisms.

  2. 2.

To quantify the ballpark of inflation-model parameter space that LISA can
probe with excellent scientific insight, i) if ESA delivers a very accurate
LISA noise model, ii) if the astrophysical community LISA:2022yao  manages to
precisely model the astrophysical foregrounds, iii) and if the data analysis
and waveform communities LISAConsortiumWaveformWorkingGroup:2023arg  achieve
binary waveform reconstruction with residuals that do not mimic too strongly a
SGWB signal.

  3. 3.

To estimate the accuracy at which LISA can reconstruct the SGWB templates,
provided the three “ifs” in 2) are met, as a function of their parameters.

The latter is particularly important for the theory community, as it gives
some accuracy maps that the SGWB theoretical predictions should aim at.
Indeed, theoretical uncertainties above the estimated accuracies would risk
introducing biases, and then jeopardising the reconstruction of both the
cosmological SGWB and all the other GW sources.

The outline of the paper is as follows. In section 2 we propose a list of
templates that approximate wide classes of SGWB originated by inflationary
mechanisms. For each template, we single out some specific inflationary setups
that illustrate concrete examples of the physics and their fundamental-
parameter versus template-parameter maps. Section 3 describes how we extend
the SGWBinner code Caprini:2019pxz ; Flauger:2020qyi  in order to perform
searches and parameter estimations for the SGWB templates. In section 4 we use
the extended SGWBinner to identify the parameter space that LISA can probe for
the considered templates in the presence of galactic and extragalactic
foregrounds. We also forecast the striking accuracy at which the inflationary
SGWB can be reconstructed, provided sufficient knowledge on three LISA areas
is reached (specifically, the modelling of the LISA noise, the foregrounds and
the primordial SGWB). In section 5 we map the forecast reconstructions of the
template parameters to some of the parameter space of inflationary setups. We
devote section 6 to our conclusions. In appendix A we show how changes in the
parameters of each template influence the spectral shape of the SGWB. In
appendix B we include further technical details on SGWB sourced at second
order by scalar fluctuations.

##  2 SGWB template databank for inflationary mechanisms

The inflationary mechanism can be realised in a variety of different models.
In the following subsections, we introduce a selection of templates that
describe well-motivated inflationary scenarios able to amplify the SGWB
spectrum in the LISA frequency band. The resulting list of templates does not
exhaust all possibilities, but aims at initiating the first LISA template
databank for SGWBs from inflation.

It is worth noting that the templates we present do not have a unique
parameterisation. However, for efficient numerical analysis, it is convenient
to choose parameterisations that minimize degeneracies and correlations. One
way to achieve this is by focusing on parameters that better characterise the
frequency shape and the ‘geometrical’ features (such as _e.g.,_ the high or
low frequency tilt, the position of a peak, etc) of the SGWB signal: we follow
this guiding principle in our selection of templates below. In appendix A we
show how the frequency profile of the different templates is affected by the
change in their parameters. This is aimed at offering a visual understanding
of how the template depends on its parameters.

In this paper, we will express the SGWB power spectrum in terms of
$\Omega_{\rm GW}$, defined as the GW fractional energy density per unit
frequency logarithmic interval,

$\Omega_{\rm GW}\left(f\right)\equiv\frac{1}{\rho_{0,\rm
crit}}\,\frac{d\rho_{\rm GW}}{d\ln f}\,,$ (1)

where $\rho_{0,\rm crit}\equiv 3H_{0}^{2}/8\pi G$ is the present time critical
energy density of a flat universe, and $H_{0}\equiv h\times 100\,{\rm
km}\,{\rm s}^{-1}\,{\rm Mpc}^{-1}$ is the present time value of the Hubble
parameter. As customary, we multiply $\Omega_{\rm GW}\left(f\right)$ by the
dimensionless Hubble parameter squared $h^{2}$ so as not to propagate
uncertainties related to its estimate in our analysis.

###  2.1 Power law

The first template we consider is for a power-law (PL) signal

$h^{2}\Omega^{\rm PL}_{\textrm{GW}}(f,\vec{\theta}_{\rm
cosmo})=h^{2}\Omega_{*}\left(\frac{f}{f_{*}}\right)^{n_{t}}\;,$ (2)

parameterised in terms of the parameters $\vec{\theta}_{\rm
cosmo}=\\{\Omega_{*},f_{*},n_{t}\\}$. Here, $n_{t}$ is the tilt of the power
law, $f_{*}$ is a pivot frequency, $\Omega_{*}$ parameterises the amplitude of
the spectrum at the frequency $f_{*}$. The two parameters $f_{*}$ and
$\Omega_{*}$ are degenerate so one of them should be fixed in the analysis. We
fix $f_{*}=1\,$mHz. This choice minimizes the correlation between $\Omega_{*}$
and $n_{t}$ in the analysis around $f_{*}$, _i.e.,_ close to the frequency at
which LISA reaches its best sensitivity ($\simeq 3\,$mHz).

The PL template has a simple parameter dependence and can be analyzed
relatively quickly. For this reason, we see no need for setting a stringent
prior to optimize the search. We thus take a log-uniform prior on
$\Omega_{*}\in[10^{-30},\,10^{-5}]$ and a flat prior on $n_{t}\in[-10,\,10]$,
keeping however in mind that signals with
$h^{2}\Omega_{\textrm{GW}}\gtrsim\mathcal{O}\left(10^{-6}\right)$ would be in
tension with the big-bang nucleosynthesis (BBN) bound666 In order for our
priors to be as uninformative as possible, hereafter we will always provide a
broad prior range on the amplitude of the SGWB, including values that are
excluded by the BBN+CMB bound $h^{2}\,\int_{f_{\rm min}}^{f_{\rm max}}\,{\rm
d}\ln f\,\Omega_{\rm GW}(f)\,<\,1.2\times 10^{-6}$ at 95% CL, see
Pagano:2015hma . In practice, however, we will never adopt benchmarks that are
ruled out by this bound and therefore we will not discuss this issue anymore.
LISACosmologyWorkingGroup:2022jok . The prior on the tensor spectral index is
motivated by the fact that $|n_{t}|\lesssim{\cal O}(1)$ in all models we
consider. Although, strictly speaking, for our purposes the parameterisation
in eq. 2 needs to be valid only in the LISA observational window (so that
$n_{t}$ might be negative within this window) we stress that an overall
growing trend of the spectrum with increasing frequency is needed to have an
inflationary signal visible in LISA, while being compatible with CMB bounds
LISACosmologyWorkingGroup:2022jok .

In general, we expect that any SGWB resembles a PL at frequencies far away
from those corresponding to the typical scales involved in the processes
sourcing the GWs, or those diluting their energy density throughout cosmic
history. In the following, we link explicitly the PL template and some
illustrative, well-motivated inflationary mechanisms.

Axion inflation.

In this setup, the inflaton $\phi$ is an axion coupled to a gauge field
through the interaction $\phi F{\tilde{F}}/(4f_{\phi})$, with $f_{\phi}$ being
the axion decay constant. The rolling axion strongly amplifies the gauge
field, which in turn generates a SGWB with large amplitude Barnaby:2010vf ;
Sorbo:2011rz . The production is exponentially sensitive to the parameter
$\xi\equiv{|\dot{\phi}|}/({2f_{\phi}H})$ (from now on, a dot denotes
differentiation with respect to cosmic time) whose growth during inflation
strongly tilts the GW spectrum to blue. The overall frequency shape of the
SGWB from CMB to interferometer scales is not a single power-law. Note that,
in the limit in which the two slow-roll parameters $\epsilon_{*}$ and
$\eta_{*}$ are hierarchical, eq. (6) provides a one-to-one correspondence
between $(\Omega_{*},\,n_{T})$ and the relevant input parameters:

$\displaystyle\xi_{*}\simeq-\frac{3}{2\pi}\,W\left(-\frac{0.015}{\sqrt[6]{\frac{h^{2}\Omega_{*}M_{\text{pl}}}{H_{*}}}}\right)\;,\;\;\left\\{\begin{array}[]{l}\eta_{*}\simeq-\frac{n_{T}}{2\left(2\pi\xi_{*}-3\right)}\;\;,\;\;\eta_{*}\gg\left|\epsilon_{*}\right|\\\
\epsilon_{*}\simeq\frac{n_{T}}{2\left(2\pi\xi_{*}-5\right)}\;\;,\;\;\left|\epsilon_{*}\right|\gg\eta_{*}\end{array}\right.\;,$
(5)

where $W(x)$ is the Lambert function. But the PL template is adequate for the
part of the signal falling in the LISA frequency window, with amplitude and
spectral tilt given by Barnaby:2010vf ; Sorbo:2011rz ; Bartolo:2016ami

$\displaystyle h^{2}\Omega_{*}$ $\displaystyle\simeq$ $\displaystyle 1.5\times
10^{-13}\,\frac{H_{*}^{4}}{M_{\rm Pl}^{4}}\,\frac{{\rm
e}^{4\pi\xi_{*}}}{\xi_{*}^{6}}\;,$ $\displaystyle n_{T}$ $\displaystyle\simeq$
$\displaystyle-4\epsilon_{*}+\left(4\pi\xi_{*}-6\right)\left(\epsilon_{*}-\eta_{*}\right)\;.$
(6)

Here $\epsilon\equiv-\dot{H}/H^{2}$ and $\eta\equiv-\ddot{\phi}/(H\dot{\phi})$
are the first and second slow-roll parameter respectively, and the star
indicates that they are evaluated when the mode of pivot frequency $f_{*}$
leaves the horizon during inflation. Namely, at $k_{*}\equiv{2\pi
f_{*}}=a_{*}H_{*}$, with $a_{*}$ being the scale factor of the universe at
that time. The expression for the spectral tilt in eq. (6) assumes a steady-
state evolution of the parameter $\xi$, which holds in the regime of weak
backreaction of the produced gauge fields on the inflaton dynamics, but also
at the transition between the weak and the strong backreaction regime, during
which the SGWB spectrum might rise to observable levels.777The result for
modes that leave the horizon during the strong backreaction regime is still
under investigation Garcia-Bellido:2023ser ; Caravano:2022epk ;
Figueroa:2023oxc  due to the nontrivial evolution of $\xi$ in this regime
Peloso:2022ovc . Depending on the parameters of the model, it is possible that
this uncertain regime manifests itself at high frequencies outside the LISA
window. Note that, in the limit in which the two slow-roll parameters
$\epsilon_{*}$ and $\eta_{*}$ are hierarchical, eq. (6) provides a one-to-one
correspondence between $(\Omega_{*},\,n_{T})$ and the relevant input
parameters:

$\displaystyle\xi_{*}\simeq-\frac{3}{2\pi}\,W\left(-\frac{0.015}{\sqrt[6]{\frac{h^{2}\Omega_{*}M_{\text{pl}}}{H_{*}}}}\right)\;,\;\;\left\\{\begin{array}[]{l}\eta_{*}\simeq-\frac{n_{T}}{2\left(2\pi\xi_{*}-3\right)}\;\;,\;\;\eta_{*}\gg\left|\epsilon_{*}\right|\\\
\epsilon_{*}\simeq\frac{n_{T}}{2\left(2\pi\xi_{*}-5\right)}\;\;,\;\;\left|\epsilon_{*}\right|\gg\eta_{*}\end{array}\right.\;,$
(9)

where $W(x)$ is the Lambert function.

Massive graviton during inflation.

The breaking of space diffeomorphisms can give rise to a massive graviton
during the inflationary epoch; see _e.g.,_ refs. Endlich:2012pz ;
Ricciardone:2016lym ; Bartolo:2015qvr ; Cannone:2014uqa ; Cai:2014uka ;
Lin:2015cqa ; Cannone:2015rra ; Graef:2015ova ; Domenech:2017kno ;
Ricciardone:2017kre ; Celoria:2017idi ; Dimastrogiovanni:2018uqy ;
Fujita:2018ehq . Such a graviton mass $m_{h}$ has the effect of tilting the
SGWB spectrum towards the blue, making the signal potentially detectable with
LISA, while below the current bound at CMB frequencies. Such a tilt is related
to the graviton mass by

$\displaystyle n_{T}\,=3-2\sqrt{\frac{9}{4}-\frac{m_{h}^{2}}{H^{2}}}\,\,,$
(10)

where we have focused on the vanishing $\epsilon$ limit (see _e.g.,_ ref.
Ricciardone:2016lym ). A connection between the amplitude of the SGWB and the
model parameter is, however, more complicated to draw. Indeed, while the tilt
is robustly given by eq. (10), the overall normalization of the SGWB is model
dependent and controlled by several free parameters Endlich:2012pz ;
Ricciardone:2016lym ; Bartolo:2015qvr ; Cannone:2014uqa ; Cai:2014uka ;
Lin:2015cqa ; Cannone:2015rra ; Graef:2015ova ; Domenech:2017kno ;
Ricciardone:2017kre ; Celoria:2017idi ; Dimastrogiovanni:2018uqy ;
Fujita:2018ehq .

Finally, let us also emphasize that, in order for the square root in eq. (10),
not to be imaginary the graviton mass only needs to satisfy $m_{h}/H<3/2$. We
also point out that values much smaller than $3/2$ are apparently in tension
with the Higuchi bound $m_{h}/H>\sqrt{2}$ Higuchi:1986py . Nevertheless, the
Higuchi bound relies on exact de Sitter isometries: the latter can be
spontaneously broken in the aforementioned scenarios where scalars acquire
spatial vacuum expectation values, hence allowing for the Higuchi bound to be
violated Ricciardone:2017kre .

Time-dependent tensor sound speed.

The sound speed of gravitational tensor modes (or other helicity-2 degrees of
freedom coupled to it) can exhibit a time dependence in inflationary setups
with non-minimally coupled scalars Kobayashi:2011nu ; Creminelli:2014wna  or
with extra helicity-2 component of spin-2 fields Iacconi:2020yxn ;
Dimastrogiovanni:2021mfs . These phenomena have implications for CMB
polarisation experiments Raveri:2014eea , but they can also have interesting
consequences for GW detections, leading _e.g.,_ to a growth of the GW spectrum
at small scales (see _e.g.,_ ref. Bartolo:2016ami  and references therein).

In these scenarios, it is assumed that the tensor speed becomes equal to the
speed of light at the end of inflation. For the case of single-field slow-roll
inflation plus heavy helicity-2 components, the produced SGWB in the LISA band
approximately follows eq. (2) with $\Omega_{*}$ being a non-trivial function
of the sound speed Iacconi:2020yxn ; Dimastrogiovanni:2021mfs  and $n_{t}$
reading

$n_{t}\simeq|s_{2}|(2\nu-1)\;,$ (11)

where $s_{2}=\dot{c}_{2}/(Hc_{2})$, a parameter which we assume constant
during inflation Iacconi:2020yxn , and $\nu=\sqrt{9/4-m_{h}^{2}/H^{2}}$. We
recall that here $m_{h}$ is the graviton mass, while $c_{2}$ is the time-
dependent sound speed of the additional helicity-2 field.

###  2.2 Log-normal bump

For a spectral shape with a log-normal bump (LN), we adopt the template

$h^{2}\Omega^{\rm LN}_{\textrm{GW}}(f,\vec{\theta}_{\rm
cosmo})=h^{2}\Omega_{*}\,{\rm
exp}\left[-\frac{1}{2\rho^{2}}\,\log_{10}^{2}\left(\frac{f}{f_{*}}\right)\right]\;,$
(12)

where the parameters $\vec{\theta}_{\rm cosmo}=\\{\Omega_{*},f_{*},\rho\\}$
control, respectively, the height, the position, and the width of the bump.

Due to the simplicity of the template, we assume the log-uniform priors
$h^{2}\Omega_{*}\in[10^{-30},\,10^{-5}]$, $\rho\in[0.01,\,10]$ (ranging from a
narrow to a broad bump, including values that have been considered in the
literature; see below for one specific example), and
$f_{*}\in[10^{-5},\,10^{-1}]$ Hz, so that the bump is detectable in the LISA
window.

The LN template typically covers inflationary phenomena where the GW
production is maximal at the scale $f_{*}$, which then persists for some time
during inflation. The bump is narrow (respectively, broad) for small
(respectively, large) values of $\rho$, corresponding to a shorter
(respectively, longer) stage of GW production. We now discuss a concrete model
motivating the LN template.

Axion spectator.

While the inflaton slowly rolls throughout inflation, there might also exist a
‘spectator’ axion field $\chi$, which undergoes rolling for a certain number
of e-folds $\Delta N$ of inflation Namba:2015gja . During this phase, the
axion can excite gauge fields through a $\chi F{\tilde{F}}$ interaction: in
turn, the excited gauge fields can source a sizeable SGWB. This production
mechanism is similar to the axion inflation mechanism described above, with
the distinction that now the SGWB exhibits a peak at the frequency $f_{*}$,
corresponding to the typical wavenumber scale $k_{*}$ that left the horizon
during the period of fastest rolling of the field $\chi$. The SGWB spectrum
can be approximated by the shape (12), where the height of the peak is
exponentially sensitive to $\dot{\chi}$, while the width is an increasing
function of $\Delta N$ Namba:2015gja ; Campeti:2022acx .888References
Dimastrogiovanni:2016fuu ; Thorne:2017jft ; Putti:2024uyr  studied variations
of the same class of models with a GW signal that can be well described by
this template.

In particular, based on the calculations and conventions of ref. Namba:2015gja
, if we fix the maximum parameter $\xi_{*}=5$ and the inflaton slow roll
parameter $\epsilon_{\phi}=10^{-3}$, we obtain
$\left\\{\log_{10}\Omega_{*},\,\log_{10}\rho\right\\}\simeq\left\\{-9.9,-0.13\right\\}$
for the broader bump studied in ref. Namba:2015gja , while
$\left\\{\log_{10}\Omega_{*},\,\log_{10}\rho\right\\}\simeq\left\\{-11.8,-0.37\right\\}$
for the narrower one. The moment at which the axion experiences the maximum
speed, sets $f_{*}$, and the case $f_{*}=10^{-3}$ Hz can be easily achieved.

###  2.3 Broken power law

The template

$\displaystyle h^{2}\Omega_{\rm GW}^{\rm BPL}(f,\vec{\theta}_{\rm cosmo})$
$\displaystyle=$ $\displaystyle
h^{2}\Omega_{*}\frac{\left(\frac{f}{f_{*}}\right)^{n_{{t},1}}}{\left\\{\frac{1}{2}\left[1+\left(\frac{f}{f_{*}}\right)^{1/\delta}\right]\right\\}^{(n_{t,1}-n_{t,2})\delta}}\;,$
(13)

is adequate for SGWBs with frequency shapes looking like a smooth broken power
law (BPL). Its parameters are $\vec{\theta}_{\rm
cosmo}=\\{\Omega_{*},f_{*},n_{t,1},n_{t,2},\delta\\}$, with $\Omega_{*}$
parameterising the amplitude of the peak located at frequency $f_{*}$, while
$n_{t,1}$ and $n_{t,2}$ are, respectively, the tilt of the signal tail at
$f\ll f_{*}$ and $f\gg f_{*}$, with the final parameter $\delta$ controlling
the sharpness of the transition between the two regimes.

For sufficiently large and positive $n_{t,1}$, and a negative and sufficiently
large (in absolute value) $n_{t,2}$, the amplitude of the spectrum can be
localized within the LISA frequency band, and be negligible at the frequencies
probed by other GW detectors. We take the log-uniform prior
$h^{2}\Omega_{*}\in[10^{-30},\,10^{-5}]$. For the other parameters, we take
the flat priors $n_{t,1}\in[-10,\,10]$ and $n_{t,2}\in[-10,\,10]$, and the
log-uniform prior $f_{*}\in[10^{-5},\,10^{-1}]\,$Hz. We finally take flat
priors $\delta\in[0.1,10]$ to describe both the case in which the slope
changes softly, or abruptly.

The template is motivated by, _e.g.,_ the concrete examples we discuss next.

Second slow-roll phase.

When GWs are induced at second-order by curvature perturbations sourced by a
broad and flat primordial curvature power spectrum, they lead to an SGWB
featuring a BPL frequency shape with $n_{t,2}\simeq 0$ (see appendix B for
details). Such a scenario can be achieved, for example, through an ultra-slow-
roll (USR) phase followed by a second slow-roll regime generating the required
plateau Wands:1998yp ; Leach:2000yw ; Leach:2001zf ; Biagetti:2018pjj ;
Franciolini:2022pav . This model, originally suggested in ref. DeLuca:2020agl
, could potentially connect the pulsar timing array frequency range (where
recent data provide evidence for the existence of a SGWB NANOGrav:2023gor ;
EPTA:2023fyk ; Reardon:2023gzh ; Xu:2023wog  of yet-to-be determined nature)
to signatures in the LISA frequencies associated with the formation of
primordial black holes Sasaki:2018dmp  with masses in the range
$(10^{-15}\,-\,10^{-11})M_{\odot}$, where they may comprise the totality of
the dark matter in the universe Vaskonen:2020lbd ; Franciolini:2022pav .

Hybrid inflation with mild waterfall stage.

The template above can also describe a broad bump, which is asymmetric around
$f_{*}$ if the primordial curvature power spectrum takes the shape of a broad
bump. Popular models within this category are multifield setups such as hybrid
inflation Linde:1993cn ; Garcia-Bellido:1996mdl ; Clesse:2015wea  or other
types of models featuring a slow-turn in the field space together with a
tachyonic instability of isocurvature perturbations Fumagalli:2020adf ;
Braglia:2020eai ; Braglia:2020taf . We focus on hybrid inflation for
illustration. In this setup, there are two dynamical scalar fields, the
inflaton and a so-called hybrid field. The potential of these fields is
designed to induce a second stage of inflation during which isocurvature
perturbations grow tachyonically on superhorizon scales, and are converted
into adiabatic perturbations, leading to a peak in the primordial curvature
power spectrum Garcia-Bellido:1996mdl ; Clesse:2015wea . The induced SGWB
typically takes the form of a BPL with $n_{t,1}>0$ and $n_{t,2}<0$
Clesse:2018ogk ; Spanos:2021hpk ; Braglia:2022phb . The relationship between
the BPL parameters and the parameters of the two-fields potential cannot
generally be expressed in a closed form, but some features can be generically
related to the underlying dynamics during inflation. For example, let us
consider the potential of the $\alpha$-attractor case Braglia:2022phb . We
have the following features: the linear term of the hybrid field controls
$\Omega_{*}$; the value of the hybrid field at the minimum of the potential
sets $f_{*}$; the first derivative along the hybrid direction modulates
$n_{t,2}$; and the time profile of the squared mass of the hybrid field
controls how much isocurvature modes with $f<f_{\rm peak}$ are amplified, and
hence the size of $n_{t,1}$. For instance, for the parameter values in Eq.
(5.3) of ref. Braglia:2022phb , one numerically obtains a SGWB resembling the
BPL with $\log_{10}(h^{2}\Omega_{*})=-9.3$, $f_{*}=1\,{\rm mHz}$,
$n_{t,1}=2.65$, $n_{t,2}=-2.1,\,\delta=5.3$.

###  2.4 Double peak

The template

$\begin{aligned} h^{2}\Omega^{\rm DP}_{\rm GW}(f,\vec{\theta}_{\rm
cosmo})=&\,\,h^{2}\Omega_{*}\Bigg{[}\beta\,\left(\frac{f}{\kappa_{1}f_{*}}\right)^{\\!n_{p}}\left[\frac{c_{1}-f/f_{*}}{c_{1}-\kappa_{1}}\right]^{\\!\\!\frac{n_{p}}{\kappa_{1}}\left(c_{1}-\kappa_{1}\right)}\Theta\\!\left(c_{1}-\frac{f}{f_{*}}\right)\\\
&+\exp\\!\left[-\frac{1}{2\rho^{2}}\log_{10}^{2}\\!\left(\frac{f}{\kappa_{2}f_{*}}\right)\right]\left\\{1+{\rm
erf}\left[-\gamma\log_{10}\\!\left(\frac{f}{\kappa_{2}f_{*}}\right)\right]\right\\}\Bigg{]}\end{aligned}\;,$
(14)

is aimed at reconstructing a double-peak (DP) SGWB signal. The free parameters
controlling the shape are $\vec{\theta}_{\rm
cosmo}=\\{h^{2}\Omega_{*},f_{*},\beta,\kappa_{1},\kappa_{2},\rho,\gamma\\}$,
whereas $c_{1}=\sqrt{2/3}$ and $n_{p}\simeq 2.5$ are kept constant, thanks to
universal infrared properties of the GW background, as discussed in section
B.2. The parameter $f_{*}$ sets the frequency of the transition from the low-
frequency to the high-frequency peak via the Heaviside step function $\Theta$.
The high-frequency peak is a skewed log-normal with its maximal amplitude
$h^{2}\Omega_{*}$ occurring at $f\simeq\kappa_{2}f_{*}$, with $\kappa_{2}>1$.
Its skewness and width are set by $\gamma$ and $\rho$, respectively. In
contrast, the low-frequency peak reaches its maximum amplitude $\beta
h^{2}\Omega_{*}$, with $0<\beta<1$, at frequency $f\simeq\kappa_{1}f_{*}$ with
$\kappa_{1}<1$. At $f\ll\kappa_{1}f_{*}$, the template is a power law $\propto
f^{n_{p}}$. A minimum between the two peaks arises at $f\simeq c_{1}f_{*}$
when $\rho$ is sufficiently small.999In principle, a DP frequency shape can be
described by functional forms different from the one in eq. (14). As
alternatives, we considered a ratio of polynomials and a sum of power laws
times exponentials of polynomials. We discard these alternatives in favor of
eq. 14. By performing a Bayesian analysis, both with a uniform search (equally
weighted simulated data points) and a weighted search (data weighted with the
LISA sensitivity curve), we indeed found that our choice (14) better describes
the DP SGWB signal predicted in the inflationary models that motivate the
template.

The DP SGWB profile is expected from second-order emission triggered by
enhanced scalar perturbations (see appendix B for details). A signal described
by the DP template represents the characteristic SGWB profile induced by large
(compared to CMB scales) inflationary scalar fluctuations re-entering the
horizon during a radiation-dominated phase Tomita:1975kj ; Matarrese:1993zf ;
Matarrese:1997ay ; Acquaviva:2002ud ; Mollerach:2003nq ; Nakamura:2004rm ;
Ananda:2006af ; Baumann:2007zm ; Espinosa:2018eve ; Kohri:2018awv ;
Inomata:2019yww ; Domenech:2019quo . This enhancement of scalar fluctuations
can be generated in single- and multi-field inflationary models featuring a
short period of non-attractor evolution, which, in turn, amplifies the SGWB in
the corresponding frequency range.

In addition, if the support of the scalar spectrum is narrow in momentum
space, a resonance between the transfer function of the scalar modes and the
Green function in the tensor evolution, occurring when the frequency of the
scalar source coincides with the frequency of the GWs, introduces a second
pronounced peak on the SGWB profile Ananda:2006af ; Saito:2008jc .

This broad class of scenarios can be generated by one of the following
approximated primordial curvature power spectra:

$\displaystyle\mathcal{P}_{\zeta}^{\rm ln}(k)$ $\displaystyle=$
$\displaystyle{\cal
A}_{s}\exp\\!\left[-\frac{1}{2\Delta^{2}}\ln^{2}\\!\left(\frac{k}{k_{*}}\right)\right]\qquad\textrm{(log-
normal)}\,,$ (15) $\displaystyle\mathcal{P}_{\zeta}^{\rm bpl}(k)$
$\displaystyle=$ $\displaystyle\frac{{\cal
A}_{s}(p_{1}+p_{2})}{\left[p_{2}\left(\frac{k}{k_{*}}\right)^{-p_{1}}+p_{1}\left(\frac{k}{k_{*}}\right)^{p_{2}}\right]}\qquad\textrm{(broken
power-law)}\;,$ (16)

where $\Delta,\,p_{1},\,p_{2}>0$. We prove in section B.2 and section B.3 that
both spectra lead to a SGWB frequency shape that is well fitted by the DP
template, and we use them to set the priors on the DP template. Specifically,
we determine the overall ranges of DP template parameters that allow us to fit
the SGWBs derived from $\mathcal{P}_{\zeta}^{\rm ln}(k)$ when
$\Delta\in[0.1,1]$ and $\mathcal{P}_{\zeta}^{\rm bpl}(k)$ when
$p_{2}\in[0.5,4]$ and $p_{1}=4$ (see motivations below) within the ranges
$h^{2}\Omega_{*}\in[10^{-30},10^{-5}]$,
$f_{*}\in[10^{-5},10^{-1}]\,\textrm{Hz}$, $\beta\in[0.05,0.6]$,
$\kappa_{1}\in[0.05,0.9]$, $\kappa_{2}\in[1.1,3]$, $\rho\in[0.04,0.5]$, and
$\gamma\in[0.02,20]$. We then assume log-uniform priors on $h^{2}\Omega_{*}$
and $f_{*}$, and flat priors on the other DP parameters.

Let us now motivate the two spectral shapes in eqs. (15) and (16) in terms of
concrete scenarios.

Log-normal $\mathcal{P}_{\zeta}^{\rm ln}$.

One example of peaks in the primordial curvature power spectrum that can be
parameterised by a narrow log-normal shape is provided by single-field models
where the enhancement of scalar fluctuations is given by a resonant mechanism.
This is for instance realised in scenarios with an oscillating speed of sound,
with amplitude and frequency parameterised by $(\xi,k_{*})$, that modifies the
standard evolution equation for the perturbations Cai:2018tuh ; Cai:2019jah ;
Chen:2019zza ; Chen:2020uhe ; Addazi:2022ukh . The Mukhanov–Sasaki equation
can thus be approximated as a Mathieu equation which, in turn, presents a
parametric instability for certain ranges of modes located around the
oscillatory frequency $k_{*}$ with width $\sim\xi k_{*}$. In principle, this
feature provides a direct link between the position and width of the peak in
the primordial power spectrum, and the frequency and amplitude of the speed-
of-sound oscillations.

Broken power-law $\mathcal{P}_{\zeta}^{\rm bpl}$.

A period of USR in single-field models of inflation produces a primordial
curvature power spectrum with a shape that can be approximated as in eq. (16).
The powers are determined by the second slow-roll parameter
$\eta\equiv-\ddot{H}/(2H\dot{H})$, where $H$ is the Hubble parameters and dot
indicates derivatives with respect to cosmic time. Denoting it $\eta_{1}$
during the slow-roll period preceding USR, and $\eta_{3}$ during the constant-
roll period typically following the USR, one finds Karam:2022nym :
$p_{1}=5-|1-2\eta_{1}|$ and $p_{2}=2\eta_{3}$. Typically, the conditions
$|\eta_{1}|\ll 1$ and $\eta_{1}<0$ arise, so that the growth of the spectrum
is $\mathcal{P}_{\zeta}\propto k^{4}$ Byrnes:2018txb , though a steeper growth
can also be realised in some scenarios Carrilho:2019oqg ; Fumagalli:2020adf ;
Tasinato:2020vdk ; Braglia:2020taf ; Davies:2021loj . The peak amplitude of
the curvature spectrum is determined by the second slow-roll parameter
$\eta_{2}$ during USR and the duration $\Delta N_{\rm USR}$ of the USR period
as $\mathcal{A}_{s}\sim e^{2\eta_{2}\Delta N_{\rm USR}}\mathcal{P}_{\rm CMB}$,
where $\mathcal{P}_{\rm CMB}$ denotes the amplitude at the CMB scales.
Assuming Wands duality Wands:1998yp  between the USR and the final constant-
roll, which commonly holds for smooth potentials, the second slow-roll
parameter during the USR period is $\eta_{2}=3-\eta_{3}$. Assuming the
enhanced spectrum is produced by an USR phase of inflation, one can perform a
reverse engineering procedure to connect the spectral shape to the
inflationary dynamics and the inflaton potential Franciolini:2022pav . The
location of the peak $k_{*}$ is determined by the number of e-folds between
the epoch of Hubble crossing of CMB modes and the onset of USR dynamics.

The broken power-law spectrum in eq. (16) is also realised in thermal
inflation models where the fluctuations become large around the time when the
inflaton potential turns from convex to concave. In this case, the growth is
$\mathcal{P}_{\zeta}\propto k^{3}$ and the decay and the peak amplitude are
determined by the tachyonic mass of the thermal inflaton in the false vacuum
and by the potential energy of the false vacuum Dimopoulos:2019wew .

###  2.5 Excited states

The template

$h^{2}\Omega_{\mathrm{GW}}^{\rm ES}(f,\vec{\theta}_{\rm
cosmo})=\frac{h^{2}\Omega_{*}}{0.052}\frac{1}{x^{3}}\left[1-\frac{x^{2}}{4\gamma_{\rm
ES}^{2}}\right]^{2}\left[\sin(x)-4\frac{\sin^{2}(x/2)}{x}\right]^{2}\Theta\\!\left(x_{\rm
cut}-x\right),$ (17)

where $x\equiv(f\,\omega_{\rm ES})/2$, parameterises a SGWB generated in the
presence of a scalar excited state (ES) during inflation. Its parameters are
$\vec{\theta}_{\rm cosmo}=\\{\Omega_{*},\gamma_{\rm ES},\omega_{\rm ES}\\}$,
where $\Omega_{*}$ sets the amplitude of the primary peak at frequency $f_{\rm
max}\simeq 6/\omega_{\rm ES}$, $\omega_{\rm ES}$ sets the periodicity of the
subsequent peaks (at $f>f_{\rm max}$), and $\gamma_{\rm ES}$ determines the
largest frequency $f_{\rm cut}=2/3\gamma_{\rm ES}f_{\rm max}$, corresponding
to $x_{\rm cut}=2\gamma_{\rm ES}$, where the template is applicable. As
explained below, theoretical constraints leads us to expect $\gamma_{\rm
ES}\gg 1$, but not arbitrarily large. We thus set a flat prior
$\log_{10}(\gamma_{\rm ES})\in[0.699,2]$. We take the flat prior
$\log_{10}(h^{2}\Omega_{*})\in[-30,-5]$, while we assume
$\log_{10}(\omega_{\rm ES}\,{\rm Hz})\in[2,5]$ because this range of values
provides a signal peaked in the LISA frequency band (for signals with
$\omega_{\rm ES}$ outside this range, simpler templates than the ES would be
more suitable).

A SGWB following the ES template is expected in inflationary setups where a
large number of scalar particles are produced dynamically with momenta peaked
around a narrow range of scales Fumagalli:2021mpc . Such an excited state
dynamically arises in the presence of a transient non-adiabatic evolution
during inflation. The energy-momentum tensor of the produced particles sources
the tensor modes, giving rise to an ES SGWB. The typical time of particle
production, or, equivalently, the associated frequency $f_{\rm out}$, sets the
frequency of the periodicity, with $\omega_{\rm ES}\simeq 2/f_{\rm out}$,
while the deeper inside the horizon particle production took place, the larger
$\gamma_{\rm ES}$, _i.e.,_ the larger frequencies the signal extends to.
Moreover, the higher the number of scalars in an excited state, the larger the
amplitude of the signal.

A dynamical mechanism to generate an excited state can be found, for instance,
in models of multifield inflation in which the background trajectory deviates
strongly and for a brief period from a geodesic of the field-space manifold
Palma:2020ejf ; Fumagalli:2020adf . Other realisations can be found in single-
field models with a feature in the potential Inomata:2021uqj ; Inomata:2021tpx
, parametric resonances due to a periodic structure in the potential
Peng:2021zon ; Inomata:2022yte , or with a spectator field coupled to the
inflaton via a periodic function Cai:2021wzd . Quite generically, tensor modes
with spectral shape as in eq. 17 can be generated whenever scalar fluctuations
are amplified on sub-horizon scales by some mechanism during inflation
Fumagalli:2021mpc  (see ref. Inomata:2021zel  for upper bounds). Let us delve
into this topic further by exploring a specific scenario.

Strong turn in two-field inflation.

A transient non-adiabatic evolution is provided by a brief period during which
the turn of the trajectory in field space is large. This leads to the
dynamical appearance of an excited state, and consequently to a signal
following the ES template101010This setup also predicts the formation of
primordial black holes Palma:2020ejf ; Fumagalli:2020adf .. Starting from the
generic multifield non-linear sigma models Lagrangian, one can derive the
effective action for the fluctuations at second order around a given
homogeneous background. In the latter, adiabatic and entropic degrees of
freedom are coupled by a term proportional to $\eta_{\perp}\equiv-
V_{N}/(M_{\mathrm{Pl}}H^{2}\sqrt{2\epsilon})$, which is the dimensionless
parameter measuring the turning rate of the trajectory, while $V_{N}$ is the
projection of the first derivative of the potential over the entropic
direction. A top-hat profile for the time dependence of $\eta_{\perp}$, with
strength $\eta_{\perp}\gg 1$ and duration less than one e-folds, provides a
useful idealization of a strong and sharp turn. Note that, in order to satisfy
backreaction and perturbativity bounds, $\eta_{\perp}$ cannot be arbitrarily
large Fumagalli:2020nvq , motivating our previous priors choices. In this
example, the ES template parameters are related to the fundamental ones as

$\gamma_{\rm ES}\simeq\eta_{\perp},\quad\omega_{\rm ES}\simeq\frac{2}{f_{\rm
out}}\,,$ (18)

where $f_{\rm{out}}$ is the frequency corresponding to the mode leaving the
horizon at the end of a turn of strength $\eta_{\perp}$. Furthermore, the
overall amplitude can be expressed as

$h^{2}\Omega_{*}=0.03\,r_{i}\,\mathcal{N}^{4}\,\gamma_{\rm
ES}^{5}\,|\beta|^{4}\,H_{*}^{4}/(\pi M_{\rm Pl})^{4},$ (19)

with $r_{i}$ a redshift factor, $H_{*}$ the value of the Hubble parameter at
the turn, $\mathcal{N}$ the number of scalars affected by the excited state,
and $|\beta|^{2}\gg 1$ the large occupation number of particles
Fumagalli:2021mpc .

###  2.6 Linear oscillations

The linear oscillations (LO) template

$h^{2}\Omega_{\textrm{GW}}^{\textrm{LO}}(f,\vec{\theta}_{\rm
cosmo})=\Big{[}1+\mathcal{A}_{\textrm{lin}}\cos\big{(}\omega_{\textrm{lin}}f+\theta_{\textrm{lin}}\big{)}\Big{]}h^{2}\Omega_{\textrm{GW}}^{\,\rm
env}(f,\vec{\theta}_{\rm env})\;,\\\ $ (20)

describes oscillations periodic in $f$ that modulate a smooth envelope. Its
parameters are $\vec{\theta}_{\rm cosmo}=\\{\vec{\theta}_{\rm
env},\mathcal{A}_{\mathrm{lin}},\omega_{\mathrm{lin}},\theta_{\mathrm{lin}}\\}$,
with $\vec{\theta}_{\rm env}$ describing the overall envelope
$\Omega_{\textrm{GW}}^{\,\rm env}(f)$, and the other three controlling the
relative amplitude, the period and the phase of the oscillations. Since the
oscillations and the envelope signal are theoretically associated with a
common sourcing mechanism (see below), we consider a profile for
$\Omega_{\textrm{GW}}^{\,\rm env}(f)$ exhibiting a dominant, narrow peak. In
particular, linear oscillations arise in the context of GWs induced at second
order by scalar fluctuations, when the enhancement of the latter is narrowly
peaked around a given scale. Thus, the natural envelope to use would be the
double peak template of section 2.4 augmented with the fact that linear
oscillations would appear only on the highest peak Fumagalli:2020nvq . The
latter being given by a log-normal bump (modulo the skewness parameter
$\kappa_{2}$ that we neglect here for simplicity) we restrict ourselves, in
the analysis below, to a log-normal envelope, _i.e.,_
$\Omega_{\textrm{GW}}^{\,\rm env}(f)=\Omega_{\textrm{GW}}^{\rm LN}(f)$, with
$\rho\equiv 10^{\Delta}$, although other envelopes can be considered. As
discussed in section 2.4, $\rho$ has values within an order of magnitude, so
there is no need for a log prior for the width of the envelope.

For the envelope, we make use of the log-normal bump template of section 2.2
and we choose the flat priors $\Omega_{*}\in[-30,-5]$, $\log_{10}(f_{*}/{\rm
Hz})\in[-5,-1]$, and $\rho\in[0.05,1]$. For the oscillatory part we set the
flat priors $\mathcal{A}_{\rm lin}\in[0.05,1]$ and $\omega_{\rm
lin}\,\in[0.1,10^{3}]\,\mathrm{mHz}^{-1}$.

The template (20) well describes the SGWB generated at horizon re-entry of
primordial density fluctuations with power spectrum Fumagalli:2020nvq

$\displaystyle\mathcal{P}_{\zeta}(f)=\overline{\mathcal{P}}_{\zeta}(f)\bigg{[}1+A_{\textrm{lin}}\cos\Big{(}\omega
f+\phi_{\textrm{lin}}\Big{)}\bigg{]}\,,$ (21)

where we traded wavenumbers $k$ for frequencies $f$. Despite this phenomenon
being inherently nonlinear, the sinusoidal modulations in the power spectrum
get processed into corresponding sinusoidal oscillations in
$\Omega_{\textrm{GW}}$ with $\omega_{\textrm{lin}}=c_{s}^{-1}\omega$
Fumagalli:2020nvq ; Witkowski:2021raz , where $c_{s}$ is the propagation speed
of density fluctuations, equal to $1/\sqrt{3}$ in the conventional scenario of
radiation domination at horizon reentry. The oscillations do get averaged out,
though, so that the relative amplitude $\mathcal{A}_{\textrm{lin}}$ in eq.
(20) is typically $\mathcal{O}(20\%)$, even for $A_{\textrm{lin}}\simeq 1$
(with a nonstandard expansion history, the amplitude can be larger up to
$\mathcal{O}(40\%)$ Witkowski:2021raz ). Finally, the envelope
$\overline{\Omega}_{\textrm{GW}}(f)$, is determined by the envelope
$\overline{\mathcal{P}}_{\zeta}(f)$, whose details differ between different
inflationary models.

Sharp features.

The power spectrum (21) is characteristic of inflationary models with a sharp
feature, _i.e.,_ a transition during inflation that occurs on a timescale
approximately smaller than one e-fold Slosar:2019gvt . We can unequivocally
assign a scale $k_{\textrm{f}}$ to these sharp events, the momentum of modes
crossing the Hubble radius at that time, or equivalently the corresponding
frequency $f_{\textrm{f}}$. Generically, sharp features generate particle
production. When the latter is significant, it leads to an enhancement of the
primordial power spectrum at the corresponding scales, modulated by
$\mathcal{O}(1)$ oscillations ($A_{\textrm{lin}}\simeq 1$) with frequency
$\omega=2/f_{\textrm{f}}$ Fumagalli:2020nvq . As the particle production is
effective for scales inside the Hubble radius at the time of the features
($f\gg f_{\textrm{f}}$), the range of scales for which the envelope
$\overline{\mathcal{P}}_{\zeta}(f)$ is significant spans several periods of
the oscillations.

Sharp features are not associated with a single model of inflation, but can
rather be realised through a variety of mechanisms Slosar:2019gvt . Large
features leading to significant particle production and to the power spectrum
(21), and hence the template (20), can occur in single-field inflation,
_e.g.,_ caused by a step in the inflaton potential Inomata:2021uqj ;
Dalianis:2021iig , or in multifield settings, _e.g.,_ when the inflationary
trajectory exhibits a sharp turn Palma:2020ejf ; Fumagalli:2020adf ;
Aragam:2023adu . It is difficult, in general, to give analytic expressions of
the power spectrum in terms of model parameters and one has to resort to
numerical computations. Hence, the parameters appearing in eq. (21) can be
thought as the ‘fundamental’ ones that one may be interested in
reconstructing, and chief amongst them is the frequency $\omega$, indicative
of the time at which the feature arises.

The example in section 2.5 of a strong sharp turn in two-field inflation
provides a useful illustration. With the same notations, one finds
Palma:2020ejf ; Fumagalli:2020nvq

$\displaystyle\overline{\mathcal{P}}_{\zeta}(f)$ $\displaystyle={\cal
P}_{0}\frac{e^{2\sqrt{(2-\kappa)\kappa}\,\eta_{\perp}\delta}}{4(2-\kappa)\kappa},\qquad
A_{\textrm{lin}}=1,\quad$ (22) $\displaystyle\omega$
$\displaystyle=\frac{2e^{-\delta/2}}{f_{\textrm{f}}},\qquad\phi_{\textrm{lin}}=2\arctan\left[\frac{\kappa}{\sqrt{(2-\kappa)\kappa}}\right]+\pi\;,$

which is valid for $\kappa\equiv f/(f_{\textrm{f}}\,\eta_{\perp})<2$, and
where ${\cal P}_{0}$ denotes the amplitude of the power spectrum without turn
(a large enhancement limit has already been taken in these formulæ). In this
example, the peak of the $\Omega_{\rm GW}$ envelope is given by
$f_{*}=2c_{s}\eta_{\perp}f_{\textrm{f}}=4c_{s}\eta_{\perp}e^{-\delta/2}/\omega$.
Let us recall that the frequency $\omega_{\rm lin}$ of the oscillations in
$\Omega_{\rm GW}$ depends generically on the one in the primordial power
spectrum $\omega$ through $\omega_{\rm lin}=c_{s}^{-1}\omega$. We thus have
the following relation between the peak frequency of the envelope and
$\omega_{\rm lin}$:

$\omega_{\rm lin}=4\eta_{\perp}e^{-\delta/2}/f_{*}.$ (23)

In general, the phenomenon of sharp feature leads to $\omega_{\rm
lin}=\mathcal{O}(10-100)/f_{*}$.

###  2.7 Logarithmic resonant oscillations

The logarithmic resonant oscillations (RO) template is given by

$\displaystyle h^{2}\Omega_{\textrm{GW}}^{\textrm{RO}}(f,\vec{\theta}_{\rm
cosmo})$ $\displaystyle=$
$\displaystyle\Big{\\{}1+\mathcal{A}_{1}(A_{\log},\omega_{\rm
log})\cos\big{[}\omega_{\textrm{log}}\ln(f/\textrm{Hz})+\theta_{\textrm{log},1}\big{]}$
$\displaystyle\hskip 15.07993pt+\mathcal{A}_{2}(A_{\log},\omega_{\rm
log})\cos\big{[}2\omega_{\textrm{log}}\ln(f/\textrm{Hz})+\theta_{\textrm{log},2}\big{]}\Big{\\}}h^{2}\Omega^{\rm
env}_{\textrm{GW}}(f,\vec{\theta}_{\rm env})\,,$

with Fumagalli:2021cel

$\displaystyle\mathcal{A}_{1}$
$\displaystyle=\frac{A_{\textrm{log}}\mathcal{C}_{1}(\omega_{\textrm{log}})}{1+A_{\textrm{log}}^{2}\mathcal{C}_{0}(\omega_{\textrm{log}})},\qquad\theta_{\textrm{log},1}=\phi_{\textrm{log}}+\theta_{\textrm{log},1}(\omega_{\textrm{log}}),$
(25) $\displaystyle\mathcal{A}_{2}$
$\displaystyle=\frac{A_{\textrm{log}}^{2}\mathcal{C}_{2}(\omega_{\textrm{log}})}{1+A_{\textrm{log}}^{2}\mathcal{C}_{0}(\omega_{\textrm{log}})},\qquad\theta_{\textrm{log},2}=2\phi_{\textrm{log}}+\theta_{\textrm{log},2}(\omega_{\textrm{log}})\,,$

where $\mathcal{C}_{0,1,2}(\omega_{\textrm{log}})$ and
$\theta_{\textrm{log},1,2}(\omega_{\textrm{log}})$ are numerical functions
that depend on the cosmic expansion at the time the SGWB was produced
Witkowski:2021raz ; these functions are shown in fig. 1 under the assumption
of a conventional radiation-domination era. For reasons explained below, we
choose a flat envelope for $h^{2}\Omega^{\rm env}_{\textrm{GW}}$. The model
parameters of the RO template are thus
$\vec{p}_{\textrm{RO}}=\\{\Omega_{*},A_{\textrm{log}},\omega_{\textrm{log}},\phi_{\textrm{log}}\\}$,
with the latter three parameters controlling the relative amplitude, the
period, and the phase of the oscillations.

The RO template describes log-periodic oscillations, of frequencies
$\omega_{\textrm{log}}$ and $2\omega_{\textrm{log}}$, that modulate a smooth
envelope (scenarios with more general logarithmic oscillations exist
Calcagni:2023vxg ). A noteworthy aspect of the template is the qualitative
change in the oscillatory structure of section 2.7 as a function of the
frequency $\omega_{\textrm{log}}$. For values smaller than the critical
frequency $\omega_{\textrm{log,c}}\simeq 4.77$, the oscillation with frequency
$\omega_{\textrm{log}}$ dominates over the one with double frequency, while
the situation is reversed for large frequency values, with the precise cross-
over value ${\cal O}(1)\omega_{\textrm{log,c}}$ depending on
$A_{\textrm{log}}$ Fumagalli:2021cel . We set flat priors as follows:
$\log_{10}(h^{2}\Omega_{*})\in[-30,-5]$, $\log_{10}(A_{\rm log})\in[-3,0]$,
$\log_{10}(\omega_{\rm log})\in[0,2]$ and $\phi_{\rm
log}\in[-\pi,\pi]$.111111Note that $\omega_{\textrm{log}}$ is well below the
upper bound $\omega_{\textrm{log}}^{\rm max}=1/\Delta\ln f\simeq f_{*}/\Delta
f\approx 10^{4}$ due to the frequency resolution $\Delta f$ of LISA around
$f_{*}\approx 2\times 10^{-3}\,{\rm Hz}$ Calcagni:2023vxg . For the linear
oscillations of section 2.6, one has $\omega_{\textrm{lin}}^{\rm max}=1/\Delta
f\approx 5\times 10^{6}\,{\rm Hz}^{-1}$, again much larger than the
theoretical prior reported above eq. 21.

Figure 1: Functions $\mathcal{C}_{0,1,2}(\omega_{\textrm{log}})$ and
$\theta_{\textrm{log},1,2}(\omega_{\textrm{log}})$ appearing in the template
(25) and as computed in ref. Fumagalli:2021cel  for GWs induced during a
period of radiation domination.

The RO template captures well the properties of the SGWB induced by primordial
density fluctuations with power spectrum

$\displaystyle\mathcal{P}_{\zeta}(f)=\overline{\mathcal{P}}_{\zeta}(f)\Big{[}1+A_{\textrm{log}}\cos\big{(}\omega_{\textrm{log}}\ln(f/f_{*})+\phi_{\textrm{log}}\big{)}\Big{]}\,,$
(26)

especially for the scales where $\Omega_{\textrm{GW}}$ is maximal
Fumagalli:2020nvq ; Fumagalli:2021cel . For simplicity, here we concentrate on
a flat envelope $\overline{\mathcal{P}}_{\zeta}(f)$ (see ref.
Fumagalli:2021dtd  for motivations), in which case the mapping from eq. (26)
to eq. (2.7) is exact.

Power spectra of the form (26) are characteristic of inflationary models in
which some components of the background oscillate with a frequency
$\tilde{\omega}$ larger than the Hubble scale, inducing a resonance with the
oscillations of the quantum modes of the density perturbations, and resulting
in the spectrum (26) with $\omega_{\textrm{log}}=\tilde{\omega}/H>1$
Chen:2008wn . These so-called resonant features have been extensively studied
on cosmological scales Slosar:2019gvt , where observations mandate a small
amplitude of oscillations $A_{\textrm{log}}$, a typical example being axion
monodromy inflation Silverstein:2008sg ; Flauger:2009ab ; Parameswaran:2016qqq
; Ozsoy:2018flq . Scenarios with larger amplitude of the oscillations, and of
the power spectrum, have started to be studied only recently, see _e.g.,_
refs. Fumagalli:2021cel ; Bhattacharya:2022fze ; Mavromatos:2022yql .

Axion-type inflation.

Explicit models with large resonant features have been studied in refs.
Bhattacharya:2022fze ; Mavromatos:2022yql . Both are multifield axion models
that share the same qualitative features, namely, periodic modulations of the
inflationary potential due to subleading nonperturbative corrections, which
lead to periodic violations of slow-roll conditions and strong sharp turns of
the inflationary trajectory, resulting in the resonant amplification of
fluctuations (see also ref. Boutivas:2022qtl ). In both cases, an analytical
understanding of the link between the microphysical parameters of the
inflationary model and the power spectrum is lacking. However, just like the
linear-oscillations template, the parameters appearing in eq. (26) can be
thought of as the main ones to be reconstructed, since this expression
encompasses several realisations of the resonance mechanism. In particular,
the frequency $\omega_{\textrm{log}}$ is of prime interest, as it indicates in
a model-independent manner the existence of periodic modulations of the
inflationary Lagrangian with frequency $\omega_{\textrm{log}}$ in e-fold.

###  2.8 Deformations of the above templates, due to additional physics

The embedding of a given inflationary setup may include independent mechanisms
that act during or after inflation and deform the SGWB frequency shape that is
predicted in their absence. By testing such deformations, LISA can probe
deviations from the standard model of cosmology and its particle content.
Including such deformations in the LISA SGWB template bank would make it much
larger, with benefits and costs that should be carefully assessed. For this
reason, we here list some plausible deformations to remind us of their
existence and implications for the template bank, but we leave their detailed
analysis to the future.

####  2.8.1 Pre-BBN cosmology with matter-domination or kination eras

The PL template (2) and the BPL template (13) can be used to describe the
effect of the Hubble expansion rate after inflation. When the expansion
evolution of the universe differs from the radiation-dominated phase, the
primordial tilt $n_{t,\rm{prim}}$ is modified by the equation of state of the
universe $w$ as Seto:2003kc ; Boyle:2005se ; Boyle:2007zx ; Nakayama:2008ip ;
Kuroyanagi:2011fy ; Figueroa:2019paj ; ValbusaDallArmi:2020ifo

$n_{t}=n_{t,\rm{prim}}+\frac{2\left(3w-1\right)}{3w+1}\;.$ (27)

We consider the equation of state $w$ as an additional parameter, within a
reasonable range $w\in[-1/3,\,1]$. In fact $w<-1/3$ would correspond to an
accelerated expansion during which the behavior of GWs is different from the
one we consider, while $w>1$ is physically not motivated as it implies a sound
speed faster than the speed of light. An interesting case is when the
transition scale between two different $w$’s falls within the frequency band
in which case LISA may be able to observe the bend of the spectrum. In this
case, the spectral tilts $n_{t,1}$ and $n_{t,2}$ are described by the equation
of state after and before the transition $w_{1}$ and $w_{2}$ as in eq. 27. The
transition frequency $f_{*}$ corresponds to the temperature of the universe
$T_{*}$ when the transition takes place,

$f_{*}=2.6\times 10^{-7}\,\frac{T_{*}}{\rm
GeV}\,\left(\frac{g_{*}}{106.75}\right)^{1/2}\,\left(\frac{g_{*s}}{106.75}\right)^{-1/3}\,{\rm
Hz}\;,$ (28)

where $g_{*}$ and $g_{*s}$ are the effective number of relativistic degrees of
freedom contributing to the radiation density and to the entropy density,
respectively. We have normalized these quantities to the value they assume
when all the Standard Model degrees of freedom are relativistic. The relation
(28) indicates that LISA can probe physics near the $100$ TeV scale.

Any change in the Hubble expansion rate generically affects GWs generated
outside the horizon during inflation. Thus, in general, the SGWB can be a
probe of a non-standard Hubble expansion history after inflation. We note that
GWs generated at the horizon entry, _e.g.,_ induced GWs, show different
features Domenech:2020kqm  because the source term is also affected by a
change in the Hubble expansion rate, as mentioned above for the LO and RO
templates Fumagalli:2021dtd . The following are two examples of this
mechanism; in both cases, we take $n_{t,1}\neq n_{t,2}$, and both positive.

Early matter phase due to inflaton oscillation.

An early matter-dominated phase can occur soon after inflation through the
oscillation of the inflaton field at the bottom of the quadratic potential,
where the field decays into particles and reheats the universe. Such matter
phase can also be realised by any types of scalar field Kuroyanagi:2013ns ;
DEramo:2019tit ; Kuroyanagi:2020sfw  including the curvaton field
Nakayama:2009ce  and moduli particles in string theory Durrer:2011bi . The GW
modes which reenter the horizon during the matter-dominated phase have an
$f^{n_{t,\rm{prim}}-2}$ dependence. If we assume a perturbative decay of the
scalar field into fermionic particles, the bend of the spectrum can be
described by the following fitting function Kuroyanagi:2014nba :

$\frac{\Omega_{\rm GW}(f)}{\Omega_{\rm
GW,0}}=(1-0.22x_{R}^{1.5}+0.65x_{R}^{2})^{-1}\,,$ (29)

where $x_{R}\equiv f/f_{R}$ and $\Omega_{\rm GW,0}$ is the GW spectrum
obtained assuming a radiation-dominated era. The transition frequency $f_{R}$
corresponds to the energy scale at the end of reheating, more precisely
expressed akin to eq. (28) through substitution of $T_{*}$ with the
temperature $T_{R}$ of the thermal bath at the end of reheating.

Early kination phase.

In some scenarios of the early universe, the radiation-dominated epoch is
followed by a kination epoch, in which the kinetic energy of a scalar field
dominates the energy density of the universe. Examples of this type of model
are quintessential inflation Giovannini:1998bp ; Peebles:1998qn ;
Giovannini:1999bh ; Giovannini:1999qj ; Tashiro:2003qp ; Giovannini:2008tm ;
Ahmad:2019jbm  and particle-physics motivated scenarios Gouttenoire:2021jhk .
During the kination epoch, the energy density of the scalar field scales as
$\rho_{\phi}\propto a^{-6}$ and gives a Hubble expansion rate $H\propto
a^{-3}$, which results in the frequency dependence $\propto
f^{n_{t,\rm{prim}}+1}$. In this case, the SGWB spectrum is given by

$\frac{\Omega_{\rm GW}(f)}{\Omega_{\rm GW,0}}=(1-0.5x_{\rm
kin}^{2/3}+1.27\,x_{\rm kin})\,,$ (30)

where $x_{\rm kin}\equiv f/f_{\rm kin}$ and we have assumed an instantaneous
transition Figueroa:2019paj ; Duval:2024jsg . The transition frequency $f_{\rm
kin}$ is also expressed analogously to eq. (28) by replacing $T_{*}$ with
$T_{\rm kin}$, which is the temperature of the universe at the transition from
the kination to the radiation-dominated epoch. Note that the detailed shape of
the transition depends on the mechanism to end the kination phase. The second
term of Eq. (30) controls the details of the transition and should be modified
depending on the transition model.

####  2.8.2 Varying degrees of freedom

When we consider a contribution of one particle with mass $m$ and degrees of
freedom $g$, the feature can be fitted by a hyperbolic tangent function
Caldwell:2018giq :

$\Omega_{\rm GW}(f)=\Omega_{\rm GW,0}(f)F(f,g,m)\,,$ (31)

where

$F(f,g,m)=\frac{1-\epsilon(g)\tanh[\ln f/f_{0}(m)]}{1+\epsilon(g)}\,,$ (32)

where $\epsilon(g)=(1-\Delta)/(1+\Delta)$ and $\Delta\simeq(1+g/g_{\rm
SM})^{-1/3}$. The characteristic frequency is given by the particle mass as
$2\pi f_{0}(m)=H_{*}a_{*}/a_{0}|_{T\simeq m/b}$ with $b=2.2/\Delta$, which has
been determined empirically. Here $H_{*}$ and $a_{*}$ correspond to the Hubble
rate and the scale factor when the particle became non-relativistic and can be
expressed as $H(T_{*})=\sqrt{\pi^{2}g_{*}/90}~{}T_{*}^{2}/M_{\rm Pl}$ and
$a(T_{*})=(11g_{*}/43)^{-1/3}T_{0}/T_{*}$ with $g_{*}=g_{\rm SM}+g$, where
$M_{\rm Pl}$ is the reduced Planck mass.121212Recently, the detectability of
the impact of a hypothetical smooth crossover in the early universe beyond the
Standard Model of particle physics on the the scalar-induced gravitational
wave was reported in ref. Escriva:2024ivo .

##  3 Template-based reconstruction with the SGWBinner pipeline

This section summarizes the data analysis methods we employ in this work. Most
of the analyses rely on the SGWBinner code Caprini:2019pxz ; Flauger:2020qyi ,
which, for this work, has been modified to perform template-based analyses.
After briefly recalling some features of the LISA SGWB data analysis, we
discuss the key ingredients of the algorithm and the updates on the
astrophysical foregrounds, compared to ref. Flauger:2020qyi .

LISA will provide three time-domain data streams $d_{i}$, with $i$ running
over the LISA channels, which we divide into segments of individual duration
$\tau$. Then the frequency domain data $\tilde{d}_{i}(f)$ are defined via
Fourier transform as

$\tilde{d}_{i}(f)=\int_{-\tau/2}^{\tau/2}d_{i}(t)\;\textrm{e}^{-2\pi
ift}\;\textrm{d}t\;.$ (33)

where for simplicity the central time of that segment has been set to zero. In
the following, we assume that appropriate methods, _e.g.,_ some procedures to
be integrated within the LISA global fit scheme Cornish:2005qw ;
Vallisneri:2008ye ; MockLISADataChallengeTaskForce:2009wir ;
Littenberg:2023xpl , remove all transients, including loud deterministic
signals and glitches, from the data stream, leaving ‘clean data’ consisting
only of stochastic components. Under this assumption, we express the data as a
superposition of some noises $\tilde{n}^{\nu}_{i}$ and signals
$\tilde{s}_{i}^{\sigma}$ as

$\tilde{d}_{i}(f)=\sum_{\nu}\tilde{n}^{\nu}_{i}(f)+\sum_{\sigma}\tilde{s}^{\sigma}_{i}(f)\;,$
(34)

with $\nu$, $\sigma$ running respectively over the different noise and signal
components. We make the hypothesis that each noise and signal component is
independent of the others and characterised solely by specific statistical
properties. Additionally, our analysis is limited to signals that are
Gaussian, isotropic, stationary, and non-chiral.131313In principle, the real
signal might violate all these assumptions. For the impact of anisotropies and
non-Gaussianities see, _e.g.,_ refs. Bartolo:2019oiq ; Bartolo:2019yeu ;
Contaldi:2020rht ; LISACosmologyWorkingGroup:2022kbp ;
Dimastrogiovanni:2022eir ; Malhotra:2022ply ; Schulze:2023ich ;
Mentasti:2023icu ; Mentasti:2023uyi ; Perna:2024ehx . While most early
Universe mechanisms predict stationary signals, the presence of anisotropies,
projected in the data by a time-varying and sky-dependent response function,
will induce time modulations in the measurements. This effect is well-known to
be present, _e.g.,_ for the astrophysical foreground (see discussion around
eq. 40) due to the incoherent superposition of signals from compact binaries
in our galaxy Adams:2013qma . While including this effect in the analysis
might ease the separation of the different components, we ignore it in the
present work. Finally, the problem of detecting chirality with a planar
interferometer, like LISA, is highly non-trivial. By construction, planar
interferometers are insensitive to chirality Seto:2007tn ; Seto:2008sr ;
Smith:2016jqs . However, cross-correlating the measurements of different, and
non-coplanar, detectors Seto:2007tn ; Seto:2008sr ; Crowder:2012ik ;
Orlando:2020oko ; ValbusaDallArmi:2023ydl , or using the dipole induced by the
detector motion with respect to the SGWB frame Seto:2006hf ; Seto:2006dz ;
Domcke:2019zls , might help to overcome this limitation. Stationarity and
Gaussianity are assumed to hold for noise, too.141414As for the signal, these
hypotheses might be violated by real data. Transients (_e.g.,_ glitches) and
other effects, such as modulations of the noise due to instrumental component
degradation, might induce non-stationarities and non-Gaussianities in the
noise. While transients will be systematically modeled and removed from the
data stream (see, _e.g.,_ , refs. Robson:2018jly ; Baghi:2021tfd ), the real
analysis will consistently keep track of long-term modulations in the noise
model. Thus, assuming both components to have zero mean, we obtain

$\langle\tilde{d}_{i}(f)\rangle=0\;,\qquad\langle\tilde{d}_{i}(f)\tilde{d}^{*}_{j}(f^{\prime})\rangle=\frac{\delta(f-f^{\prime})}{2}\left[\sum_{\nu}P^{\nu}_{N,ij}(f)+\sum_{\sigma}P^{\sigma}_{S,ij}(f)\right]\;,$
(35)

where $P^{\nu}_{N,ij}$ and $P^{\sigma}_{S,ij}$ denote the noise and signal
power spectra respectively and the brackets denote an ensemble
average.151515The Dirac delta in frequency arises from stationarity, and, in
reality, it would only be an exact Dirac delta in the limit of infinite
observation time $T_{\rm obs}$. For finite observation time, a sinc$[T_{\rm
obs}(f-f^{\prime})]$ function appears in eq. 37. For simplicity, we restrict
ourselves to the limit where $T_{\rm obs}f\ll 1$, where the sinc can be
replaced by a Dirac delta. Notice that, in reality, each noise and signal
component is a combination of some ‘physical’ spectrum and a response (or
transfer) function, which projects the spectrum onto the data stream. In the
following, we will respectively denote with $T^{\nu}_{ij,lk}$ and
$\mathcal{R}_{ij}$ the transfer functions for the noise and signal components
(see refs. Flauger:2020qyi ; Hartwig:2021mzw ; Nam:2022rqg  for details).
While each noise component describes a different physical effect that
propagates differently through Time-Delay Interferometry (TDI), leading to a
different transfer function, all the signal components are GW signals, thus
sharing a common response function.

###  3.1 TDI, signal, and noise description

The laser frequency contribution LISA_performance  is the dominant source of
instrumental noise for LISA. To suppress this large noise component and to
allow any GW detection, Time-Delay Interferometry (TDI) is employed
Tinto:2020fcc . TDI is a post-processing technique that, combining
measurements performed at different times, produces synthesized data streams
representing laser-noise-free virtual interferometers. In the following, we
denote with $\eta_{ij}(t)$ the phase measurement performed in spacecraft $i$
at time $t$ of a signal emitted from spacecraft $j$ at time $t-L_{ij}$, where
$L_{ij}$ is the distance between the two spacecrafts. Moreover, we denote with
$D_{ij}$ the delay operator defined as $D_{ij}x(t)\equiv x(t-L_{ij})$. In
practice, TDI consists in defining variables as a linear combination of
single-link measurements and delay operators. As shown in refs. Armstrong_1999
; Prince:2002hp ; Shaddock:2003bc ; Shaddock:2003dj ; Tinto:2003vj ;
Vallisneri:2005ji ; Muratore:2020mdf ; Muratore:2021uqj , different TDI
combinations lead to laser-noise suppression, each with distinct sensitivities
to GW signals and instrumental noise. The most common choice for LISA data
analysis is to use three Michelson-like variables161616In this work, we
consider “first-generation” TDI variables only. These variables achieve laser-
noise cancellation in a scenario that respects our working assumptions. More
realistic investigations, which involve, _e.g.,_ , time-evolving unequal arms,
would require “second-generation” TDI variables Vallisneri:2005ji ;
Muratore:2020mdf ; Muratore:2021uqj ; Hartwig:2021mzw ., denoted as X, Y, and
Z, with X defined as

${\rm
X}=(1-D_{13}D_{31})(\eta_{12}+D_{12}\eta_{21})+(D_{12}D_{21}-1)(\eta_{13}+D_{13}\eta_{31})\;,$
(36a) and Y and Z variables are obtained through cyclic permutations of the
indexes. The XYZ variables are often recombined into (quasi-)orthogonal
channels, typically referred to as A, E, and T defined as Prince:2002hp  ${\rm
A}=\frac{{\rm Z}-{\rm X}}{\sqrt{2}}\;,\qquad{\rm E}=\frac{{\rm X}-2{\rm
Y}+{\rm Z}}{\sqrt{6}}\;,\qquad{\rm T}=\frac{{\rm X}+{\rm Y}+{\rm
Z}}{\sqrt{3}}\;.$ (36b)

It is known that for configurations with equal arm lengths and the same noise
levels for all spacecrafts, the AET combination is perfectly diagonal.
Moreover, under these assumptions, the T channel strongly suppresses GW
signals compared to instrumental noise, effectively providing a noise monitor.
Since in this work, we do not violate these assumptions171717See, _e.g.,_ ,
ref. Hartwig:2021mzw  for a recent analysis for a non-equilateral geometry and
unequal noise in the three spacecraft using different sets of TDI variables.,
we focus on this particular TDI basis, significantly simplifying the
computations required, _e.g.,_ for likelihood evaluation. This feature makes
the AET channels particularly appealing for our analysis.

Assuming TDI suppresses laser noise, the main residual noise sources for LISA
(also known as secondary noises) are the Test Mass (TM) noise, representing
deviations from free-fall in the TM trajectories, and the Optical Metrology
System (OMS) noise, representing uncertainties in the determination of the TM
positions. The total noise power spectrum for any TDI variable can thus be
expressed as

$P_{N,ij}(f)\equiv\sum_{\nu}P^{\nu}_{N,ij}(f)=\left[T^{\textrm{TM}}_{ij,lk}(f)S^{\textrm{TM}}_{lk}(f)+T^{\textrm{OMS}}_{ij,lk}(f)S^{\textrm{OMS}}_{lk}(f)\right]\;,$
(37)

where $T^{\textrm{TM}}_{ij,lk}$ and $T^{\textrm{OMS}}_{ij,lk}$ are the TM and
OMS transfer functions, projecting the individual TM and OMS noise components
onto the TDI variable. $S^{\textrm{TM}}$ and $S^{\textrm{OMS}}$ are estimated
as LISA:2017pwj

$\displaystyle S^{\text{TM}}_{lk}(f)$
$\displaystyle=A_{lk}^{2}\;\left(1+\left(\frac{0.4\textrm{mHz}}{f}\right)^{2}\right)\left(1+\left(\frac{f}{8\textrm{mHz}}\right)^{4}\right)\left(\frac{1}{2\pi
fc}\right)^{2}\;\left(\frac{\mathrm{fm}^{2}}{\mathrm{s}^{3}}\right)\;,$ (38)
$\displaystyle S^{\text{OMS}}_{lk}(f)$
$\displaystyle=P_{lk}^{2}\;\left(1+\left(\frac{2\times
10^{-3}\textrm{Hz}}{f}\right)^{4}\right)\times\left(\frac{2\pi
f}{c}\right)^{2}\;\times\left(\frac{\mathrm{pm}^{2}}{\mathrm{Hz}}\right)\;.$
(39)

Under the assumption of equal arms and noise levels, we can simplify the above
coefficients as products of Kronecker $\delta$s and two constant parameters,
$A_{lk}^{2}=A^{2}\,\delta_{lk}$ and $P_{lk}^{2}=P^{2}\,\delta_{lk}$. It is
worth noting that, for the A and E channels, the TM noise dominates at lower
frequencies, while the OMS noise dominates at higher frequencies. On the
contrary, the OMS noise prevails across all frequencies for the T channel. In
the analyses performed in this paper, we apply Gaussian priors to $A$ and $P$,
which are centered on their nominal values, 3 and 15, and have a width of 20%.

Beyond cosmological signals, the SGWB signal in the LISA band will have two
guaranteed contributions from astrophysical sources. At low frequency, there
will be a stochastic contribution from many unresolved Compact Galactic
Binaries (CGBs) Nissanke:2012eh . At larger frequencies, there will be a
signal due to the incoherent superposition of Stellar Origin Black Hole
Binaries (SOBHBs) and Binary Neutron Stars (BNSs) Regimbau:2011rp ;
KAGRA:2021kbb  (see also Perigois:2020ymr ; Babak:2023lro ; Lehoucq:2023zlt ).
As a consequence, the signal part of eq. 35 can be expressed as

$P_{S,ij}(f)\equiv\sum_{\sigma}P^{\sigma}_{S,ij}(f)=\mathcal{R}_{ij}(f)\left[S_{\textrm{Gal}}(f)+S_{\textrm{Ext}}(f)+S_{\textrm{Cosmo}}(f)\right]\;,$
(40)

where $S_{\textrm{Gal}}(f)$ and $S_{\textrm{Ext}}(f)$ represent the galactic
and extragalactic foreground contribution respectively, and the last
contribution is given by all the signal templates listed in section 2 and
related to the underlying inflationary mechanism. The GW strain power spectral
density (PSD) is related to $\Omega_{\textrm{GW}}(f)$, energy density per
logarithmic frequency interval normalized by the critical density, defined in
eq. (1), through

$h^{2}\Omega_{\textrm{GW}}(f)\equiv\frac{4\pi^{2}f^{3}}{3H_{0}^{2}}S_{h}(f)\;.$
(41)

This provides a rescaling between $S$ and $\Omega$, which we employ also for
the noise.

The two foreground models we adopt in this paper are the state-of-the-art
spectral models. In particular, for the CGB we use the model of ref.
Karnesis:2021tsh :

$S_{\textrm{Gal}}(f)=A_{\textrm{Gal}}\left(\frac{f}{1\,\textrm{Hz}}\right)^{-\frac{7}{3}}\times
e^{-(f/f_{1})^{\alpha}}\times\frac{1}{2}\left[1+\tanh{\frac{f_{\textrm{knee}}-f}{f_{2}}}\right]\;,$
(42)

where the first factor comes from the superposition of many inspiraling
signals Phinney:2001di , the second factor is due to the loss of stochasticity
at higher frequencies, and the last factor, which produces a sharp cut-off in
the spectrum, models the complete removal of binaries emitting at sufficiently
large frequency. For what concerns the values of the parameters used for the
injection, we use the time-dependent parameterisation introduced in ref.
Karnesis:2021tsh :

$\displaystyle\log_{10}(f_{1})$ $\displaystyle=$ $\displaystyle
a_{1}\log_{10}(T_{\mathrm{o}bs})+b_{1}\,,$
$\displaystyle\log_{10}(f_{\textrm{knee}})$ $\displaystyle=$ $\displaystyle
a_{k}\log_{10}(T_{\mathrm{o}bs})+b_{k}\,,$ (43)

and set $T_{\mathrm{o}bs}=4$ years with $100\%$ duty cycle,
$a_{1}=-0.15;~{}b_{1}=-2.72;~{}a_{k}=-0.37;~{}b_{k}=-2.49$, together with
$A_{\textrm{Gal}}=1.15\cdot 10^{-44};~{}\alpha=1.56;f_{2}=6.7\times
10^{-4}$Hz. In principle, all these should be measured together with the
parameters of the SGWB of cosmological origin. In practice, we restrict our
analysis by varying only the amplitude parameter $A_{\textrm{Gal}}$. In
reality, since we work in $\Omega_{\rm{GW}}$ units, the parameter we use in
the analysis is $\log_{10}(h^{2}\Omega_{\rm
Gal})\equiv\log_{10}[4\pi^{2}A_{\textrm{Gal}}/(3H_{0}^{2})]$. Finally, we
impose a Gaussian prior on $\log_{10}(h^{2}\Omega_{\rm Gal})$ with central
value $\simeq-7.8412$ and standard deviation $\simeq 0.21$.

For the extragalactic contribution, the background signal can be adequately
described by a power-law model with a fixed value for the tilt:

$h^{2}\Omega_{\textrm{Ext}}=10^{\log_{10}(h^{2}\Omega_{\rm
Ext})}\left(\frac{f}{0.001{\mathrm{H}z}}\right)^{2/3}\;.$ (44)

As for the galactic template in eq. 42, the tilt comes from the superposition
of many signals in the inspiral phase Phinney:2001di . While the subtraction
of sufficiently loud events might lead to deviations from this behavior, this
effect is expected to be small Babak:2023lro . For this reason, in the present
work, we assume the parameterisation in eq. 44, controlled by the amplitude
parameter only, to suffice. Recent observations by the LVK collaboration
KAGRA:2021duu  suggest that the SGWB due to SOBHB and BNS should have
$\Omega_{\textrm{SOBHB+BNS}}(25\textrm{Hz})=7.2^{+3.3}_{-2.3}\times 10^{-10}$
which, rescaled at LISA frequency $(f_{*}=10^{-3}\textrm{Hz})$ implies
$\log_{10}(h^{2}\Omega_{\rm Ext})\simeq-12.38$. In our analyses, we impose a
Gaussian prior on this parameter, centered around such a value, and with a
standard deviation equal to $\simeq 0.17$.

###  3.2 Data generation, likelihood, and Fisher matrix

LISA is scheduled to function for at least 4.5 years up to a maximum of 10
years. Operations such as antenna repointing introduce (periodic) data gaps in
the schedule, resulting in a duty cycle of about $82\,\%$ Colpi:2024xhw . We
assume an intermediate scenario, setting the SGWBinner code to work with
$N_{d}=126$ data segments of duration $\tau=11.4$ days each. This sums up to
$T_{\rm obs}=4$ effective years of data. Then, we perform a Fourier transform
in each of the segments to get the frequency domain data
$\tilde{d}^{s}_{i}(f_{\textmd{k}})$, where $s$ indexes segments, k indexes
frequencies within the detector range and $i$ indexes the TDI channels. As
mentioned above, we assume a segment duration of approximately 11.5 days,
corresponding to $\Delta f=N_{d}/T_{\rm obs}\sim 10^{-6}$ Hz. Under these
assumptions, we generate $N_{d}$ Gaussian realisations for the signal and all
noise components, with zero mean and variances defined by their respective
power spectral densities. To lower the numerical complexity of the problem we
perform two operations. Firstly, we define a new set
$\bar{D}^{\textmd{k}}_{ij}\equiv\sum_{s=1}^{N_{d}}\tilde{d}^{s}_{i}(f_{\textmd{k}})\tilde{d}^{s}_{j}(f_{\textmd{k}})/N_{d}$,
by averaging over segments. Then, we down-sample these data using the coarse-
graining procedure introduced in Caprini:2019pxz ; Flauger:2020qyi . By
applying these techniques, we obtain a new data set $D^{k}_{ij}$, where $k$
now indexes a sparser set of frequencies $f^{k}_{ij}$. These frequencies are
weighted according to $w^{k}_{ij}$, which corresponds to the number of points
averaged over during the coarse-graining procedure. The down-sampled data set
retains statistical properties similar to those of
$\bar{D}^{\textmd{k}}_{ij}$, while being computationally more manageable.

The full likelihood employed in our analyses reads

$\ln\mathcal{L}(\vec{\theta})=\frac{1}{3}\ln\mathcal{L}_{\mathrm{G}}(\vec{\theta}|D^{k}_{ij})+\frac{2}{3}\ln\mathcal{L}_{\mathrm{LN}}(\vec{\theta}|D^{k}_{ij})\;,$
(45)

given by the sum of a Gaussian

$\ln\mathcal{L}_{\mathrm{G}}(\vec{\theta}|D^{k}_{ij})=-\frac{N_{d}}{2}\sum_{k}\sum_{i,j}w^{k}_{ij}\left[1-D^{k}_{ij}/D^{\mathrm{T}h}_{ij}(f^{k}_{ij},\vec{\theta})\right]^{2}\;,$
(46)

and of a log-normal component

$\ln\mathcal{L}_{\mathrm{LN}}(\vec{\theta}|D^{k}_{ij})=-\frac{N_{d}}{2}\sum_{k}\sum_{i,j}w^{k}_{ij}\ln^{2}\left[D^{\mathrm{T}h}_{ij}(f^{k}_{ij},\vec{\theta})/D^{k}_{ij}\right]\;.$
(47)

The latter is included to take into account the mild non-Gaussianity
introduced by the data generation, and to avoid biased results Bond:1998qg ;
Sievers:2002tq ; WMAP:2003pyh ; Hamimeche:2008ai . Here
$\vec{\theta}=\\{\vec{\theta}_{s},\vec{\theta}_{n}\\}$ is the vector of signal
and noise parameters

$\vec{\theta}_{n}=\\{A,P\\}\;,\qquad\quad\vec{\theta}_{s}=\\{\vec{\theta}_{\textrm{fg}},\vec{\theta}_{\textrm{cosmo}}\\}\;,$
(48)

with

$\displaystyle\vec{\theta}_{\textrm{fg}}=\\{h^{2}\Omega_{\rm
Gal},h^{2}\Omega_{\rm Ext}\\}\;,$
$\displaystyle\vec{\theta}_{\textrm{cosmo}}\qquad\textrm{(template
dependent)}\;,$ (49)

while $D^{\mathrm{T}h}_{ij}(f_{k},\vec{\theta})$ is the theoretical model for
the data (containing both signal and noise). As explained above, under the
assumptions discussed in the previous section, the AET TDI basis is diagonal,
so that $D^{k}_{ij}\neq 0$ only for $i=j$, which further simplifies the
analysis. The priors for the signal, foreground, and noise parameters are
added to eq. 45 to get the posterior distribution. While the noise and
foreground priors are discussed in section 3.1, the priors for the signal
parameters of each template are set according to the discussion in section 2.
Finally, we sample the parameter space using the nested sampling algorithm
implemented in Polychord Handley:2015vkr ; Handley:2015fda , via its Cobaya
Torrado:2020dgo  interface, and visualise results using GetDist Lewis:2019xzd
.

We conclude this section by recalling the main ingredients of the Fisher
Information Matrix (FIM) formalism. The log-likelihood for Gaussian and zero
mean data $\tilde{d}_{i}^{k}$, with $i$ running over TDI channels and $k$
running over frequency bins $f_{k}$, described only by their variance, say
$C_{ij}(f_{k},\vec{\theta})$, can be written as

$-\ln\mathcal{L}(\tilde{d}^{k}_{i}|\vec{\theta})\propto\sum_{k}\left\\{\ln\left\\{\det[C_{ij}(f_{k},\vec{\theta})]\right\\}+\tilde{d}^{k}_{i}\,C^{-1}_{ij}(f_{k},\vec{\theta})\,\tilde{d}_{j}^{k*}\right\\}\;,$
(50)

which is also known as Whittle likelihood. The FIM $F_{\alpha\beta}$,
representing the information on the model parameters, is defined as

$F_{\alpha\beta}\equiv-\left.\frac{\partial^{2}\ln\mathcal{L}}{\partial\theta^{\alpha}\partial\theta^{\beta}}\right|_{\vec{\theta}=\vec{\theta}_{0}}=\sum_{k}\textrm{Tr}\left[C^{-1}\frac{\partial
C}{\partial\theta^{\alpha}}C^{-1}\frac{\partial
C}{\partial\theta^{\beta}}\right]_{\vec{\theta}=\vec{\theta}_{0}}\;,$ (51)

where, $\vec{\theta}_{0}$ represent the best-fit parameter(s), determined by
solving

$\left.\frac{\partial\ln\mathcal{L}}{\partial\theta^{\alpha}}\right|_{\vec{\theta}=\vec{\theta}_{0}}\propto\sum_{k}\frac{\partial
C_{ji}}{\partial\theta^{\alpha}}\left[C^{-1}_{ij}-C^{-1}_{im}\;\tilde{d}_{m}^{k}\tilde{d}_{n}^{k*}\;C^{-1}_{nj}\right]=0\;,$
(52)

under the assumption
$C_{ij}(f_{k},\theta_{0})=\tilde{d}^{k}_{i}\tilde{d}_{j}^{k*}$. In practice,
the discrete sum over finite frequencies can be replaced with a continuous
integral over the frequency range as

$F_{\alpha\beta}\equiv T_{\rm obs}\sum_{i\in\\{{\rm
AET}\\}}\int_{f_{\mathrm{min}}}^{f_{\mathrm{max}}}\frac{\partial\ln
C_{ii}}{\partial\theta^{\alpha}}\frac{\partial\ln
C_{ii}}{\partial\theta^{\beta}}\,\textrm{d}f\;,$ (53)

where $f_{\mathrm{min}}$, $f_{\mathrm{max}}$ are the minimal and maximal
frequencies measured by the detector, which we assume to be
$f_{\mathrm{min}}=3\times 10^{-5}$ Hz and $f_{\mathrm{max}}=0.5$ Hz
Colpi:2024xhw , while we remind that $T_{\rm obs}$ is the total observation
time and that we have exploited the fact that the AET basis is diagonal. If
non-trivial (log-)priors are included in the analysis, their derivatives
should be consistently added to eq. 53 to get the full FIM. Finally, the
covariance matrix $C_{\alpha\beta}$, which provides estimates on the
determination and (on the covariance) of the model parameters, is computed by
inverting the FIM. When we provide the errors on a given parameter, we
marginalize over all the other ones, which amounts in taking the square root
of the diagonal element of the inverse of the correlation matrix for that
parameter as its error.

Given its computational efficiency, in the following section we employ the FIM
approach to scan the parameter space of the templates introduced in section 2,
and to assess the prospect of reconstructing the template parameters with some
level of accuracy. It is worth stressing 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. Typically, a rule of
thumb to assess the quality of the FIM approximation is through the Signal-to-
Noise ratio (SNR), defined as Romano:2016dpx

$\textrm{SNR}\equiv\sqrt{T_{\rm obs}\;\sum_{i\in\\{{\rm AET}\\}}\int_{f_{\rm
min}}^{f_{\rm
max}}\left(\frac{S_{i,\mathrm{GW}}}{S_{i,\mathrm{N}}}\right)^{2}\;\textrm{d}f}\;.$
(54)

The SNR scales linearly with the signal amplitude and with the square root of
the observation time. Given the limitations of the FIM, as already mentioned,
we will test the validity of the FIM results by directly sampling the
likelihood using Nested sampling.

##  4 Reconstruction forecasts with the SGWBinner

In this section we apply the methodology outlined in section 3 to forecast the
capabilities of LISA to reconstruct Gaussian, isotropic SGWB signals
characterised by the templates discussed in section 2. We focus our analysis
on the frequency structure of the templates and the parameters characterising
them. In section 5 we then discuss the physical consequences of our findings
for inflationary models.

The results we present are based on several assumptions (some of them
anticipated in section 1). We assume that any potential discrepancy between
instrumental noise and the noise model, as well as between the SGWB signal and
the chosen template, introduces systematic errors below the level of
statistical uncertainties.

The analysis of each template follows the same rationale. We initially employ
a FIM formalism to forecast the reconstruction errors on parameters of each
inflationary template, taking into account the reconstruction uncertainties of
the instrumental noise and foregrounds. In this way, for every inflationary
template parameter $\theta_{\textrm{cosmo},i}$, we compute its Fisher
reconstruction error marginalized over all the reconstruction uncertainties on
the remaining signal and noise parameters, namely $\vec{\theta}_{n}$,
$\vec{\theta}_{\textrm{fg}}$ and $\theta_{\textrm{cosmo},j}$ with every $j\neq
i$ (see eq. 48). We then present our Fisher results in two-dimensional color
maps displaying how each error changes when varying the injected values of a
pair of inflationary template parameters and leaving the others parameters at
some given fixed injected values (see, _e.g.,_ fig. 2). To facilitate the
interpretation we also highlight specific contour lines for the errors $1\%$
and $30\%$, or $0.01$ and $0.3$, depending on whether the relative or absolute
error is reported in the maps. To investigate the impact of the astrophysical
foregrounds, such contour lines are determined both in the case where the
foreground parameters are known a priori and in the case they are
reconstructed together with the inflationary and noise template parameters.

The Fisher analysis has a drawback since it assumes the likelihood to be
Gaussian in the inflationary signal, foregrounds, and noise parameters. Due to
numerical precision limitations, it also struggles in dealing with strong
multidimensional correlations among parameters differing by many orders of
magnitude. In our analysis we test the Fisher approximation on signals
following our templates. In most of our tests, we find that the Fisher
approximation reliably estimates the order of magnitude of the errors when the
signals have ${\rm SNR}\gtrsim 10$ and the estimated errors are much smaller
than the injected values. Therefore, the results of our Fisher analysis must
be interpreted with caution when these two conditions are not met.

We then pursue the analysis of each template focusing on benchmarks. For each
benchmark, we perform a complementary Bayesian template-based analysis by
means of the SGWBinner Caprini:2019pxz ; Flauger:2020qyi , and we present
triangle plots showing 1D and 2D posterior distributions for the parameters of
the primordial signal obtained by directly sampling the likelihood using
nested sampling (see _e.g.,_ fig. 3). For these benchmarks, in order to
elucidate the Fisher analysis accuracy, we also plot the 1D and 2D contours
from the Fisher results. This helps capture some features that the Fisher
analysis does not catch. For our purposes, we consider one simulation per
benchmark, but it is worth noting that finer details of the results can
exhibit some realisation dependencies. In addition to the triangle plots, we
use the samples from the nested sampling analyses to plot the functional
posterior distribution of noise, foregrounds and primordial signal. The
latter, often called predictive posterior distribution, provides complementary
insights into the ability of LISA to reconstruct the benchmark signals. For
these plots, we use the publicly available fgivenx code fgivenx ;
Handley:2019fll .

The anisotropic nature of the galactic foreground gives rise to a time-
dependent modulation of the signal, as the detector moves through space
Allen:1996gp ; Cornish:2003tz . The improvement is however limited
Mentasti:2023uyi , and our analysis does not leverage this feature, opting
instead for considering the average signal integrated over the mission
duration. Therefore, our work may be suboptimal in this respect. For this
reason, we present comparisons between the reconstruction forecasts with and
without astrophysical components. In the latter case, the results can serve as
reference for the separation of the astrophysical and cosmological components
of the SGWB. In general, we expect that the presence of foregrounds impacts
the reconstruction in two ways: foregrounds can introduce degeneracies among
parameters, and the galactic (extragalactic) foreground can cover the
primordial signal in a large (small) fraction of the low-frequency (high-
frequency) sensitivity region. We will meet examples of these phenomena in our
analysis below.

###  4.1 Forecasts for the power law template

Figure 2: Fisher forecast for the PL template. The color map shows the 68 % CL
errors on the amplitude (left) and spectral index (right) as a function of the
value of the injected parameters. SNR contour lines are plotted in white. The
pairs of dashed (dotted) lines mark the $\sigma=0.3$ ($\sigma=0.01$) contours,
respectively in the absence (black) and in the presence (purple) of
foregrounds. The black and gray crosses display the benchmarks PL-BNK_1 and
PL-BNK_2, respectively. We also highlight the line $n_{T}=2/3$, which
corresponds to the spectral index of extra-galactic foregrounds, with a black
dashed line.

Figure 3:  _Left panel_ : 1D and 2D posterior distributions derived from the
PL template reconstruction of the benchmark PL-BNK_1
($\\{\log_{10}(h^{2}\,\Omega_{*}),n_{T}\\}=\\{-12.5,\,2.085\\}$) in the
presence of galactic and extragalactic foregrounds. The posteriors of
foregrounds and noise parameters are omitted for clarity. In the corner plot,
the gold and red dots, and corresponding vertical lines, show the injected
parameters and their reconstructed mean values. The 68% C.L. (95% C.L.)
reconstruction region is displayed in dark (light) gray. The equivalent
regions and 1D posteriors obtained with the Fisher approximation are in blue.
The top-right inset visualises the injected and reconstructed signals, with 68
and 95 % CL error bands. The error bands on the galactic foreground and
instrumental noise are too small to be visible. The LISA PLS is plotted in
solid black. _Right panel_ : Like the left panel but for PL-BNK_2
($\\{\log_{10}(h^{2}\,\Omega_{*}),n_{T}\\}=\\{-11,\,0.77\\}$). In the inset
plots of both panels, the error bands on the galactic foreground and
instrumental noise are too small to be visible.

In fig. 2 and fig. 3 we present the forecasts for the PL template discussed in
section 2.1. The color maps show that the errors on both the amplitude and the
tilt of the spectrum decrease as the amplitude $\Omega_{*}$ increases, as the
power law becomes steeper, __i.e.,__ $\lvert n_{T}\rvert$ increases, and, more
in general, as the SNR of the primordial signal becomes larger. Figure 2 shows
that the SNR contours of the primordial signal generally (although not
perfectly) follows the reconstruction error contours.181818The SNR is
evaluated with respect to the nominal LISA sensitivity, which is the one we
inject, and takes into account the existence of three TDI channels and the
foregrounds as an additional source of nuisance. Several studies have adopted
the criterion SNR $\gtrsim 10$ as a proxy for the condition of SGWB
detectability and reconstruction Caprini:2019pxz ; Flauger:2020qyi . Comparing
the black and purple dashed/dotted lines, we learn how the presence of
foregrounds degrades the measurements of the two signal parameters. For
example, for a flat signal (__i.e.,__ $n_{T}=0$), an accuracy
$\sigma\simeq\mathcal{O}(0.01)$ on the logarithm of the amplitude requires
$h^{2}\Omega_{*}\simeq 10^{-12}$ without foregrounds, and is only slightly
larger when foregrounds are included. Achieving the same level of accuracy on
the tilt requires slightly larger values for the signal amplitude. Notice a
peculiar behaviour along the line at $n_{T}=2/3$, where the primordial SGWB is
degenerate with the foreground due to the extragalactic compact binaries. The
possibility to separate a primordial signal with $n_{T}\simeq 2/3$ from
extragalactic foregrounds crucially depends on our prior knowledge about the
amplitude of the latter, which we have in practice implemented through a
Gaussian prior, as explained in the previous section.

Benchmarks. — We consider two benchmarks: PL-BNK_1 and PL-BNK_2. In the
former, the PL parameters are set as
$\\{\log_{10}(h^{2}\,\Omega_{*}),n_{T}\\}=\\{-12.5,\,\,2.085\\}$, while in the
latter they are set as
$\\{\log_{10}(h^{2}\,\Omega_{*}),n_{T}\\}=\\{-11,\,0.77\\}$. Both benchmarks
can be produced within the axion inflation scenario, while the first of them
is consistent with models of inflation with broken space diffeomorphisms, but
does not respect the so called Higuchi bound for massive gravitons
Higuchi:1986py  (more on this in section 5.1).

We run the PL-template-based SGWBinner analysis on our two benchmarks, and
display the obtained 1D and 2D posteriors in Figure 3. As the corner plots in
the figure show, the injected values of both benchmarks are reconstructed well
within the 68 % CL contours (the foreground and noise reconstruction are
omitted for clarity). In each panel, the inset plot highlights the injected
and reconstructed benchmark signal, noise the foregrounds (injected and
modelled as explained in section 3.1) with their 68 and 95 % CL error bands.
For the galactic foreground, the reconstruction is very accurate, with the
reconstructed amplitude within the 68 % CL error band (recall that we vary
only the amplitude in our analysis, keeping the spectral shape of foregrounds
fixed) while the error bands on the extragalactic foreground are larger, but
still within the 68 % CL error band. For the PL-BNK_2, the degeneracy between
the PL signal and the extragalactic foreground leads to a slightly less
accurate reconstruction. For our two injected signals, and in particular for
PL-BNK_1, we notice a very good agreement with the Fisher analysis, which
captures very well the shape of the posterior distribution of the tilt and
amplitude of the power law. We conclude that for a PL template with a
sufficiently large injected signal with respect to the LISA sensitivity, the
signal reconstruction degrades in the presence of foregrounds, but it is still
very accurate. As expected, the impact of the foregrounds is less pronounced
when the amplitude (or the tilt) of the PL is large.

###  4.2 Forecasts for the log-normal template

We now investigate the LN template given in eq. 12. Figure 4 displays the
Fisher-approximation reconstruction errors obtained in the parameter plane
with peak frequency $f_{*}=1\,$mHz (left column) and the parameter plane with
bump width $\log_{10}\rho=-0.37$ (right column). When fixing the width of the
bump, the impact of the foregrounds on the reconstruction is less important
with respect to the case where we fix the peak frequency, at least for what
concerns the accuracy on the error on the amplitude. Overall, an accuracy of
$\sigma\sim 0.1$ can be obtained with a signal characterised by SNR
$\gtrsim\,10$ in the absence of foregrounds. When foregrounds are included, to
reach the same accuracy, an order of magnitude larger amplitude is required.
The minimum error on $\Omega_{*}$, $f_{*}$, and $\rho$ is reached for
$10^{-3}\lesssim f_{*}/{\rm Hz}\lesssim 10^{-2}$, where the peak of the signal
is in the best sensitivity region of LISA. When we fix the frequency peak of
the bump (say at $1$ mHz) and we vary the width and the amplitude, the error
on the latter is at percent level for $h^{2}\Omega_{*}\gtrsim 10^{-12}$, while
for $f_{*}$ the accuracy improves for signals with smaller width (at fixed
amplitude). Finally, when we fix the frequency peak, a better accuracy on the
width requires larger amplitude signals with smaller width. Overall, the
constraints on all parameters degrade when the signal peak approaches the
borders of the LISA frequency band.

Figure 4: Fisher forecasts for the LN template. The panels show the 68 % CL
reconstruction error on $h^{2}\Omega_{*}$ (top panels), $f_{*}$ (central
panels) and $\rho$ (bottom panels) as a function of the injected values of
$f_{*}$, $h^{2}\Omega_{*}$ and $\rho$ set as specified in the axes and the
title of each panel. SNR countour lines are plotted in white. The pairs of
dashed (dotted) lines mark the $\sigma=0.3$ ($\sigma=0.01$) contours,
respectively in the absence (black) and in the presence (purple) of
foregrounds. The gray and black cross display the benchmarks LN-BNK_1 and LN-
BNK_2 respectively.

Figure 5:  1D and 2D posterior distributions derived from the LN template
reconstruction of the benchmarks LN-BNK_1 (left panel) and LN-BNK_2 (right
panel). Lines’ styles and color codes are as in fig. 3. LN-BNK_1 and LN-BNK_2
are defined as in eq. 12 with
$\\{\log_{10}(h^{2}\,\Omega_{*}),\,\log_{10}\rho,\,f_{*}/{\rm mHz}\\}$ set to
$\\{-9.9,\,-0.13,\,11\\}$ and $\\{-11.8,\,-0.37,\,1\\}$, respectively.

Benchmarks. — We consider the benchmarks LN-BNK_1 and LN-BNK_2, which are
defined as in eq. 12 with
$\\{\log_{10}(h^{2}\,\Omega_{*}),\,\log_{10}\rho,\,f_{*}\,{\rm[Hz]}\\}$ set at
$\\{-9.9,\,-0.13,\,10^{-3}\\}$ and $\\{-11.8,\,-0.37,\,10^{-3}\\}$,
respectively. These parameter choices are based on the setups detailed in
section 2.2.

In fig. 5 we plot the 1D and 2D posterior distributions of the LN parameters
obtained by the template-based SGWBinner analysis run on LN-BNK_1 (left panel)
and LN-BNK_2 (right panel). For the LN-BNK_1, the reconstructed mean values
for the amplitude, peak frequency, and width parameters are within 68 % CL .
The marginalized posterior distribution for all the parameters appears
Gaussian with the reconstructed parameters all within 1% accuracy, and in good
agreement with the Fisher analysis presented above. Therefore signal, noise,
and foregrounds are all well reconstructed for this specific benchmark. Only
for the extragalactic foreground the reconstruction is less accurate, as it is
partly obscured by the loud signal and might be prior-dominated. The posterior
distributions of the signal parameters agree very well with the Fisher
analysis; we only notice an almost negligible shift of the mean values
compared to the injected signal.

The situation is different for LN-BNK_2 (right plot), which is characterised
by an amplitude smaller than the foregrounds, but still large enough for the
bump to be detected. However, the bump being below the foreground level, its
frequency shape reconstruction has wide error bars. This can also be noticed
by the fact that the best fits are different from the injected signal, but
still fall within the 95 % CL contours of the posterior distributions.
Constraints are definitely looser than those for LN-BNK_1, but the signal
parameters can still be constrained. Despite the high SNR ($\sim 50$), we
observe a difference between the posterior distributions obtained through the
nested sampling and the contours obtained using the Fisher approach. Such a
difference may arise both from the non-Gaussian posteriors and from a certain
level of realisation dependence. As a result, the 68 % CL region tends to
shift towards amplitudes slightly smaller than the injected values, leading to
worsened constraints. Nevertheless, the constraints on parameters remain
nearly of the same order of magnitude as those obtained through the Fisher
approach, demonstrating its reliability in providing a reasonable good
estimate of error order of magnitudes. Interestingly, the foregrounds being
larger than the cosmological signal, they are reconstructed with very high
accuracy.

###  4.3 Forecasts for the broken power law template

We now discuss the BPL template defined in eq. (13). As the template depends
on numerous parameters, we only report on the analysis covering some
parameter-space slices crossing the benchmarks BPL-BNK_1 and BPL-BNK_2 for
brevity. The template parameters $\log_{10}(h^{2}\,\Omega_{*}),f_{*}/{\rm
mHz},\,n_{t,1},\,n_{t,2},\,\delta$ are defined to be equal to
$-10.5,1,\,4,0,\,1$ for BPL-BNK_1 and equal to $-9.3,1,\,2.65,-2.1,\,5.3$ for
BPL-BNK_2. These benchmarks are motivated by the theory scenarios discussed in
section 2.3.191919The BPL template also describes the SGWB that cosmological
first-order phase transitions produce in some regimes Caprini:2015zlo ;
Caprini:2019egz . Reference Caprini:2024hue  analyses the LISA BPL-template-
based reconstruction dedicated to such a cosmological SGWB source. That
analysis differs from the present one as theoretical knowledge of the first-
order phase transition dynamics permits fixing several BPL parameters a
priori. This suppresses many degeneracies that arise in the present BPL
reconstruction where we keep all template parameters a priori unknown to cover
the huge variety of inflationary mechanisms leading to a BPL signal. Since the
two scenarios lead to qualitatively different parameter choices – despite
being described by the same template – the discussion in this section will
differ from the other sections, and we will address each scenario
individually.

Figure 6: Fisher forecasts for the BPL template. The panels show the 68 % CL
reconstruction error on $h^{2}\Omega_{*}$ (top left), $f_{*}$ (top right),
$n_{t,1}$ (central left), $n_{t,2}$ (central right) and $\delta$ (bottom) as a
function of the injected values of $f_{*}$, $h^{2}\Omega_{*}$ specified in the
axes and $n_{t,1}$, $n_{t,2}$ and $\delta$ specified in the title of each
panel. SNR contour lines are plotted in white. The pairs of dashed (dotted)
lines mark the $\sigma=0.3$ ($\sigma=0.01$) contours, respectively in the
absence (magenta) and in the presence (red) of foregrounds. The green crosses
display the benchmarks BPL-BNK_1.

Figure 7:  1D and 2D posterior distributions derived from the BPL template
reconstruction of the benchmark BPL-BNK_1. Lines’ styles and color codes are
as in fig. 3. In addition, the purple contours display the 1D and 2D
posteriors recovered by assuming the injected parameter $\delta=1$ to be known
a priori. BPL-BNK_1 is defined as in eq. (13) with
$\log_{10}(h^{2}\,\Omega_{*}),f_{*}/{\rm mHz},\,n_{t,1},\,n_{t,2},\,\delta$
equal to $-10.5,1,\,4,0,\,1$, respectively.

Figure 8:  1D and 2D posterior distributions derived from the BPL template
reconstruction of the benchmark BPL-BNK_2. Lines’ styles and color codes are
as in fig. 3. BPL-BNK_2 is defined as in eq. (13) with
$\log_{10}(h^{2}\,\Omega_{*}),f_{*}/{\rm mHz},\,n_{t,1},\,n_{t,2},\,\delta$
equal to $-9.3,1,\,2.65,-2.1,\,5.3$, respectively.

BPL-BNK_1. Figure 6 presents the Fisher reconstruction forecast maps in the
case that the template parameters $f_{*}$ and $h^{2}\Omega_{*}$ are injected
with values varying in the intervals $[10^{-5},10^{-1}]$ and
$[10^{-14},10^{-9}]$, respectively, while the other template parameters are
set as for BPL-BNK_1. It turns out that, in this parameter-space slice,
requiring all parameters to be reconstructed with $\sigma\lesssim 0.3$
typically implies SNR $\gtrsim\\!\\!\,10^{3}$. An error of the order
$\sigma\sim 0.01$ on the parameters $\log_{10}(h^{2}\,\Omega_{*})$, $f_{*}$
and $n_{t,1}$ requires $\log_{10}(h^{2}\,\Omega_{*})\gtrsim-9$, while for
$n_{t,2}$ an amplitude $\log_{10}(h^{2}\,\Omega_{*})\gtrsim-11$ is sufficient.
We notice that, for a given signal amplitude, the SNR is independent of
$f_{*}$ for $f_{*}\lesssim 10^{-4}{\rm Hz}$. This is due to the fact that for
$f_{*}\lesssim 10^{-4}{\rm Hz}$ the part of the signal on the left of $f_{*}$
is outside the LISA sensitivity, whereas the one on the right is just a
constant. As soon as $f_{*}$ enters the LISA frequency band, the SNR starts
showing a dependence on both the amplitude and the other parameters. The
errors on the amplitude and peak frequency have their minima when the
amplitude has a peak in correspondence with the LISA maximum sensitivity. A
similar dependence is shown for the error on the low frequency tilt $n_{t,1}$,
although we note that the order of magnitude of this error is generically
quite large. Constraints on this parameter are bad (only a loose bound
$0<n_{t,\,1}<10$ at 95% CL, where the upper bound is set by the prior) mainly
because of the particular choice $n_{t,1}=4$. The growth of the spectrum being
quite steep, bounds on this parameter are driven by a limited number of bins
in frequency, making it harder to constrain large values of $n_{t,1}$ than
small values. Instead, the error on the high frequency tilt is smaller when
the amplitude has a peak at low frequency ($f_{*}<10^{-4}$ Hz). In fact, in
this case, the flat part of the spectrum would entirely lie within the LISA
band.

The impact of the foregrounds on the reconstruction seems less pronounced with
respect to the PL, except for a choice of $f_{*}$ close to the LISA
characteristic frequency. Such impact is even more reduced for $\delta$ and
$n_{t,2}$. The corner plot in fig. 7 shows the reconstructed mean value
parameters for the BPL-BNK_1. All mean values result within the 95% C.L.,
which reflects the accuracy of the reconstruction for the BPL case. In
particular, the high-frequency spectral tilt $n_{t,2}$ and the width $\delta$
can be measured with very high accuracy.
We notice a banana-shaped degeneracy among some parameters, in particular for
2D posteriors in the planes containing the smoothing parameter $\delta$. For
comparison, we also show the constraints obtained by fixing $\delta=1$, which
corresponds to the injected value (see purple contours in the figure). The
constraints are significantly better, although the example is for illustration
purposes only: we would need to develop theoretical arguments to fix such
(otherwise free) parameter. Nevertheless, the degeneracies do not impact the
reconstruction of the signal, of the noise, and of foregrounds. Importantly,
despite the non-Gaussian shape of the posterior distribution, the blue
contours show that our estimates for the errors, as obtained with the Fisher
analysis, are in good agreement with the MC analysis. Finally, notice that the
68 and 95 % CL error bands on the galactic foreground and instrumental noise
are small, while for the extragalactic foreground the errors are larger and
prior-dominated.

BPL-BNK_2. The benchmark BPL-BNK_2 represents a broad bump, asymmetric around
its peak. We find that the BPL recostrunction of this kind of signal
introduces several complexities that the Fisher approximation does not capture
unless the SNR is exceedingly large. For this reason, we opt not to show the
Fisher reconstruction maps but directly inspect the corner plot for our nested
sampling analysis run on BPL-BNK_2 (see fig. 8). We can notice that LISA can
accurately reconstruct the frequency shape of the injected signal. However,
the posterior distributions are extremely non-Gaussian in several parameter
directions. As expected, the Fisher formalism does not capture such complex
degeneracies.

In general, BPL-BNK_2 exhibits posteriors that are challenging to sample.
Their complexity likely require strategies that go beyond the scope of this
paper. For instance, broad bump profiles within the sensitivity window of LISA
may be over-parameterised by our BPL template, and might potentially be better
characterised using other templates.202020Although model selection
considerations go beyond the scope of this work, we checked that the BPL and
LN bump templates both fit the BPL-BNK_2 signal with similar $\chi^{2}$. An
alternative may be to expand the BPL template in eq. (13) around $f=f_{*}$,
and build combinations of $\left\\{\delta,\,n_{{t},1},\,n_{{t},2}\right\\}$
that control linear and quadratic departures from $f_{*}$. This procedure
would likely reduce the degeneracy between parameters, but the new
parameterisation would not necessarily be optimal for the search of benchmarks
like BPL-BNK_1.

###  4.4 Forecasts for the double peak template

Figure 9:  Fisher forecasts for the DP template. The panels show the 68 % CL
reconstruction error on $h^{2}\Omega_{*}$ (top left), $f_{*}$ (top right),
$\beta$ (next-to-top left), $\kappa_{1}$ (next-to-top right), $\rho$ (next-to-
bottom left), $\kappa_{2}$ (next-to-top right) and $\gamma$ (bottom) as a
function of the $\mathcal{P}_{\zeta}^{\rm ln}$ parameters $\Delta$ and
$\mathcal{A}_{s}$ (and $k_{*}$ adjusted to yield $f_{*}=5\,$mHz). SNR contour
lines are plotted in white. Depending on the panel, the pairs of dashed lines
mark the absolute error $\sigma=0.3$ ($\sigma=0.01$) and relative error 30%
($1\%$) contours, respectively in the absence [white] and in the presence
[black] of foregrounds. The white crosses display the benchmarks DP-BNK_1 and
DP-BNK_2.

Figure 10: Like in fig. 9 but with errors as a function of the
$\mathcal{P}_{\zeta}^{\rm ln}$ parameters $k_{*}$ (expressed in terms of
$f_{*}$) and $\mathcal{A}_{s}$ with fixed $\Delta=0.2$.

We now discuss the Fisher forecasts for a SGWB characterised by a DP template,
presented in section 2.4. The DP template is given by eq. 14 and is
characterised by seven parameters, whose geometrical meaning is explained in
section 2.4. Their impact on the signal shape is visualised in fig. 22 in
appendix A.

As previously stressed, SGWB profiles featuring a double peak structure can be
produced by curvature fluctuations at second order in perturbations if their
spectrum is sufficiently narrow. This implies that a given narrow curvature
power spectrum $\mathcal{P}_{\zeta}$ (with all its parameters fixed)
univocally defines the values of the seven DP parameters and, in turn, varying
the parameters of that given $\mathcal{P}_{\zeta}$ defines a hypersurface in
the DP parameter space. Thus, unlike the Fisher forecast maps in the previous
sections, for which we scan the space of the template parameters by varying
them directly, here we find it instructive to slice the DP template parameter
space by assuming the log-normal curvature power spectrum
$\mathcal{P}_{\zeta}^{\rm ln}$ of eq. (15), which only depends on three
parameters, _i.e.,_ the overall amplitude $\mathcal{A}_{s}$, the peak width
$\Delta$ and the peak position moved via $k_{*}$.212121Note that the pivot
scale $f_{*}$ in the DP template marks the transition between the two peaks
and it is not exactly equal to the peak of the primordial power spectrum.
However, the two are trivially related through the rescaling
$f_{*}=2/(\kappa_{2}\sqrt{3})f_{\rm PS}$ (where $f_{\rm PS}$ is the frequency
associated with the peak of the power spectrum and we have assumed standard
radiation domination after inflation). Section B.2 provides the numerical
mapping relating the three $\mathcal{P}_{\zeta}^{\rm ln}$ parameters to the
seven parameters of the DP template.

In fig. 9 we report the Fisher reconstruction uncertainties on the DP template
parameters as a function of the values of $\\{\mathcal{A}_{s},\Delta\\}$ when
$f_{*}=5\,\mathrm{mHz}$. To produce the figure, we fix the template parameter
$f_{*}$ and for each pair of values of $\\{\mathcal{A}_{s},\Delta\\}$ we
compute the six left DP template parameters. Subsequently, we input the values
of these seven parameters into our Fisher code that determines the 68 % CL
reconstruction errors on the DP parameters. With this choice of $f_{*}$, the
first peak happens to be close to the one of the galactic foreground (around
$1\,\mathrm{mHz}$) while the second and higher peak is always above the
galactic foreground and comparable to the extragalactic one only for
$\log_{10}h^{2}\Omega_{*}\lesssim-12$, approximately corresponding to
$\log_{10}\mathcal{A}_{s}\lesssim-3$ (cf. top panel in fig. 11). As a
consequence, including (or excluding) the foreground reconstruction in the
analysis largely influences the estimate of the parameters associated with the
first and second peak, as discussed in detail below.

The top-left panel of fig. 9 highlights that errors associated with the
overall amplitude only weakly depend on $\Delta$, and the 0.01 and 0.3 error
contour lines computed in the presence (dashed white) or the absence (dashed
black) of foregrounds slightly differ from each other for
$\log_{10}\mathcal{A}_{s}\lesssim-3.3$. This can be explained by the fact
that, in the vicinity of such a value of $\mathcal{A}_{s}$, the extragalactic
foreground starts covering the second and higher peak of the template. For the
same reason, the parameters associated with the second peak all display a
similar behaviour (see below).

The other panels show the errors on the parameters controlling the double peak
shape. For a given amplitude $\mathcal{A}_{s}$, they reach minimal values at
$\Delta\simeq 0.2$ – $0.4$. Departing from this range of $\Delta$ implies
larger errors. Indeed, the SNR decreases towards small $\Delta$, while the two
peaks tend to merge towards high $\Delta$, so that the overall structure
becomes more problematic to reconstruct (and several degeneracies among
parameters are expected). Only the error on $\gamma$ has a mild dependence on
$\Delta$. The parameter determines the behaviour of the very steep ultraviolet
tail which fast falls below the LISA sensitivity. Only a very loud signal with
$\log_{10}\mathcal{A}_{s}\gtrsim-2$ (SNR $\gtrsim 10^{5}$) would allow us to
reconstruct $\gamma$ below the percent level (for this case, $1\,\%$ error
lines end outside the panel for the plotted range).

Figure 9 moreover shows how foregrounds deteriorate the reconstruction by
including error lines from an analysis done in absence of them. As a rule of
thumb, for all parameters but $\gamma$, a reconstruction error below $30\%$ in
the range $\Delta\simeq 0.35$ requires SNR $\gtrsim\,50$ when foregrounds are
absent (_i.e.,_ perfectly known a priori). The inclusion of foregrounds
significantly affects the errors on the parameters associated with the first
peak, _i.e.,_ $f_{*},\beta,\kappa_{1}$, already for
$\log_{10}\mathcal{A}_{s}\lesssim-2.0$ (roughly corresponding to
$h^{2}\Omega_{*}$ in the range $[10^{-10},10^{-9}]$) and in particular for
values of $\Delta\geq 0.5$. As already mentioned, for $f_{*}=5\,\mathrm{mHz}$
and large $\Delta$, the first peak of the template competes with the peak of
the galactic foreground, resulting in a more challenging reconstruction for
this part of the primordial signal. In contrast, error lines on the parameters
associated with the second peak, _i.e.,_ $\rho,\kappa_{2},\gamma$, feel the
presence of the foregrounds in the analysis only when the overall amplitude is
decreased, _i.e.,_ $\log_{10}\mathcal{A}_{s}\lesssim-3.3$, and the second peak
amplitude becomes comparable with the one of the extragalactic foreground.

To complete our analysis on the Fisher forecast, in fig. 10 we show the
uncertainties on the DP template parameters as a function of the values
$\\{\mathcal{A}_{s},f_{*}\\}$; namely we fix the template shape by setting the
log-normal power spectrum parameter to $\Delta=0.2$ and we change the overall
amplitude of the signal and the scale setting its position in the LISA band.
Thus, in all panels, error lines follow the behaviour of the ones of constant
SNR.222222Note the similarity between the panels in fig. 10 and the first
three panels in fig. 4. As in the current case, there we scan overall
amplitude and pivot scale of a template with fixed shape. Furthermore, for all
parameters, the error has its minimum, for a fixed overall amplitude, around
$\log_{10}f_{*}\simeq-2.5$ with a small shift towards higher $f_{*}$ if the
parameter is related to the first peak (note that $\log_{10}(h^{2}\Omega_{*})$
is determined by the second peak), the reason being that thanks to this shift,
the first peak is centered where LISA has its maximum sensitivity. The
parameters which are easily reconstructed are the overall amplitude and the
pivot scale, while $\gamma$ remains the most elusive, although large errors on
$\gamma$ do not compromise the reconstruction of the overall structure of the
template. Errors on the parameters $\beta,\kappa_{1},\kappa_{2},\rho$ have
similar behaviours; for the sweetest spot in $f_{*}$ errors are below $30\%$
for $\log_{10}\mathcal{A}_{s}\gtrsim-3.3$. Finally, the effects of foregrounds
are relevant for $\log_{10}(f_{*}/\mathrm{Hz})\lesssim-2.4$ and for
$-3.4\lesssim\log_{10}\mathcal{A}_{s}\lesssim-2.5$. As before, for this range
of $f_{*}$ the extragalactic and (mainly) the galactic foreground can bury the
primordial signal at some LISA frequencies, depending on the amplitude of the
primordial signal.

Benchmarks.— Inflationary scenarios with a narrowly peaked primordial power
spectrum are among the main motivations for an SGWB following the DP template.
We focus on the narrow log-normal power spectrum $\mathcal{P}_{\zeta}^{\rm
ln}$ with $\Delta=0.2$ to set our DB benchmarks. The choice $\Delta=0.2$ fixes
five out of seven parameters of the DP template, namely
$\beta=0.242,\,\kappa_{1}=0.456,\,\kappa_{2}=1.234,\,\rho=0.08$ and
$\gamma=6.91$. We then consider the case $\mathcal{A}_{s}=10^{-2}$
($\mathcal{A}_{s}=10^{-3}$) as a regime where the curvature power spectrum
seeds a PBH abundance around (much below) the observed dark matter relic
density (see section 5.2). These two values of $\mathcal{A}_{s}$, combined
with $\Delta=0.2$, set the DP parameter $h^{2}\Omega_{*}$ to $-9.35$ and
$-11.35$. In addition, we consider the cases $f_{*}=1\,\mathrm{mHz}$ and
$f_{*}=5\,\mathrm{mHz}$ to illustrate how the galactic foreground, peaking
around $1\,\mathrm{mHz}$, impacts the parameter reconstruction. We denote as
DP-BNK_{1, 2, 3, 4} the DP benchmarks with
$\\{\log_{10}(h^{2}\,\Omega_{*}),\,f_{*}/\mathrm{mHz}\\}=\\{(-9.35,\,5),\,(-11.35,\,5),\,(-9.35,\,1),\,(-11.35,\,1)\\}$,
and $\beta,\kappa_{1},\kappa_{2},\rho$ and $\gamma$ fixed as above.

Figure 11:  1D and 2D posterior distributions derived from the DP template
reconstruction of the benchmarks DP-BNK_1 (top left), DP-BNK_2 (top right),
DP-BNK_3 (bottom left) and DP-BNK_4 (bottom right) defined as in eq. 14 with
$\\{\log_{10}h^{2}\Omega_{*},\,f_{*}/\textrm{Hz}\\}$ equal to $\\{-9.35,5\\}$,
$\\{-11.35,5\\}$, $\\{-9.35,1\\}$ and $\\{-11.35,1\\}$, respectively, and
$\\{\beta,\,\kappa_{1},\,\kappa_{2},\,\rho,\,\gamma\\}=\\{0.242,\,0.456,\,1.234,\,0.08,\,6.91\\}$
fixed for all benchmarks. Lines’ styles and color codes are as in fig. 3.

In fig. 11 we present the DB-template-based reconstruction for these
benchmarks. For the DP-BNK_1 and DP-BNK_3, which have large amplitude, all the
parameters are well reconstructed and consistent with the injected values. The
consistency is at $95\%\ \rm{C.L.}$ when $f_{*}=5\,\rm{mHz}$ for DP-BNK_1 and
$68\%\ \rm{C.L.}$ for DP-BNK_3. In both cases, some degeneracies are present,
but they are expected given the complex template parameter dependencies. In
particular, the parameter $\gamma$ exhibits several degeneracies with other
parameters (mainly $\rho$ and $h^{2}\Omega_{*}$). For instance, opportune
variations of the width and skewness of the second peak, $\gamma$ and $\rho$,
can mimic a change in the overall amplitude $\mathcal{A}_{s}$ (see fig. 22 in
appendix A), and this degeneracy is reflected in the corner plot.
Nevertheless, we find that all these degeneracies do not jeopardise the
reconstruction performance. Remarkably, all the 1D posteriors are Gaussian:
this fact, combined with a large SNR signal, reflects a good agreement between
the nested sampling and the Fisher analysis. Despite the large signal, both
foregrounds appear very well reconstructed, as well as the noise parameters
(_i.e.,_ unobservable error bars). This result is likely due to the DP
template’s peculiar shape which does not resemble the one of the foregrounds.

The situation changes when we consider the benchmarks DP-BNK_2 and DP-BNK_4,
which have lower amplitudes than the previous benchmarks. For DP-BNK_2, the
pivot frequency $f_{*}$ is sufficiently large for the signal not to be covered
by the galactic foreground despite its small amplitude. As a consequence, all
the DP parameters are well reconstructed, with the mean values consistent with
the injected values at 68 % CL levels. Also, in this case, the degeneracies do
not jeopardise the likelihood sampling. Moreover, the Fisher approximation is
still able to recover the shapes of the posterior distribution of the
parameters, slightly overestimating the errors.

The reconstruction drastically deteriorates for DP-BNK_4. Due to
$f_{*}=1\,\,\mathrm{mHz}$, the signal is entirely covered by the galactic
foreground. Several posterior distributions show a non-Gaussian behaviour and
some reconstructions are loose. In this regime, the Fisher formalism is unable
to capture the features of the posterior distribution of the model parameters.

###  4.5 Forecasts for the excited state template

Figure 12: Fisher forecasts for the ES template. The panels show the 68 % CL
reconstruction error on $h^{2}\Omega_{*}$ (left), $\omega_{\rm ES}$Hz (middle)
and $\gamma_{\rm ES}$ (right) as a function of the injected values of
$h^{2}\Omega_{*}$, $\omega_{\rm ES}$Hz and $\gamma_{\rm ES}$ set as specified
in the axes and the title of each panel. SNR contour lines are plotted in
white. The pairs of dashed (dotted) lines mark the $\sigma=0.3$
($\sigma=0.01$) contours, respectively in the absence (white) and in the
presence (black) of foregrounds. The crosses display the benchmarks ES-BNK_1,
ES-BNK_2 and ES-BNK_3.

Figure 13: 1D and 2D posterior distributions derived from the ES template
reconstruction of the benchmarks ES-BNK_1 (top left), ES-BNK_2 (top right) and
ES-BNK_3 (bottom). Lines’ styles and color codes are as in fig. 3. ES-BNK_1,
ES-BNK_2 and ES-BNK_3 are defined as in eq. 17 with
$\\{\log_{10}h^{2}\Omega_{*},\gamma_{\mathrm{ES}},\omega_{\mathrm{ES}}\mathrm{Hz}\\}$
set at $\\{-11,10,6\cdot 10^{3}\\}$, $\\{-11.5,10,6\cdot 10^{3}\\}$ and
$\\{-11,10,1.2\cdot 10^{3}\\}$, respectively.

The ES template described in eq. 17 is characterised by the three parameters
$h^{2}\Omega_{*},\gamma_{\mathrm{ES}}$ and $\omega_{\mathrm{ES}}$. Our Fisher
reconstruction forecast for this template is presented in fig. 12.
Reconstruction errors on $h^{2}\Omega_{*}$ and $\omega_{\mathrm{ES}}$ below 1
– 10% typically require SNR $\gtrsim$ 50 when
$3\leq\log_{10}{\omega_{\mathrm{ES}}}\leq 4$. We recall that the peak of the
template is related to the frequency of the oscillations in the ultraviolet
tail by $f_{\mathrm{max}}\simeq 6/\omega_{\mathrm{ES}}$ (see section 2.5).
This relationship explains why the relative errors on the various parameters
are minimized at $3\leq\log_{10}(\omega_{\rm ES}\,\mathrm{Hz})\leq 4$, which
range corresponds to signals peaking at
$f_{\mathrm{max}}\simeq(0.6-6)\,\mathrm{mHz}$ where the sensitivity of LISA is
optimal. Reconstructing the amplitude $h^{2}\Omega_{*}$ with the same
precision, on the other hand, requires larger amplitude values compared to
$\omega_{\mathrm{ES}}$. Finally, the parameter $\gamma_{\mathrm{ES}}$, which
controls the extension of the oscillations, is harder to measure: at $f_{\rm
max}\sim 1\,{\rm mHz}$, the error on $\log_{10}\gamma_{\mathrm{ES}}$ is below
0.3 only if $\log_{10}\,h^{2}\Omega_{*}\gtrsim-11.5$. Indeed, if the amplitude
is not sufficiently large, the ES oscillatory features are not detectable as
they fall outside the LISA sensitivity. As for the foregrounds, they affect
the errors on the parameters primarily in the region $\log_{10}(\omega_{\rm
ES}\mathrm{Hz})\gtrsim 3.5$, for which the signal peaks at a frequency
$\lesssim 2\,\mathrm{mHz}$, overlapping with the galactic foreground.

Benchmarks.— On theoretical grounds the parameter $\gamma_{\mathrm{ES}}$ is
expected to be much larger than one, but cannot be arbitrarily large. In view
of the microphysical setups motivating the template in section 2.5, we fix
$\gamma_{\mathrm{ES}}=10$ for all benchmarks. To test ES scenarios peaking at
frequencies lower and higher than the maximum of the galactic foreground, we
consider the choices $\omega_{\mathrm{ES}}\,=6\,\mathrm{mHz}^{-1}$ and
$\omega_{\mathrm{ES}}\,=6/5\,\mathrm{mHz}^{-1}$, respectively leading to an ES
main peak at $1\,\mathrm{mHz}$ and $5\,\mathrm{mHz}$. For the former choice,
in order to investigate the impact of the galactic foreground when the ES
signal is just or much below it, we consider two amplitude scenarios, namely
$\log_{10}(h^{2}\Omega_{*})=-10.5$ and $\log_{10}(h^{2}\Omega_{*})=-11.5$.
This leads us to consider benchmarks ES-BNK_1, ES-BNK_2 and ES-BNK_3 whose
parameters
$\\{\log_{10}(h^{2}\Omega_{*}),\gamma_{\mathrm{ES}},\,\omega_{\mathrm{ES}}\,\mathrm{mHz}\\}$
are respectively set at $\\{-11,\,10,\,6/5\\}$, $\\{-10.5,\,10,\,6\\}$ and
$\\{-11.5,\,10,\,6\\}$.

We first discuss the reconstruction of ES-BNK_1 (see top left panel in fig.
13). The signal peaks at $5\,{\rm mHz}$, where the galactic foreground has
small amplitude, and stands above the extragalactic one. The parameters
$h^{2}\Omega_{*}$ and $\omega_{\rm{ES}}$ are well constrained with the
reconstructed mean values within 95 % CL error bars. Their 2D posterior is