The calculation of the thermal effective potential, starting from an underlying particle physics model, is itself a nontrivial task. Due to the infrared (IR) sensitivity of bosonic fields at high temperature, perturbation theory must be re-summed. As a consequence, the effective expansion parameter increases $g^2 \to g$ (or even $g^2 \to \sqrt{g}$), and the expansion converges more slowly. Typical one-loop calculations utilise daisy re-summation (Parwani 1992; Arnold and Espinosa 1993) to include the leading $O(g^2)$ and $O(g^3)$ contributions to the potential, while methods to include still higher order contributions have been developed (Parwani 1992; Arnold and Espinosa 1993; Farakos et al. 1994; Braaten and Nieto 1995; Kajantie et al. 1996a, 2003; Croon et al. 2021; Ekstedt and Löfgren 2020; Gould 2021; Ekstedt et al. 2023). In addition to the slow convergence of the perturbative expansion, there are true IR divergences at $O(g^6)$ (four-loop order) which re-summation does not resolve (Linde 1980), underlying the importance of nonperturbative calculations (Gould et al. 2019; Kainulainen et al. 2019; Niemi et al. 2021; Halverson et al. 2021; Huang et al. 2021).

The bubble nucleation rate can be expanded around the percolation time $t_*$ as

$$ \Gamma \sim T^4 e^{-S_3/T} = \Gamma_* e^{-\beta(t-t_*)}, \quad (38) $$

and adiabaticity of the expansion of the universe yields $dT/dt = -TH$, so

$$ \frac{\beta}{H_*} = T_* \frac{d(S_3/T)}{dT} \bigg|_{T_*} . \quad (39) $$

This parameter determines whether the transition will complete mostly by the expansion of a few nucleated bubbles or by the nucleation of new bubbles everywhere in space. Indeed, the larger $\beta/H_*$ is, the faster the nucleation rate increases in time, the more bubbles will nucleate inside a Hubble horizon before the completion of the phase transition, and the smaller will their radii be when they collide, since there had not been much time for them to expand before meeting a neighbour. Thus one can expect an inverse relation between $\beta$ and the amplitude of the GW spectrum. Moreover, $\beta$ enters in the determination of the SGWB peak frequency. Note that, for very strong transitions, the definition above may become inappropriate and should be reformulated (Huber and Konstandin 2008b; Jinno et al. 2017).

As with $\alpha$, the calculation of the bubble nucleation rate, $\Gamma$, is affected by the IR sensitivity of bosonic fields at high temperatures. The spatial inhomogeneity of bubbles, as well as the time dependence of their creation, raises additional challenges. While much effort has gone into overcoming these challenges (Gleiser et al. 1993; Bodeker et al. 1994; Berges et al. 1997; Surig 1998; Strumia and Tetradis 1999; Moore and Rummukainen 2001; Garbrecht and Millington 2015; Gould and Hirvonen 2021), the theoretical uncertainties present in calculations of the bubble nucleation rate are less well understood than those for the effective potential.

The determination of the wall velocity from first principles is a much more daunting task. The passage of the wall perturbs the equilibrium of particles in the plasma, so that the temperature, velocity and chemical potential of each species are different from the background. The distribution function of each species can be written as $f_i = f_i^{eq} + \delta f_i$,
and the perturbations can be calculated by solving the corresponding Boltzmann equation (Moore and Prokopec 1995; Konstandin et al. 2014),

$$ (p^\mu \partial_\mu + K_\mu \partial_{p_\mu}) f + C[f] = 0, \qquad (40) $$

where $K_\mu$ is (related to) the force term felt by the particles due to the passage of the bubble, and $C[f]$ are collision terms that must be computed from the microphysics of the model under consideration, taking into account the mutual interactions of particles present in the plasma at the electroweak phase transition (EWPT); see e.g. Moore and Prokopec (1995), Kozaczuk (2015). The Higgs equation reads (Konstandin et al. 2014)

$$ \Box \phi + \frac{dV(\phi, T)}{d\phi} + \sum_i \frac{dm_i^2}{d\phi} \int \frac{d^3 p}{(2\pi)^3 2E_p} \delta f_i = 0, \qquad (41) $$

and the term with fluid perturbations $\delta f_i$ acts as a friction countering the bubble expansion. From this equation one obtains the wall velocity and the wall width, typically by using an ansatz for the Higgs profile

$$ \phi(z) = \frac{\phi_0}{2} \left[ 1 - \tanh \left( \frac{z}{L_w} \right) \right] \qquad (42) $$

and solving for the width $L_w$ and the velocity $v_w$. Although the wall velocity has been computed from first principles for a few models, such as in SUSY extensions (John and Schmidt 2001; Kozaczuk et al. 2015), Higgs plus singlet setups (Kozaczuk 2015; Lewicki et al. 2022) and two Higgs doublet models (Dorsch et al. 2017), the determination of $v_w$ as a function of the parameters of the underlying theory remains an open issue for most of the interesting models discussed in the literature. Crucially, the collision terms appearing in the Boltzmann equation (40) have been calculated only for a few models, such as the SM (Moore and Prokopec 1995) and singlet extensions (Kozaczuk 2015), and even then are known only at a leading-log approximation. For most models, additional collision terms involving interactions of new particles among themselves and with the SM sector have not yet been computed. Altogether, it is fair to say that our current tools allow only for an estimate of the wall velocity, and there is still a large room for improvement on this front.<fnref n="9" />

## 6.3 GW sources

A FOPT in the early universe occurs when at least one scalar field gets trapped at the symmetric metastable vacuum. As the universe cools down, the scalar field tunnels through or thermally fluctuates over the barrier into the lower-energy vacuum, and at a certain temperature bubbles in the broken phase start to nucleate in the sea of the symmetric phase. The latent heat released in the transition drives these bubbles to

<footnote n="9">⁹ In this respect, a promising alternative direction is the use of a string theory-inspired method (Maldacena 1998) known as gauge/string duality, or holography. The power of this tool is that it allows for a first-principle determination of out-of-equilibrium observables in a four-dimensional quantum field theory. It was used in (Bea et al. 2021b; Bigazzi et al. 2021; Bea et al. 2022) to compute the bubble wall velocity in a family of strongly-interacting, four-dimensional gauge theories.</footnote>
expand further, and the transition completes when these bubbles collide and merge with each other. During bubble expansion, the latent heat is converted into the energy of the bubble wall and the thermal and kinetic motion of the surrounding plasma. However, due to spherical symmetry of the bubbles, no GWs are generated during this expansion period. It is only when the bubbles collide, at the end of the FOPT, that GWs are produced from the kinetic energy of the scalar field and the fluid. Since there are numerous independent collision regions today within the typical resolution of GW detectors, the resulting spectrum is stochastic.

The amount of GWs generated during a cosmological FOPT is quantified by the fractional energy density of these waves in the universe, $\rho_{GW}$, compared to the critical energy density, $\rho_c \equiv 3H^2/(8\pi G)$, where $H$ is the Hubble constant. One can rephrase the GW density parameter in Eq. (29) as

$$ \Omega_{GW} \equiv \frac{1}{\rho_c} \frac{d\rho_{GW}}{d\ln k}. \qquad (43) $$

For a stochastic source, the spectrum is proportional to the two-point correlation of the metric perturbations $\langle \dot{h}(\mathbf{k}, t)\dot{h}(\mathbf{k}', t) \rangle \equiv (2\pi)^3 \delta^{(3)}(\mathbf{k} - \mathbf{k}') P_{\dot{h}}(k)$, as

$$ \Omega_{GW} = \frac{k^3}{24\pi^2 H^2} P_{\dot{h}}(k). \qquad (44) $$

Since sufficiently distant regions could not have exchanged information with one another during the FOPT, the correlation function for low $|\mathbf{k}|$ corresponds to white noise, i.e. $P_{\dot{h}}(\mathbf{k}, t) \sim$ constant, and the spectrum grows as $k^3$. This growth cannot continue forever, so this pattern must eventually be broken and the spectrum must decrease at some point, for the total energy in GWs must be finite. Therefore, the SGWB spectrum must contain a peak. A typical shape is illustrated in Fig. 6 created using data from PTPlot version 1.0.1.<fnref n="10" />

There are three main sources of GWs from bubble collisions: (i) the kinetic energy of the scalar field along the bubble walls, (ii) coherent motion of the plasma generated by the bubble expansion, also known as “sound waves”, and (iii) turbulent motion of this fluid, typically expected to occur at late times after the fluid develops vortical motion. The total contribution to GW energy from FOPTs in the early universe is therefore the sum

$$ \Omega_{GW} = \Omega_{GW}^{(\text{coll})} + \Omega_{GW}^{(\text{sw})} + \Omega_{GW}^{(\text{turb})}. \qquad (45) $$

How these three contributions arise in principle depends on the detail of the model. However, as discussed in the previous section, it is known that a few parameters determined from microphysics are enough to discuss cosmological consequences: transition strength $\alpha$, nucleation speed $\beta$, wall velocity $v_w$, and the transition temperature $T_*$.

<footnote n="10">¹⁰ Since the release of PTPlot version 1.0.0 there have been several modifications of PTPlot. This includes taking into account the recent erratum from Hindmarsh et al. (2017a) which led to the amplitude of the GW signal being overestimated by a factor of 10. For a more complete list of changes in version 1.0.1, see ptplot.org.</footnote>
<figure><img src="image_10.png" /><figcaption>**Fig. 6** An example of the GW power spectrum from a FOPT, along with the LISA power-law-sensitivity (PLS) curve of SNR = 10. The model parameters used in this example are $v_w = 0.9$, $\alpha = 0.1$, $\beta/H_* = 50$, $T_* = 200$ GeV, $g_* = 100$</figcaption></figure>

and the energy budget of the transition (Espinosa et al. 2010; Wang et al. 2021a; Giese et al. 2020).

The transition dynamics can vary dramatically for different $\alpha$. In typical transitions, in which $\alpha$ is moderate, the friction on the wall balances against the pressure from the released energy before bubbles collide, and as a result the walls approach a constant terminal velocity. In this case most of the energy released in the transition is transferred to the surrounding plasma, and only a negligible fraction remains in the scalar field. Macroscopically the energy transferred to the plasma can be classified into heat and bulk kinetic motion of the fluid. The bulk motion takes the form of compression waves, or sound waves. As long as the linearity of the fluid equation holds, i.e. the kinetic energy in play is sufficiently low, sound waves from different bubbles continue to propagate even long after the transition completes, and they overlap to create random velocity field with typical length scale of the sound shell $l_{\text{shell}}$. This phenomenon is found in numerical simulations (Hindmarsh et al. 2014, 2015, 2017a), and (Hindmarsh 2018; Hindmarsh and Hijazi 2019) provide a corresponding analytic model based on the picture of an ensemble of overlapping shells. Interestingly, the velocity field sources GWs at the constant scale $l_{\text{shell}}$ even well after the completion of the transition (provided that the transition is not too strong). Therefore, it works as an efficient source of GW production. The resulting GW spectrum $\Omega_{\text{GW}}^{(\text{sw})}$ contains two important characteristic scales of the bubble size and the fluid shell thickness. A hybrid scheme to calculate the GW spectrum from this phenomena has been proposed recently (Jinno et al. 2021).

The sound waves are expected to decay through the formation of non-linearities in the flow, after which the fluid enters a turbulent phase. The non-linearities develop on a time scale $\tau_{\text{sh}}$, which is inversely proportional to the kinetic energy in the fluid. It has been shown that, even for weak transitions, $\tau_{\text{sh}}$ is typically smaller than a Hubble time (Ellis et al. 2020b), with $\tau_{\text{sh}}$ becoming smaller as $\alpha$ increases. The decay of sound waves produced during a FOPT and the subsequent development of a turbulent flow has not been studied so far using direct numerical simulations and remains a major source of uncertainty when determining the GW spectrum. This holds in particular for strongly FOPTs, when the fluid kinetic energy is expected to be high, and the contribution of kinetic turbulence to the SGWB signal is expected to be relevant. Within the sound wave regime, typically the larger the transition strength the greater the production of GWs. FOPTs with smaller $\beta$ which produce fewer,
larger bubbles are also preferred. For a detailed discussion on the dependence of the SGWB spectrum from sound waves on $\alpha$ and $\beta$ (see Caprini et al. 2020).

For deflagrations with $\alpha \sim \mathcal{O}(0.1)$ the sound shell picture breaks down (Cutting et al. 2020). In such transitions, substantial vorticity can be generated during the initial collision of sound shells and bubbles. Furthermore, hot high pressure regions of the metastable phase can form during the transition. These hot regions suppress further bubble nucleation and decrease the wall velocity at the initial stage of the transition (Konstandin and No 2011; Mégevand and Ramírez 2018), while they form droplets at the later stage, extending the duration of the transition and suppressing the production of kinetic energy. Both the generation of vorticity and the formation of droplets can modify the GW spectrum, either by speeding up the decay of sound waves or by reducing the fluid kinetic energy and thus GW spectrum amplitude. Further study is required in this regime to map out the production of vorticity and study the evolution of flows with mixed modes, as well as to find the suppression of kinetic energy from droplets in a wide range of parameter space. It is however worth noting that the wall velocity still needs to be computed as a function of parameter space for many BSM models. It is therefore unclear if there are many realistic models that would give rise to the strong transitions with subsonic wall velocities which would be affected by these results.

Another interesting scenario is extremely strong transitions, in which $\alpha \gg 1$ is realized (Randall and Servant 2007; Espinosa et al. 2008; Konstandin and Servant 2011a; Hambye and Strumia 2013; Jaeckel et al. 2016; Jinno and Takimoto 2017b; Marzola et al. 2017; Iso et al. 2017; Chiang and Senaha 2017; von Harling and Servant 2018; Bruggisser et al. 2018a, b; Hambye et al. 2018; Baldes and Garcia-Cely 2019; Hashino et al. 2019; Prokopec et al. 2019; Brdar et al. 2019b; Marzo et al. 2019; Breitbach et al. 2019; Baratella et al. 2019; Fairbairn et al. 2019) but the bubble walls still reach a terminal velocity because of the higher-order friction terms on the wall (Bodeker and Moore 2017; Höhe et al. 2021). In this case the relevant hydrodynamical solution is strong detonation. While the dominant GW source comes from fluid kinetic motion, we must note several things: first, the onset of turbulence is earlier in this type of model (Ellis et al. 2020b); second, the overlap of sound shells, which is one of the requirements for the linear growth of the GW spectrum, may be somewhat delayed (Jinno et al. 2019b). Finally, in typical models there exists a maximum $\alpha$ for which transitions complete. Above this, vacuum energy dominates and the metastable state begins to inflate, prohibiting the percolation of the bubbles of the true vacuum (Ellis et al. 2019b).<fnref n="11" />

When $\alpha$ exceeds a certain threshold, or simply in vacuum transitions, the bubble walls continue to accelerate until they collide (“run away”). While recent studies of the next-leading friction terms on bubble walls (Bodeker and Moore 2017; Höhe et al. 2021) suggest that the parameter space for runaway is much smaller than previously thought (Bodeker and Moore 2009), such behaviour can still occur in extremely strong transitions (Ellis et al. 2019c, 2020c). In this case the only

<footnote n="11">¹¹ In cases where the phase transition through thermal bubble nucleation fails other mechanisms such as vacuum fluctuations (Lewicki et al. 2021) or spinodal instability (Bea et al. 2021a) can lead to completion of the transition which would have a very different phenomenology and clearly distinguishable GW spectra.</footnote>
contribution is from the bubble wall $\Omega_{\text{GW}}^{(\text{coll})}$. The expanding and colliding bubble walls are highly relativistic and much thinner compared to the typical bubble size. This observation motivates the so-called envelope approximation (Kosowsky et al. 1992; Kosowsky and Turner 1993), in which the shear stress is approximated to be localised in an infinitesimally thin shell at the bubble wall and disappears upon collision (Huber and Konstandin 2008a; Jinno and Takimoto 2017a). However, recent studies revealed that for runaway and vacuum transitions the collided region of intersecting bubbles has a rich structure and cannot be neglected (Weir 2016; Cutting et al. 2018, 2021; Jinno et al. 2019a). The GW spectrum that takes collided regions into account has been actively studied numerically for two bubble (Lewicki and Vaskonen 2020a, b) and many bubble collisions (Cutting et al. 2018, 2021) as well as (semi)-analytically (Jinno and Takimoto 2019; Konstandin 2018). One of the implications is that, while the peak of the spectrum does not grow in amplitude after collisions as occurs for sound waves, the spectrum may have a growing structure towards the IR, which enhances the detection prospect by LISA.

In addition to the aforementioned sources, magnetic fields could be present prior to or generated during a thermal FOPT. In the presence of a magnetic field, the turbulent motion of the primordial plasma would become dynamically coupled to the magnetic field, leading to MHD turbulence Brandenburg et al. (2017), Brandenburg et al. (1996), Christensson et al. (2001). Moreover, the hydrodynamic turbulent motion induced by the expansion of the bubbles arises in an ionised plasma, and this fact in itself can lead to MHD turbulence. The dynamical coupling between the velocity and the magnetic field plays an important role in determining the shape of the GW signal. The latter is affected by the presence of the magnetic field both through the aforementioned dynamical coupling, and because the magnetic anisotropic stresses produce GWs on their own, as studied in early analytic works (Deryagin et al. 1986; Kosowsky et al. 2002; Caprini and Durrer 2006; Gogoberidze et al. 2007; Caprini et al. 2009b). So far, the hydrodynamical simulations of thermal FOPTs have neglected magnetic fields. The early analytic works on the subject require one to make assumptions on the temporal correlation functions of the turbulence velocity and the magnetic field, which strongly affect the SGWB spectral shape (Caprini et al. 2009a). Previous analytical works were extended with updated modelling of the MHD turbulence in Niksa et al. (2018). Numerical simulations avoid resorting to analytical modelling of the MHD turbulence evolution, since they compute directly the solution to MHD equations (Roper Pol et al. 2020; Kahniashvili et al. 2021). Recent numerical simulations and previous analytic estimates agree on the shape of the GW spectrum at high frequencies, scaling as $f^{-8/3}$. However, the dynamical evolution of the magnetic field during the GW production can affect this slope (Roper Pol et al. 2020). At intermediate frequencies, the numerical simulations of Roper Pol et al. (2020) predict a spectrum $\Omega_{\text{GW}}(f) \propto f$, that eventually is expected to turn to $f^3$ at frequencies in the super-horizon range, due to causality (Caprini et al. 2009a). $\Omega_{\text{GW}}(f)$ evaluated from the numerical simulations grows as $f^3$ at early times, shifting to a linear increase $\propto f$ as the time progresses and GWs are building up. The detail of the transition from $f$ to $f^3$ is an active topic of research. Combining numerical simulations with our theoretical understanding of the MHD turbulence
dynamics will help to obtain an accurate prediction of the spectral shape of the SGWB produced by MHD turbulence in the future. As is the case for other sources, the SGWB by MHD turbulence also depends on the parameters of the FOPT. The strength of the FOPT and the size of the bubbles can be related to the kinetic energy density, and to the characteristic scale of turbulence, respectively.

This picture can get more complicated in the case of magnetic field production during the FOPT, and it would require an appropriate treatment of the magnetic field in both the false and true vacuum regions, which are dynamically converting through the bubble-driven FOPT. During the FOPT, inhomogeneities of the Higgs field give rise to the production of magnetic fields (Vachaspati 1991; Ahonen and Enqvist 1998; Stevens et al. 2008; Vachaspati 2021; Zhang et al. 2019; Yang and Bian 2022).<fnref n="12" /> The previous mechanism does not require CP violations, however these are expected in relation to baryogenesis through spontaneous lepton number symmetry breaking at a FOPT (Cohen et al. 1991). This parity violation leads to the production of helical magnetic fields (Vachaspati 1991; Cornwall 1997; Vachaspati 2021; García-Bellido et al. 2003). Such parity-violating turbulent sources lead to the production of circularly polarised GWs, studied analytically in Kahniashvili et al. (2005), Kahniashvili et al. (2008), Kisslinger and Kahniashvili (2015), Ellis et al. (2020a), and computed numerically in Kahniashvili et al. (2021). Further studies are required in this direction for a clear understanding of the GW and the polarisation signals produced from MHD turbulence as a function of the FOPT parameters, and a computation of turbulence in the plasma from first principles.

## 6.4 Discovery prospects beyond the standard model

GWs are a unique probe of the physics of the early universe. Unlike photons, gravitons were not in thermal equilibrium at early times, and may therefore directly encode information about events long before the time of last scattering. The characteristic broken power-law spectrum from a FOPT peaks at a frequency which can be related to a temperature and a time in the early universe via its redshift. Figure 7 shows this relation using the frequency peak of the acoustic spectrum, described in Caprini et al. (2020). It is seen that for appropriate fiducial values, a FOPT at the electroweak scale temperature, $T \sim 10^2$ GeV, sources a SGWB spectrum peaking at frequencies to which LISA has its best sensitivity (Grojean and Servant 2007; Caprini et al. 2016, 2019). The sensitivity to weak scale physics implies an opportunity to study FOPTs of several different kinds. In particular, several proposals for FOPTs that supply the out-of-equilibrium circumstance for baryogenesis are anchored to the weak scale. However, LISA is also sensitive to energy scales much beyond the electroweak scale: for exceptionally strong FOPTs the peak can fall outside the LISA sensitivity region but the tail of the signal still has a high SNR. In this case, FOPTs even at the MeV or PeV scale can be detected. This feature opens up potential complementary access in other experiments, as will be described in the next section.

<footnote n="12">¹² Current studies consider vacuum FOPTs. The extension to a thermal FOPT is not straightforward, and it requires further analysis.</footnote>
<figure><img src="image_11.png" /><figcaption>Fig. 7 LISA sensitivity to the SGWB frequency peak of FOPTs in the early universe. Here the fiducial values $v_{IW} = 1$ and $g = 106$ are used to relate the frequency peak of the acoustic spectrum (Caprini et al. 2020) to a time and temperature in the early universe, as a function of the transition rate parameter $\beta/H$.</figcaption></figure>

### 6.4.1 First-order phase transitions in the LISA window

FOPTs are studied in the context of a variety of models, which can be divided into two main categories: the spontaneous breaking of a gauge symmetry, and confinement due to the strong coupling of a gauge symmetry (and the related chiral symmetry breaking). The EWPT, which is particularly interesting in light of electroweak baryogenesis (see Sect. 6.4.2), is an example of the former. The transition temperature associated with electroweak symmetry breaking implies that the peak frequency of the GW spectrum falls within the LISA sensitivity window, as seen in Fig. 7. As is well known, the EWPT is a crossover transition in the SM. Therefore, new physics is required to induce a FOPT that can produce a SGWB. FOPTs can be studied in a model-independent effective-field-theory approach, adding non-renormalisable operators to the Higgs potential. These analyses show that for the EWPT to be sufficiently strongly first-order, the cut-off scale must be rather low, below the TeV scale (Zhang 1993; Grojean et al. 2005; Delaunay et al. 2008; Chala et al. 2018), therefore studies of the EWPT in SM extensions generically require the introduction of at least one new degree of freedom at the $O(100)$ GeV scale.

A weak-scale hidden sector featuring the spontaneous breaking of a gauge symmetry could also lead to a SGWB in the LISA window. An observable amplitude of the signal implies a significant fraction of the radiative degrees of freedom must have participated, hinting at a connection with DM (Schwaller 2015; Jaeckel et al. 2016; Croon et al. 2018; Bertone et al. 2020). The spectral form of such a signal is unlikely to give away much of its origin, though it can be related to particle spectra in specific models (Croon et al. 2018). Complementary search strategies may be applied if such a spectrum is detected—more about this in Sect. 6.5.

Confining gauge theories are thought to feature a FOPT if the confined degrees of freedom either only include the gauge fields themselves, or if confinement implies the breaking of a chiral symmetry $SU(N_f) \times SU(N_f)$ and the number of chiral fermions (dynamical at the transition temperature) exceeds $N_f \ge 3$. The foundation for this insight dates back to an argument based on the linear sigma model (Pisarski and Wilczek 1984) and has since been supported in part by lattice simulations
(Iwasaki et al. 1996). Studying the phenomenology of the transition is challenging, as effective theories either apply to the low-temperature (confined) phase and the high-temperature free phase. Moreover, an important role is played by the instantons of the strongly coupled theory. Initial explorations based on the linear sigma model indicate stronger SGWB are found for gapped spectra (Croon et al. 2019b). However, the predictions depend on the chiral symmetry breaking model used (Helmboldt et al. 2019).

Another alternative is solitosynthesis, a FOPT resulting from the growth of non-topological solitons stabilised by a conserved global charge (Griest and Kolb 1989). Unlike phase transitions resulting from the nucleation of critical bubbles, stable sub-critical populations of such solitons—Q-balls—exist, and accumulate charge until they reach a critical size and grow. This slow process may imply supercooled transitions still complete (Kusenko 1997; Croon et al. 2020). This scenario typically relies on the existence of a global symmetry under which the universe is asymmetrically charged. It has therefore been studied in the context of supersymmetry (Kusenko 1997; Postma 2002; Pearce 2012) and asymmetric DM (Huang and Li 2017; Croon et al. 2020).

As described in the previous section, the SGWB due to a cosmological FOPT can be described in terms of a small number of thermal parameters: characteristic temperature, bubble wall velocity, latent heat and a dynamical parameter such as a nucleation rate. These thermal parameters then predict a spectrum with a shape in most cases primarily defined by just two parameters: the peak frequency and peak amplitude. Therefore, the SGWB-inverse problem typically features degeneracies. However, resolution of the full spectrum and in particular the peak gives a unique picture of a crucial stage of the early universe, and may play an important role in answering some of the most fundamental questions. Further resolution of the underlying BSM physics can also be accomplished through complementary experimental observation.

## 6.4.2 Connection with baryogenesis

An (electroweak) FOPT may also be responsible for the generation of the observed baryon asymmetry. This provides an appealing connection between baryogenesis and the direct probe of the electroweak epoch in the early universe through GWs especially since the corresponding signals would peak in the LISA band. At such a first order EWPT, satisfying simultaneously the three Sakharov conditions for baryogenesis (baryon number violation, CP violation and a departure from thermal equilibrium) becomes possible: baryon number is not an exact symmetry of the SM, violated via non-perturbative processes involving gauge and Higgs fields.¹³ In the presence of CP violation, these processes are biased towards producing more baryons than anti-baryons, and because the FOPT introduces a departure from thermal equilibrium, the reverse mechanism is suppressed, ensuring that the generated asymmetry is not washed-out and remains to the present day.

Among these electroweak baryogenesis mechanisms, those most widely studied are non-local ones, involving charge transport (see Morrissey and Ramsey-Musolf 2012; Konstandin 2013 for reviews): CP violation in the scatterings of plasma
particles with the phase transition boundary yields an excess of fermion handedness, which diffuses into the electroweak symmetric phase and is there converted to a baryonic excess via SM baryon number violating processes known as sphalerons; the expansion of the FOPT bubbles then leads the net baryonic excess inside the bubble, where the asymmetry is frozen-out up to the present day. This last condition depends on the sphaleron rate inside the bubble, $\Gamma_{\text{sph}} \sim e^{-Esph/T}$ (with $E_{\text{sph}}$ the sphaleron energy). Sufficient suppression of $\Gamma_{\text{sph}}$ for successful baryogenesis (i.e. $\Gamma_{\text{sph}}$ much smaller than the Hubble expansion rate) yields the condition of a strong EWPT, $v/T \gtrsim 1$, where $v$ is the electroweak vacuum expectation value after the phase transition.

Up to very recently it has been assumed that these electroweak baryogenesis scenarios, as they rely on diffusion of particles in the plasma, were not effective for fast moving bubble walls (i.e. faster than the speed of sound of the plasma, with the bubbles consequently expanding as detonations). On the other hand, GWs are more favourably produced by stronger transitions resulting in faster bubble walls. This seemed to put a tension between the generation of a sizeable GW background from the FOPT (e.g. large enough to be observable by LISA) and successful baryogenesis. However, recent works (Cline and Kainulainen 2020; Laurent and Cline 2020; Dorsch et al. 2021) indicate that while diffusion is not as efficient for supersonic bubble walls, it still allows for successful electroweak baryogenesis. It has also been shown that the production of a large GW signal can be made compatible with electroweak baryogenesis for bubbles expanding as subsonic deflagrations (Dorsch et al. 2017), yet in this case the plasma friction against the moving wall must be large enough that even quite strong transitions remain subsonic, which requires the presence of sub-TeV particles with relatively large coupling to the Higgs—a tightly constrained scenario by now.

The strongest bound on non-local electroweak baryogenesis comes from the non-observation of an electric dipole moment (EDM) of the electron by highly precise experiments (Andreev et al. 2018). The additional sources of CP violation required by these transport mechanisms typically impact the electron EDM already at 2-loop order via Barr-Zee diagrams (Engel et al. 2013)—two orders above the SM CP violation contribution, which manifests itself only when all three families of quarks are involved, thus affecting EDMs at 4-loop only, being therefore highly suppressed and unconstrained by EDMs. Whether electroweak baryogenesis is still viable under such tight bounds remains to be investigated. There are, however, models in which it might be possible to avoid those constraints. One possible example involves coupling the CP violation to a vacuum expectation value of an additional scalar field which vanishes after the transition rendering the model safe from EDM constraints (Espinosa et al. 2012; Cline and Kainulainen 2013; Vaskonen 2017; Huang et al. 2018; Cline et al. 2021). Alternatively, CP violation can be secluded in a dark sector, and communicated to the visible sector in a way that the corresponding EDM contribution is suppressed or absent (Cline et al. 2017; Carena et al. 2019, 2020). Yet

<footnote n="13">¹³ More precisely, the baryonic current $J_B^\mu$ is anomalous, and the anomaly $\partial_\mu J_B^\mu = (N_c g^2)/(32\pi^2) F^{\mu\nu} \tilde{F}_{\mu\nu}$ allows for the creation of baryons through the dynamics of the non-abelian gauge fields.</footnote>
another way is to increase the temperature of the EWPT from symmetry non-restoration phenomena (Baldes and Servant 2018; Glioti et al. 2019; Matsedonskyi and Servant 2020; Biekötter et al. 2021; Carena et al. 2021). The observable effect is a GW peak position shifted towards higher frequencies.

Finally, there are also local baryogenesis mechanisms where baryon production and CP violation occur at the same point in space, involving non-trivial gauge and Higgs field dynamics (Tranberg 2011; Konstandin and Servant 2011b). However, these mechanisms are highly inefficient at electroweak scale temperatures (Lue et al. 1997), although they can still be relevant in more exotic scenarios of very supercooled FOPTs. The main advantage of such cold electroweak baryogenesis mechanisms is their ability to generate the observed baryon asymmetry with just the SM CP violation from the Cabibbo-Kobayashi-Maskawa matrix (Tranberg 2011), or from the strong CP phase (Servant 2014), thus avoiding the tight bounds from EDMs (Andreev et al. 2018). Since the thermal plasma is very diluted at the end of the transition due to the supercooling, the GWs are sourced mainly by the kinetic energy of the bubbles.

## 6.5 LISA complementarity to other experimental tests

The potential of LISA to explore FOPT is complementary to other experimental tests such as searches for beyond SM physics at colliders, tests from flavour physics including neutrinos physics and B mesons physics or searches for signals from the dark sector. LISA detection of SGWB from FOPTs can not only explore the electroweak energy scale but could explore higher energy scales which could be inaccessible to present and future colliders. The other GW observatories available at the time that LISA flies, can also strengthen and complement the LISA findings.

Figure 8 sketches the FOPT energy scales that LISA can probe with SNR >10 as well as its complementarity with other GW observatories expected in the early 2040s. The figure also displays the parameter region ruled out by current BBN, pulsar timing array (PTA) and LIGO/Virgo data.<fnref n="14" />

The EWPT in the SM is a second order phase transition. However, there is a wealth of BSM scenarios in which the EWPT is of the first order leading to sizeable SGWB signals. Without claiming to be exhaustive, one can for example mention the extension of the scalar sector of the SM with extra singlet(s) (Craig et al. 2013; Profumo et al. 2015; Craig et al. 2016; Hashino et al. 2017; Beniwal et al. 2017; Kang et al. 2018; Matsui 2018; Alves et al. 2018; Beniwal et al. 2019; Ahriche et al. 2019; Chala et al. 2019; Chiang and Lu 2020; Alves et al. 2020, 2021; Shajiee and Tofighi 2019; Morais et al. 2020; Morais and Pasechnik 2020) (although some could come without collider traces, Ashoorioon and Konstandin 2009), some having been discussed in the context of LIGO (Balazs et al. 2017), or doublet(s) (Fromme et al. 2006; Dorsch et al. 2013, 2014; Laine et al. 2017; Wang et al. 2020b; Zhou and Bian 2022) (see also Fujikura et al. 2018 for an example from twin Higgs models),

<footnote n="14">¹⁴ The figure assumes that the SGWB hints arising in the present PTA analyses, are not of primordial origin. The energy scales detectable by LISA result broader than in Fig. 7 since scenarios with the SGWB frequency outside the LISA band can still fulfil the SNR >10 requirement; see Figueroa et al. (2018), Megias et al. (2020) for details.</footnote>
**Fig. 8** The parameter reach of present and future GW observational network for phase transitions in a runaway-like regime with $v_w \simeq 1$ and $\alpha \gg 1$. In each shaded area the corresponding experiment (see labels) detects the FOPT signal with SNR $> 10$. The region labelled BBN, PTA and aLIGO O3 are ruled out. Image adapted from Figueroa et al. (2018), Megías et al. (2020)

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

composite Higgs models (Chala et al. 2016; Bruggisser et al. 2018b; De Curtis et al. 2019; Bian et al. 2019; Xie et al. 2020), so-called fermionic extensions of the SM involving for example $SU(2)_L$ doublet and/or two $SU(2)_L$ singlet vector-like leptons strongly coupled to the Higgs boson (Angelescu and Huang 2019) or the various supersymmetric extensions of the SM (Davies et al. 1996; Apreda et al. 2002; Huber and Konstandin 2008b; Huang et al. 2015; Kozaczuk et al. 2015; Huber et al. 2016; Demidov et al. 2018).

There is also a wealth of theoretical BSM setups in which a FOPT can occur at a higher energy scale than the EW energy scale. These setups for example encompass extensions of the scalar sector (Jinno et al. 2016; Alanne et al. 2020; Baldes and Servant 2018), classical conformal approaches (Brdar et al. 2019b) or nearly conformally invariant field theories in which the generation of neutrino masses is linked to spontaneous scale symmetry breaking (Agashe et al. 2020), gauge extensions of the SM such as models including $U(1)_{B-L}$ (Madge and Schwaller 2019; Ellis et al. 2019c; Hasegawa et al. 2019; Dev et al. 2019; Bian et al. 2020; Ellis et al. 2019a, 2020c), left-right symmetries (Brdar et al. 2019a; Fornal 2021) or even mbigger symmetries (Huang and Zhang 2019), grand unified theories, non-supersymmetric (Croon et al. 2019a; Huang et al. 2020; Buchmüller et al. 2020a; Corianò et al. 2020; Okada et al. 2021) or supersymmetric (Haba and Yamada 2020), extra-dimension models and in particular the one including warped spacetimes (Randall and Sundrum 1999; Goldberger and Wise 1999; Rattazzi and Zaffaroni 2001; Creminelli et al. 2002; Garriga and Pomarol 2003; Randall and Servant 2007; Nardini et al. 2007; Konstandin et al. 2010; Konstandin and Servant 2011a; Bunk et al. 2018; Megías et al. 2018a, b) (often coming together with a stabilisation mechanism such as the Goldberger-Wise mechanism (Goldberger and Wise 1999) and holographic duals of the Composite Higgs models mentioned above) or gauge-Higgs unification approaches (Adachi et al. 2020).

Many of these BSM models are being explored at current particle colliders such as the LHC. Searches performed at the LHC experiments such as ATLAS and CMS
using the data from Run1 and Run2 did not find any evidence for BSM signals and
have allowed one to put constraints on masses and couplings of their corresponding
predicted new particles.

These BSM models could possibly be further explored at future colliders such as
the high-luminosity phase of the LHC or the Future Circular Collider,<fnref n="15" /> either by
direct detection of new particles or due to their impact on precision measurements of
many observables, such as masses and couplings in the different SM sectors (scalar,
gauge or fermionic). Still, most models have parameter space corners which will
remain inaccessible to present and near-future colliders, but whose early universe
dynamics results in GWs detectable by LISA. Figure 9 illustrates this complemen-
tarity in a singlet extension, showing the values of α and β/H* for different parameter
points and their prospective detection by LISA and HL-LHC.

Another exciting possibility of complementary observations comes from astro-
physics and more specifically gamma ray observations. These are connected with
intergalactic magnetic fields, as the latter would lead to pair production γ + **B** →
**e**⁺**e**⁻ followed by scattering of the charged particles off the magnetic field:
**e**⁺**B** → **e**⁺ + γ, to produce a secondary γ. This process leads to a dilution of the
spectra observed from distant blazars, which recently enabled a lower bound on the
magnitude of intergalactic magnetic fields (Tavecchio et al. 2010; Ando and Kusenko
2010; Neronov and Vovk 2010; Essey et al. 2011; Chen et al. 2015; Ackermann et al.
2018). Due to the lack of astrophysical processes to account for magnetic fields in
intergalactic voids, explanations through a primordial origin have been put forward.
One possibility is inflationary magnetogenesis (Turner and Widrow 1988; Ratra
1992; Martin and Yokoyama 2008; Kobayashi 2014) although it has been pointed out
recently that models where magnetic fields are produced before the EWPT
generically lead to problems with baryon isocurvature constraints (Kamada et al.
2021). Another interesting possibility involves a FOPT (Vachaspati 1991; Sigl et al.
1997; Tevzadze et al. 2012). While there are many differences between the magnetic
fields produced by both scenarios, notably involving the uncertainty of their
evolution after the transition (Durrer and Neronov 2013), the magnetic field is
governed by the same basic parameters as the GW signal. In particular this means
that the transition in question has to be strong enough to produce a magnetic field
amplitude above the known lower bounds. This, of course, leads to a certain
correlation between the two possible signals, and it was found that the parameter
space predicting a GW signal observable by LISA can also facilitate generation of a
magnetic field strong enough to explain the observed diffusion of the blazar spectra
for some models predicting a first order EWPT (Ellis et al. 2019a). Moreover, the
presence of helicity in the magnetic fields produced at a FOPT plays a crucial role on
its dynamical evolution, playing a crucial role in the resulting restrictions on the
strength of the FOPT (Brandenburg et al. 2017, 2019), as well as providing circularly
polarised SGWB. Such polarisation can also be measured by LISA providing a
distinct signal of this scenario (Ellis et al. 2020a; Domcke et al. 2020b).¹⁶

<footnote n="15">¹⁵ Although many of them can be elusive to collider physics (Addazi et al. 2019).</footnote>
<figure><img src="image_13.png" /><figcaption>**Fig. 9** Typical values of $\alpha$ and $\beta/H_+$ for the general singlet model, superimposed with LISA SNR curves for $T_*=50$ GeV (solid lines). The models in blue (orange) are unlikely (likely) to be probed by the high-luminosity LHC. The dotted straight lines are the contours of the fluid turnover time quantifying the effect of turbulence. In the grey shaded region the decay of sound waves into turbulence is less important than the Hubble damping and the SNR curve reflects this effect. Image reproduced with permission from Caprini et al. (2020), copyright by IOP/SISSA</figcaption></figure>

# 7 Cosmic strings

Section coordinators: M. Lewicki, L. Sousa. Contributors: P. Auclair, J.J. Blanco-Pillado, C-F. Chang, Y. Cui, D.G. Figueroa, A. Jenkins, S. Kuroyanagi, M. Lewicki, M. Sakellariadou, K. Schmitz, L. Sousa, D. Steer, J. Wachter.

## 7.1 Introduction

Cosmic defects are stable energy configurations that can be formed in early universe phase transitions, usually due to a spontaneous symmetry-breaking process driven by some scalar field(s) acquiring a non-zero expectation value within a (topologically non-trivial) vacuum manifold. If the symmetry broken is global, all non-constant field configurations produce energy-momentum, and are loosely referred to as *global defects*. When the broken symmetry is gauged, we speak of *local defects* instead. Here we will focus on cosmic strings, which correspond to stable one-dimensional defect solutions of field theories (Nielsen and Olesen 1973), independently of whether they are global or gauge. Alternatively, cosmic strings can also be identified with cosmologically stretched fundamental strings from String Theory, formed e.g. at the end of brane inflation (Dvali and Vilenkin 2004; Copeland et al. 2004).

The energy per unit length of a string is $\mu \sim \eta^2$, with $\eta$ a characteristic energy scale. The string tension, at least for the simplest cases, is typically of order $\mu$. In the case of topological strings, $\eta$ represents the energy scale of the phase transition. A network of strings formed in the early universe emits GWs throughout most of the history of the universe, generating a SGWB from the superposition of many

<footnote n="16">¹⁶ Other proposed measures for detecting magnetic helicity in voids involve exploring the twisting of photon paths coming from multiple blazar sources (Vachaspati 2021; Tashiro and Vachaspati 2015; Kahniashvili and Vachaspati 2006).</footnote>
uncorrelated emissions. Here we forecast the constraints that LISA may put on the dimensionless combination

$$G\mu \sim 10^{-6} \left( \frac{\eta}{10^{16} \text{ GeV}} \right)^2, \qquad (46)$$

where $G = 1/M_p^2$ is Newton's constant, and $M_p = 1.22 \times 10^{19}$ GeV the Planck mass. We note that various potential observational signatures of cosmic string networks, other than GW emission, have been discussed in the literature. These include anisotropies in the CMB (Ade et al. 2014; Charnock et al. 2016; Lizarraga et al. 2016; Ringeval 2010), lensing events (Vilenkin 1984; Bloomfield and Chernoff 2014), and cosmic ray emission (Brandenberger 1987; Srednicki and Theisen 1987; Bhattacharjee et al. 1992; Damour and Vilenkin 1997; Wichoski et al. 2002; Peloso and Sorbo 2003; Sabancilar 2010; Vachaspati 2010; Long et al. 2014). (see Hindmarsh and Kibble 1995; Vilenkin and Shellard 2000; Sakellariadou 2007; Vachaspati et al. 2015 for a review.) Currently, CMB data from the Planck Satellite (Ade et al. 2014) imply $G\mu < 10^{-7}$ for NG, Abelian-Higgs, and semi-local strings. The most stringent bounds, however, come from searches for the SGWB, with PTA constraining $G\mu$ for NG strings to be $G\mu \lesssim 10^{-11}$ (Blanco-Pillado et al. 2018a; Ringeval and Suyama 2017), and LIGO-Virgo observations constraining it to be as low as $G\mu < 2 \times 10^{-14}$, depending on the string network model (Abbott et al. 2019c, 2018a). As we will show, LISA will be sensitive to string tensions with $G\mu \gtrsim 10^{-17}$ for NG strings, improving current upper bounds by $\sim 10$ orders of magnitude relative to CMB constraints, by $\sim 6$ orders of magnitude relative to current PTA constraints, and even by $\sim 3$ orders of magnitude relative to future constraints from next generation of PTA experiments. One important caveat is that these predictions do not include the impact the astrophysical background of GWs. Boileau et al. (2022) shows that, depending on different assumptions on the galactic and extragalactic astrophysical foregrounds, the forecast sensitivity of LISA should be reduced to roughly $G\mu \gtrsim O(10^{-16} - 10^{-15})$. There is ongoing work in the consortium aimed at providing a refined pipeline for cosmic strings and sharpening these predictions.

As the characteristic width $\delta \sim 1/\eta$ of a cosmic string is much smaller than the size of the horizon, in the following we mainly assume that strings are well described by the NG action, which is the leading-order approximation when the curvature scale of the strings is much larger than their thickness. We refer to such strings as NG strings (we focus on string networks without junctions). Cosmic string networks are expected to reach an attractor solution known as scaling, for which the energy density of the network remains a fixed fraction of the background energy density in the universe. When strings within the network collide, they intercommute (i.e. “exchange partners”) and reconnect after the collision (technically speaking this corresponds to an intercommutation probability $P = 1$, which we mainly assume, though we will also comment briefly on the case $P < 1$). As a result, loops are formed whenever a string self-intersects or two curved strings collide. Sub-horizon sized loops decouple from the cosmological evolution and oscillate under their own tension, emitting GWs in this process. The relativistic nature of the strings typically
leads them to form cusps, corresponding to points on the string momentarily moving at the speed of light (Turok 1984). Furthermore, intersections of strings will also generate discontinuities in their tangent vector known as kinks. Cusps and kinks generate GW burst emissions (Damour and Vilenkin 2000, 2001), which add up all throughout cosmic history to form a SGWB.

Besides sub-horizon loops, the network also contains long strings that that stretch across a Hubble volume. These are either infinite or super-horizon loops, and they are also expected to emit GWs. However, the dominant contribution is generically produced by the superposition of radiation from many sub-horizon loops along each line of sight.

Cosmic string networks are expected to create one of the largest SGWB of cosmological origin known. In the following we analyse the ability of LISA to probe the SGWB emitted by a network of cosmic strings, considering leading models of the string networks as described in the literature. The prospects for detection or, alternatively, for setting up new stringent constraints on string parameters, are very encouraging.

## 7.2 Network modelling

The SGWB generated by cosmic string networks in a given frequency *f* has, at the present time, contributions from all the loops created throughout the history of the universe that have emitted GWs that have a frequency *f* today. The number density *n*(*l*, *t*) of non-self-intersecting, sub-horizon, cosmic string loops of invariant length *l* at any cosmic time *t* is then the crucial ingredient in this computation. In this section, we review the models for the loop number density that have been proposed in the literature.

### 7.2.1 Model 1

We start by considering an analytical approach—originally developed in Kibble (1985) and later extended in Caldwell and Allen (1992), DePies and Hogan (2007), Sanidas et al. (2012), Sousa and Avelino (2013)—which is based on two assumptions. The first assumption is that the production of loops is the dominant energy-loss mechanism in the evolution of the long string network. In this case, as we shall see, the loop production function *f*(*l*, *t*) *dl*—which gives us the number density of loops of lengths between *l* and *l* + *dl* produced per unit time and per unit volume—is determined by the large-scale evolution of the long string network. Here, we shall use the approach introduced in Sousa and Avelino (2013) and use the Velocity-dependent One-Scale model (Martins and Shellard 1996, 2002) to describe its evolution. This model provides an analytical description of the evolution of the characteristic length *L* ≡ (μ/ρ<sub>∞</sub>)<sup>1/2</sup>—where ρ<sub>∞</sub> is the energy density of the long string network—and of the root-mean-squared velocity v̄ of the network (Martins and Shellard 1996, 2002):
$$ \frac{d\bar{v}}{dt} = (1 - \bar{v}^2) \left[ \frac{k(\bar{v})}{L} - 2H\bar{v} \right], \qquad (47) $$

$$ \frac{dL}{dt} = (1 + \bar{v}^2)HL + \frac{\tilde{c}}{2}\bar{v}, \qquad (48) $$

with

$$ k(\bar{v}) = \frac{2\sqrt{2}}{\pi} (1 - \bar{v}^2) (1 + 2\sqrt{2}\bar{v}^3) \frac{1 - 8\bar{v}^6}{1 + 8\bar{v}^6}. \qquad (49) $$

The curvature parameter $k(\bar{v})$ accounts (to some extent) for the effects of small-scale structure (Martins and Shellard 2002).

Here, $\tilde{c}$ is a phenomenological parameter that quantifies the efficiency of the loop-chopping mechanism, that may be calibrated with simulations (for NG strings, $\tilde{c} = 0.23 \pm 0.04$ fits both radiation and matter era simulations (Martins and Shellard 2002)). Using Eq. (48), we find that the energy lost as a result of loop production is given by

$$ \left. \frac{d\rho_{\infty}}{dt} \right|_{\text{loops}} = \tilde{c} \bar{v} \frac{\rho_{\infty}}{L} = \mu \int_{0}^{\infty} lf(l, t) dl, \qquad (50) $$

and determines the normalisation of loop production function.

The second assumption of this model (which will be somewhat relaxed later) is that all the loops are created with a length $l$ that is a fixed fraction of the characteristic length of the long string network: $l = \alpha_L L$, where $\alpha_L < 1$ is a (constant) free-parameter of the model. We then have

$$ f(l, t) = \left(\frac{\mathcal{F}}{f_r}\right) \frac{\tilde{c} \bar{v}}{l L^3} \delta(l - \alpha_L L), \qquad (51) $$

where $\bar{v}$ and $L$ are given by Eqs. (47) and (48). Here the normalisation of $f(l, t)$ is determined using Eq. (51), except for the correction factor $f_r$ and $\mathcal{F}$. The first factor, $f_r \sim \sqrt{2}$, is introduced to account for the energy lost as a result of the redshifting of the peculiar velocities of loops (Vilenkin and Shellard 2000). Also, since in general we do not expect all loops to be created with exactly the same size but to follow a distribution of lengths, a second factor $\mathcal{F}$ is included to account for the effect of the spread of the distribution (see Sanidas et al. 2012; Sousa et al. 2020 for a detailed analysis of this effect.) For NG strings, $\mathcal{F}$ is estimated to be $O(0.1)$ (Blanco-Pillado et al. 2014).

Note that Eq. (51) is valid throughout any cosmological era, even during the radiation-to-matter and matter-to-dark-energy transitions (in which the network is not in a linear scaling regime). As a result, it allows us to compute $f(l, t)$ through cosmic history in a realistic cosmological background (Sousa and Avelino 2013). The number density of loops $n(l, t)$ at all times may then be found by solving
$$n(l, t) = \int_{t_i}^{t} dt' f(l', t') \left(\frac{a(t')}{a(t)}\right)^3, \qquad (52)$$

while accounting for the decrease in the length of loops caused by the emission of GWs $l = l' + \Gamma G\mu (t' - t)$ (where $\Gamma$ is a dimensionless constant characterising GW emission efficiency—that we define more precisely in Sect. 7.3—, $l$ is the length of a loop created at a time $t'$ with a length $l'$ at a later time $t > t'$).

The loop-size parameter $\alpha_L$ may either be calibrated with numerical simulations (as we shall do in Sect. 7.4), but may be also be treated as a free parameter of the model to explore a wider variety of scenarios (as in Sect. 7.5.4). Although $\alpha_L$ is the natural parameter of this model, we will express our results in terms of $\alpha = \alpha_L \xi_r$ (where $\xi_r$ is the value of $L/t$ in the radiation era) in order to ease comparison with other loop distribution models.

### 7.2.2 Model 2

The second model we consider is a loop number density distribution, $n(l, t)$, extracted from large scale numerical simulations of the string networks (Blanco-Pillado et al. 2014). These number densities were first obtained from the integration of the loop production function $f(l, t)$ of non-self-intersecting loops computed directly from the simulation. Recently, these loop number densities were also determined directly from simulations without the intermediate use of a loop production function (Blanco-Pillado and Olum 2020). Both approaches produce the same result.

In the following we give the resulting distributions for three types of loops: those created in the radiation era, emitting GWs in the radiation era; those created in the radiation era, emitting GWs in the matter era; and those created in the matter era, emitting GWs in the matter era.

With $\Omega_{rad}$ the fraction of the critical density in radiation, and understanding the redshift $z$ to be a function of time, we can write

$$n_{r,r}(l, t) = \frac{0.18}{t^{3/2}(1 + \Gamma G\mu t)^{5/2}} \Theta(0.1 - l/t), \qquad (53a)$$

$$n_{r,m}(l, t) = \frac{0.18(2\sqrt{\Omega_{\text{rad}}})^{3/2}(1+z)^3}{(l + \Gamma G\mu t)^{5/2}} \Theta(0.09 t_{\text{eq}}/t - \Gamma G\mu - l/t), \qquad (53b)$$

$$n_{m,m}(l, t) = \frac{0.27 - 0.45(l/t)^{0.31}}{t^2(l + \Gamma G\mu t)^2} \Theta(0.18 - l/t), \qquad (53c)$$

where $\Theta(x)$ is the Heaviside step function, and the subscript “r,m” indicates “loops produced in the (r)adiation era, emitting in the (m)atter era”, and so on. For Eqs. (53a) and (53c), the cutoffs are due to the maximum size of loop which can be produced in these eras. For Eq. (53b), there is an additional term in the cutoff which models the decay due to GW emission; that is, a loop formed in the radiation era must have some minimum size in order to survive into the matter era. Equation (53)
contains all the information necessary to compute the SGWB due to decaying loops, as outlined in Sect. 7.3.

Analysis of Eq. (53c), in comparison to Eqs. (53a) and (53b), shows it to be negligible for $G\mu \lesssim 10^{-10}$ and thus for cosmic string tensions of interest to LISA. The cosmic string generated SGWB potentially to be observed by LISA will be for radiation-era loops, integrated up to the present. Of particular interest are radiation-era loops emitting in the matter era; these loops contribute to a “bump” in the SGWB spectrum which may, depending on $G\mu$, fall within the LISA band (see Sect. 7.4). For this particular reason, but also more broadly, it is important to understand how gravitational backreaction affects the evolution of loops. At present, the decay of loops into GWs is accounted for by $\Gamma G\mu$ terms, but this does not take into account potential changes to the loop population due to, e.g. loop fragmentation as a result of backreacted trajectories. While current theoretical studies (Blanco-Pillado et al. 2018b, 2019b), and numerical studies of simple models (Robbins and Olum 2019), suggest that the change to loop trajectories will be small, it is important to keep in mind that no current, simulation-inferred network model incorporates gravitational backreaction directly.

### 7.2.3 Model 3

The third model we consider is based on the analytical studies of the small-scale structure and the correlation functions of individual infinite strings (Dubath et al. 2008; Polchinski and Rocha 2007, 2006). This analysis predicts a power-law loop production function in terms of scaling units, whose slope is parametrised by a parameter $\chi$ linked to the fractal dimension of the strings. The production of cosmic string loops is suppressed by gravitational back-reaction, on scales below

$$ \gamma_c = \Upsilon(G\mu)^{1+2\chi}, \qquad (54) $$

in which $\Upsilon = O(20)$. The characteristics of the loop production function were later inferred from measurements of the loop number density in large-scale simulations of NG strings (Ringeval et al. 2007; Lorenz et al. 2010).<fnref n="17" /> The numerically inferred loop-production function is then extended from large scales down to the gravitational back-reaction scale, resulting in an extra population of small loops.

The loop number density of Model 3 can be approximated by the following expressions

<footnote n="17">¹⁷ Note that the relationship between the loop number density and the loop production function is not always one-to-one, see Auclair et al. (2019) for a review.</footnote>
$$t^4 n_{r,r}(l,t) = \begin{cases} \frac{0.08}{(\gamma + \Gamma G \mu)^{3-2\chi_r}} & \text{if } \Gamma G \mu < \gamma \\ \frac{0.08(1/2 - 2\chi_r)}{(2 - 2\chi_r)\Gamma G \mu \gamma^{2-2\chi_r}} & \text{if } \gamma_c < \gamma < \Gamma G \mu \\ \frac{0.08(1/2 - 2\chi_r)}{(2 - 2\chi_r)\Gamma G \mu \gamma_c^{2-2\chi_r}} & \text{if } \gamma < \gamma_c \end{cases} \quad (55)$$

$$t^4 n_{m,m}(l,t) = \begin{cases} \frac{0.015}{(\gamma + \Gamma G \mu)^{3-2\chi_m}} & \text{if } \Gamma G \mu < \gamma \\ \frac{0.015(1 - 2\chi_m)}{(2 - 2\chi_m)\Gamma G \mu \gamma^{2-2\chi_m}} & \text{if } \gamma_c < \gamma < \Gamma_d G \mu \\ \frac{0.015(1 - 2\chi_m)}{(2 - 2\chi_m)\Gamma G \mu \gamma_c^{2-2\chi_m}} & \text{if } \gamma < \gamma_c \end{cases} \quad (56)$$

$$t^4 n_{r,m}(l,t) = \left(\frac{t}{t_{eq}}\right)^4 \left(\frac{1+z}{1+z_{eq}}\right)^3 t_{eq}^4 n_{\text{rad}}^{(3)} \left[\frac{\gamma t + \Gamma_d G \mu (t - t_{eq})}{t_{eq}}\right] \quad (57)$$

in which $(\chi_r, \chi_m) = (0.2, 0.295)$. This extra population of small loops dominates the SBGW in the high frequency spectrum as discussed in Abbott et al. (2018a), Auclair et al. (2020), Abbott et al. (2021b) and hence can lead to very different constraints on $G\mu$ to that of the two above-mentioned models. The energy density of these small loops is very large and the question of energy balance, in the context of the one-scale model, has been raised in Blanco-Pillado et al. (2019a).

### 7.2.4 Field theory

The above modellings are all based on NG strings, which are infinitely thin. Cosmic strings appear however naturally as solitonic solutions of classical field theory models (Nielsen and Olesen 1973), so in principle they can decay not only by GW radiation but also directly into excitations of their elementary constituents. For example, it has been observed with numerical simulations that global (axionic) strings decay into the massless Goldstone modes present in the vacuum of the theory (Davis and Shellard 1989). For local strings with no long-range interactions, the excitations in the vacuum are however massive, and hence are naturally expected to be suppressed for long wavelengths comparable to the length of the strings (Martins et al. 2004; Olum and Blanco-Pillado 2000, 1999). Furthermore, recent simulations of individual loops in the Abelian-Higgs model (Matsunami et al. 2019) report that extrapolating their results to large loop sizes, the GW emission should be expected to dominate over particle emission for loops larger than a certain critical length scale.

In contrast, large-scale field theory simulations of Abelian-Higgs strings (Vincent et al. 1998; Hindmarsh et al. 2009; Daverio et al. 2016; Hindmarsh et al. 2017b) observe that the network of long strings reaches a scaling regime, thanks to energy loss into classical radiation of the scalar and gauge fields involved. As a consequence of this, sub-horizon loops formed during the network evolution decay very promptly.
This intriguing discrepancy has been under debate for over twenty years, but the
origin of this radiation is not currently understood.

The similarities and differences between field theory and NG simulations of string
networks can then be summarised as follows: the infinite strings are rather similar in
curvature radius and length density, but loops decay into field modes in the field
theory simulations. In field theory simulations, the string energy density goes into
radiated modes of the fields, which do not belong to the string network anymore. As
a consequence, the string loops decay within a Hubble time, and hence do not
continue to contribute as a source of GWs. In the NG picture, this channel does not
exist, and instead the energy of the infinite strings goes into loops, which then decay
via GWs. Our analysis in the following is mostly based on the NG classical evolution
of strings. We assume, as supported by NG simulations, that loops are formed
throughout cosmic history, and they decay into GWs. Our discussion about the
ability of LISA to measure a GW background from cosmic strings is therefore based
on this fundamental assumption.

7.3 Computation of the gravitational-wave spectrum from loops

The incoherent superposition of GWs emitted from oscillating cosmic string loops
leads to a SGWB. The calculation of the SGWB from the cosmic string networks
have been widely studied in the literature (Vachaspati and Vilenkin 1985; Blanco-
Pillado et al. 2014; Blanco-Pillado and Olum 2017; Blanco-Pillado et al. 2018a;
Ringeval and Suyama 2017; Vilenkin 1981; Hogan and Rees 1984; Accetta and
Krauss 1989; Bennett and Bouchet 1991; Caldwell and Allen 1992; Siemens et al.
2007; DePies and Hogan 2007; Olmez et al. 2010; Sanidas et al. 2012, 2013;
Binetruy et al. 2012; Kuroyanagi et al. 2012, 2013; Sousa and Avelino 2013, 2014;
Cui et al. 2019; Chang and Cui 2020; Jenkins and Sakellariadou 2018; Gouttenoire
et al. 2020a). The generic form of a GW spectrum in logarithmic intervals of
frequency, introduced in Eq. (29), also applies to the case with cosmic strings. For the
cosmic string application, Eq. (29) is implemented by integrating over the GW
emission from all the loops throughout cosmic history that contribute to a certain
frequency.<fnref n="18" /> In the following subsections, we review two methods of calculating
SGWB from string loops that have been developed in the literature along with recent
updates. There are ingredients that are common to both approaches. The first
ingredient is $n(l, t)$, the number density of non-self-intersecting, sub-horizon, cosmic
string loops of invariant length $l$ at cosmic time $t$, As shown in Sect. 7.2, $n(l, t)$ can
be estimated by analytical or simulation-based methods. Another characteristic
function is the gravitational loop power spectrum $P_{\text{GW}}(f, l)$, which may be
determined as an average based on simulations, or approximated by focusing on the
analytical high-frequency behaviour of particular events on the strings, e.g. cusps and
kinks.

<footnote n="18">¹⁸ We neglect the GW contribution from long (horizon-spanning) strings as being subdominant to the GW contribution from loops for NG strings. See Buchmüller et al. (2013), Matsui and Kuroyanagi (2019).</footnote>
### 7.3.1 Method I

The first method to calculate $\Omega_{GW}(t_0, f)$ (Vachaspati and Vilenkin 1985; Blanco-Pillado et al. 2014; Blanco-Pillado and Olum 2017; Blanco-Pillado et al. 2018a; Caldwell and Allen 1992; Siemens et al. 2007; DePies and Hogan 2007; Olmez et al. 2010; Kuroyanagi et al. 2012, 2013; Sousa and Avelino 2013; Cui et al. 2019) assumes that $P_{GW}(f, l)$ takes the form of

$$P_{GW}(f, l) = G\mu^2 lP(fl), \qquad (58)$$

where $P(f, l)$ is an averaged function that can be computed from an ensemble of loops of length $l$ (but of different shapes) based on simulations.

Integrating all emissions throughout the history of the universe and taking into account the redshift effects, the GW energy density for a particular frequency $f$ as measured today is

$$\frac{dP_{GW}}{df}(t_0, f) = G\mu^2 \int_{0}^{t_0} dt \left(\frac{a(t)}{a_0}\right)^4 \int_{0}^{\infty} dl \, l n(l, t) P\left(\frac{a_0}{a(t)} fl\right). \qquad (59)$$

The GW radiation power of an isolated loop of length $l$ can be computed using the standard formulae in the weak gravity regime (Vachaspati and Vilenkin 1985; Vilenkin and Vachaspati 1987; Garfinkle and Vachaspati 1987). As a simple approximation, we assume that the loops evolve in flat space in a periodic manner, and thus emit GW at discrete frequencies

$$\omega_n = 2\pi n/T, \qquad (60)$$

where $T = l/2$ is the oscillation period, and the harmonic modes $n = 1, 2, \dots$. As a consequence, we replace $P(fl)$ by $P_n$, a discrete function of the harmonic mode number. For an individual loop, a simple (monochromatic) power spectrum can be obtained assuming the emission is dominated by specific events (e.g. cusps, kinks) (Vachaspati and Vilenkin 1985; Burden 1985; Garfinkle and Vachaspati 1987; Binetruy et al. 2009)

$$P_n = \frac{\Gamma}{\zeta(q)} n^{-q}, \qquad (61)$$

where $\zeta(q)$ is Riemann zeta function and $\Gamma = \sum_{n=1}^{\infty} P_n$ is the total power of emission, which is found to be highly peaked around $\Gamma \sim 50$ (Blanco-Pillado et al. 2015, 2014; Blanco-Pillado and Olum 2017; Wachter and Olum 2017; Blanco-Pillado et al. 2018b; Chernoff et al. 2019; Blanco-Pillado et al. 2019b). The power law parameter $q$ is 4/3, 5/3, or 2 for GW emissions dominated by cusps, kinks, or kink-kink collisions, respectively.<fnref n="19" /> With this method, the GW energy density today can be calculated by summing over loop harmonic modes (Blanco-Pillado and Olum 2017; Blanco-Pillado et al. 2014).

<footnote n="19">¹⁹ We should keep in mind that these simple power laws may not be a good approximation at low $n$. In this case, the entire loop's structure becomes important.</footnote>
$$ \frac{d\rho_{\text{GW}}}{df}(t_0, f) = G\mu^2 \sum_{n=1}^{\infty} C_n(f)P_n, \qquad (62) $$

where

$$ C_n(f) = \frac{2n}{f^2} \int_0^\infty \frac{dz}{H(z)(1+z)^4} \mathbf{n}\left(\frac{2n}{(1+z)f}, t(z)\right). \qquad (63) $$

Note that in practice, only a finite number of modes need to be included in calculations, but the necessary number of modes to ensure a good convergence for a reliable result depends on the background cosmology. For standard cosmology $\Omega_{\text{GW}}(t_0, f)$ converges by summing over $10^3 - 10^5$ modes (depending on the value of $q$) (Sanidas et al. 2012). However, recent studies (Cui et al. 2020; Gouttenoire et al. 2020a; Blasi et al. 2020) demonstrated that more than $10^5$ modes may be necessary to ensure a convergence to the correct power law at high frequency. For instance, in the presence of an early matter dominated era, summing over a small number of modes gives $f^{-1}$ at high $f$ while an $f^{-1/3}$ relation emerges with higher modes included.

### 7.3.2 Method II

The second method analytically estimates $P_{\text{GW}}(f, l)$ based on unresolved burst events on the strings (Ringeval and Suyama 2017; Abbott et al. 2018a; Siemens et al. 2007; DePies and Hogan 2007; Kuroyanagi et al. 2012, 2013; Jenkins and Sakellariadou 2018; Auclair et al. 2020). The cosmic string GW emission may be dominated by cusps, kinks, or kink-kink collisions. The contribution to the SGWB from unresolved bursts is given by

$$ \frac{d\rho_{\text{GW}}}{df}(t_0, f) = f^2 \int_{z_{\min}}^{\infty} dz \int_{0}^{\infty} dl h^2(l, z, f) \frac{d^2 R(z, l)}{dzdl}, \qquad (64) $$

where $\frac{d^2 R(z, l)}{dzdl}$ denotes the burst rate per unit loop length $l$ and per unit redshift $z$ which is proportional to $\mathbf{n}(l, t)$. $z_{\min}$ will be defined shortly. The amplitude $h(l, z, f)$ is the Fourier transform of the waveform of the bursts (Siemens et al. 2007; Damour and Vilenkin 2000, 2001; Binetruy et al. 2009) and reads

$$ h(l, z, f) = A_q(l, z) f^{-q}, \qquad (65) $$

with

$$ A_q(l, z) = g_1^{(q)} \frac{G\mu H_0 l^{2-q}}{(1+z)^{q-1} \varphi_r(z)}, \quad \varphi_r(z) \equiv H_0 \int_0^z \frac{dz'}{H(z')} dz. \qquad (66) $$

Again the power $q$ is equal to 4/3, 5/3 and 2 for cusps, kinks and kink-kink collisions, respectively. The calibration constant $g_1^{(q)}$ accounts for loop geometry effects on different bursts. Since the cusps and kinks radiate non-isotropically, the above
waveform is only valid for directions near the cusp or kink direction. A cutoff angle $\theta_{cutoff}$ is thus introduced to account for the geometric beaming effect:

$$ \theta_{cutoff}(l, z, f) = \left( \frac{1}{g_2 f (1+z)l} \right)^{1/3}, \qquad (67) $$

with $g_2 = \sqrt{3}/4$. The burst rate is then given by (Damour and Vilenkin 2000, 2001; Olmez et al. 2010)

$$ \frac{d^2 R(z, l)}{dzdl} = 2 \frac{4\pi\varphi_r^2(z)}{(1+z)^3 H(z) H_0^2} \left( \frac{n(l, t(z))}{l(1+z)} \right) \Delta(l, z, f), \qquad (68) $$

where the quantity

$$ \Delta(l, z, f) \approx \left( \frac{\theta_{\text{cutoff}}(l, z, f)}{2} \right)^{3(2-q)} \Theta(1 - \theta_{\text{cutoff}}(l, z, f)), \qquad (69) $$

represents the fraction of observable bursts (Olmez et al. 2010).

Earlier burst events that cannot be resolved at a GW detector thereby contribute to a SGWB. Integrating over $z$ and $l$, the spectrum due to a given type of burst is

$$ \Omega_{\text{GW}}(f) = \frac{\left(g_1^{(q)}\right)^2 g_2^{q-2} N_q (G\mu)^2}{2^{4-3q} 3 H_0^2} (2\pi f)^3 \int_0^\infty dx \int_{z_{\min}(x,f)}^\infty dz \frac{(ft(x))^{-2-q}}{(1+z)^{4+q}} x^{1-q} \frac{n(x)}{H(z)}, \qquad (70) $$

where $x \equiv l/t$, $z_{\min}$ is the solution to $\theta_{cutoff}(l, z_{\min}, f) = 1$, and $N_q \sim \mathcal{O}(1)$ is the average number of bursts per oscillation in a loop (Ringeval and Suyama 2017; Jenkins and Sakellariadou 2018; Abbott et al. 2018a). An alternative calculation of $\Omega_{\text{GW}}(f)$ by integrating over $z$, $h$ is shown in Cui et al. (2020).

More recent bursts of large amplitude could be resolved individually and thus detected as transient events by GW detectors, which is a different type of signal relative to a SGWB. If a burst is to be resolved in a given frequency band $f$, it must produce a strain greater than the experimental sensitivity $h > h_{\text{exp}}$ with rate less than $f$. The rate of such events is (Siemens et al. 2006; Auclair et al. 2020)

$$ R_{\text{exp}}(f) = \int_{0}^{z_{*}} dz \int_{\max(h_{\min}, h_{\text{exp}})}^{h_{\max}} dh \frac{d^2 R}{dz dh}(h, z, f), \qquad (71) $$

where $z_*$ enforces the rate condition and is given by

$$ f = \int_{0}^{z_{*}} dz \int_{h_{\min}}^{h_{\max}} dh \frac{d^2 R}{dz dh}(h, z, f). \qquad (72) $$

For a given cosmic string model, the SGWB channel is usually more sensitive than individual bursts, providing stronger constraints on parameters such as $G\mu$ (Abbott et al. 2018a). Nevertheless in certain motivated scenarios, such as string formation
before or during inflation, a SGWB could be suppressed and burst events would become the leading channel for discovery (Cui et al. 2020).

In addition to the Gaussian SGWB and the resolvable burst events as discussed above, an unresolved string signal from a set of infrequent bursts at low redshift $z \ll 1$ emits a distinct non-Gaussian and non-continuous popcorn-like signal. This effect is negligible on the current pulsar timing limit $G\mu < 10^{-11}$Regimbau et al. (2012), Olmez et al. (2010), and has no effect on the SGWB that is emitted from a large loop size $\alpha \sim 0.1$ (Blanco-Pillado and Olum 2017; Siemens et al. 2007; Olmez et al. 2010; Binetruy et al. 2012). However, if the burst occurs in our neighbourhood with strong amplitude, then GW detectors may be able to identify the signal (Chernoff and Tye 2018; Helfer et al. 2019). Such a coherent signal may provide a complementary detection mechanism: cosmic strings may generate a non-negligible amount of GW memory (Aurrekoetxea et al. 2020) that would not be captured by SGWB probes (Jenkins and Sakellariadou 2021).

### 7.3.3 Templates for the stochastic gravitational-wave background generated by cosmic string loops

Models 1 and 2 An analytical approximation for the contribution of each of three loop populations—loops that decay in the radiation era, radiation-era loops that survive into the matter era, and loops created in the matter era—to the SGWB was derived in Sousa et al. (2020). Therein, they found that radiation-era loops give rise to a SGWB of the form

$$ \Omega_{\text{GW}}^r(f) = \frac{128}{9} \pi A_r \Omega_r \frac{G\mu}{\epsilon_r} \left[ \left( \frac{f(1+\epsilon_r)}{B_r \Omega_m / \Omega_r + f} \right)^{3/2} - 1 \right], \quad (73) $$

where we have defined $\epsilon_r = \alpha / (\Gamma G \mu)$ and

$$ A_i = \frac{\tilde{c}}{\sqrt{2}} \mathcal{F} \frac{v_i}{\xi_i^3} \quad \text{and} \quad B_i = \frac{2H_0 \Omega_i^{1/2}}{v_i \Gamma G \mu}, \quad (74) $$

and where the labels $i=r, m$ are used to refer to the values of the corresponding variables in the radiation and matter eras, respectively. In the radiation era, we have $v_r = 1/2$, $\xi_r = 0.271$, $v_r = 0.662$ and $A_r = 5.4 \mathcal{F}$. This contribution dominates the high-frequency region of the SGWB spectrum (for $f \gg B_r \Omega_m / \Omega_r$), giving rise to a plateau of amplitude

$$ \Omega_{\text{GW}}^{\text{plateau}} h^2 = \frac{128}{9} \pi A_r \Omega_{\text{rad}} h^2 \frac{G\mu}{\epsilon_r} \left[ (\epsilon_r + 1)^{3/2} - 1 \right] \approx 1.02 \times 10^{-2} \frac{G\mu}{\epsilon_r} \left[ (\epsilon_r + 1)^{3/2} - 1 \right]. \quad (75) $$

The contribution of the loops that are created in the radiation era but survive into the matter era is of the form
$$ \Omega_{\text{GW}}^{m}(f) = 32\sqrt{3}\pi(\Omega_{m}\Omega_{r})^{3/4}H_{0}\frac{A_{r}(\epsilon_{r}+1)^{3/2}}{\Gamma f^{1/2}\epsilon_{r}} \left\{ \frac{\left(\frac{\Omega_{m}}{\Omega_{r}}\right)^{1/4}}{\left(B_{m}\left(\frac{\Omega_{m}}{\Omega_{r}}\right)^{1/2}+f\right)^{1/2}} \left[ 2 + \frac{f}{B_{m}\left(\frac{\Omega_{m}}{\Omega_{r}}\right)^{1/2}+f} \right] - \frac{1}{(B_{m}+f)^{1/2}} \left[ 2 + \frac{f}{B_{m}+f} \right] \right\}, \quad (76) $$

and gives rise to a peak in the low-frequency portion of the spectrum (for $f \ll B_r \Omega_m / \Omega_r$).

Matter-era loops also give rise to a peak-like contribution in the same frequency range, albeit with a different shape:

$$ \Omega_{\text{GW}}^{m}(f) = 54\pi H_{0} \Omega_{m}^{3/2} \frac{A_{m}}{\Gamma} \frac{\epsilon_{m} + 1}{\epsilon_{m}} \frac{B_{m}}{f} \left\{ \frac{2B_{m} + f}{B_{m}(B_{m} + f)} - \frac{1}{f} \frac{2\epsilon_{m} + 1}{\epsilon_{m}(\epsilon_{m} + 1)} + \frac{2}{f} \log \left( \frac{\epsilon_{m} + 1}{\epsilon_{m}} \frac{B_{m}}{B_{m} + f} \right) \right\}, \quad (77) $$

where $\epsilon_m = \epsilon_r \xi_m / \xi_r$ and we have $v_m = 2/3$, $\xi_m = 0.625$, $v_m = 0.583$ and $A_m = 0.39\mathcal{F}$. The shape of the SGWB generated by cosmic string loops in this frequency range is then determined by the interplay between these two contributions. For larger values of $\alpha$, $\Omega_{\text{GW}}^{rm}$ dominates, although its relative importance decreases as $\alpha$ decreases (since $\Omega_{\text{GW}}^{rm} \sim \alpha^{1/2}$). So, for small enough $\alpha$, the contribution of matter-era loops dominates the SGWB spectrum in the low-frequency range.

It was demonstrated in Sousa et al. (2020) that the SGWB is well described by an approximation of the form

$$ \Omega_{\text{GW}}(f) = \Omega_{\text{GW}}^{r}(f) + \Omega_{\text{GW}}^{rm}(f) + \Omega_{\text{GW}}^{m}(f) \quad (78) $$

for $\alpha \gtrsim \Gamma G\mu$ and $f < 3.5 \times 10^{10}/(1 + \epsilon_r)$ Hz (larger frequencies are outside the sensitivity windows of the current major GW experiments; see however (Sousa et al. 2020) for a description of how to extend these results for larger frequencies). For $\alpha \lesssim \Gamma G\mu$, loops survive significantly less than a Hubble time and decay effectively immediately on cosmological times scales. In this case, no loops produced in the radiation era are expected to survive into the matter era and, as a result, the contribution $\Omega_{\text{GW}}^{rm}$ should be switched off. However, for small enough $\alpha$ (i.e. $\alpha \lesssim 0.1\Gamma G\mu$), this approximation further simplifies to (Sousa and Avelino 2014; Sousa et al. 2020)

$$ \Omega_{\text{GW}}(f) = \frac{64\pi}{3} G\mu \Omega_{r} A_{r} + 54\pi \frac{H_{0}\Omega_{m}^{3/2}}{\epsilon_{m}\Gamma} \frac{A_{m}}{f} \left[ 1 - \frac{B_{m}}{\epsilon_{m}f} \right]. \quad (79) $$

Note that, since Model 1 may be calibrated to describe Model 2 (Auclair et al. 2020), Eq. (78) can also be used as a template for this model by setting $\mathcal{F} = 0.1$ and $\alpha = 0.1$.<fnref n="20" /> Note also that, although these approximations only provide a description of the SGWB generated by the fundamental mode of emission, one may use it to

<footnote n="20">²⁰ As discussed in Sect. 7.2.2, for $G\mu \lesssim 10^{-10}$, the contribution from matter-era loops $\Omega_{\text{GW}}^{m}$ is negligible and does not need to be included.</footnote>
construct an analytical approximation up to an arbitrary number of modes of emission $n_*$ for any $q$ using

$$ \Omega_{\text{GW}}(f, q, n_*) = \sum_{j=1}^{n_*} \frac{j^{-q}}{\mathcal{E}} \Omega_{\text{GW}}(f/j), \qquad (80) $$

with $\mathcal{E} = \sum_m n_* m^{-q}$ (Fig. 10).

Model 3 A first review of the SGWB of Model 3 was conducted in Ringeval and Suyama (2017), and we present below analytical approximations derived in Auclair et al. (2020). The loop distribution of Model 3 is similar on large scales to the two above-mentioned models, but is also characterised by an extra population of small loops with invariant lengths smaller than $\Gamma G\mu t$. This extra population is treated separated from the rest of the loop distribution resulting in five contributions, i.e. the extra populations of small loops emitting GWs either during the radiation era or the matter era, and the three populations mentioned in the previous section.

The contribution from radiation-era loops is given by

$$ \Omega_{\text{GW}}^r = \frac{64\pi C_r \Omega_r}{3\Gamma(2 - 2\chi_r)} (\Gamma G\mu)^{2\chi_r} \left(1 + \frac{4H_r(1+z_{\text{eq}})}{f\Gamma G\mu}\right)^{2\chi_r - 2}, \qquad (81) $$

where $C_r = 0.08$ and $\chi_r = 0.2$. This formula is characterised by a plateau at high frequency. The loops produced during radiation-era that survive during matter era also contribute to the SGWB

$$ \Omega_{\text{GW}}^{rm} = \frac{54\pi C_r H_m \Omega_m}{\Gamma f (\Gamma G \mu)^{1-2\chi_r}} (1+z_{\text{eq}})^{3(2\chi_r-1)/2} \left[ \frac{x^{2-6\chi_r}}{2-6\chi_r} {}_2F_1(3-2\chi_r, 2-6\chi_r; 3-6\chi_r; -\frac{3H_m x}{f\Gamma G\mu}) \right]_1^{\sqrt{1+z_{\text{eq}}}}. \qquad (82) $$

Matter-era loops give rise to a SGWB given by

$$ \Omega_{\text{GW}}^{m} = \frac{2 \times 3^{2\chi_m} \pi C_m \Omega_m}{H_m^{2-2\chi_m} \Gamma f^{2\chi_m-2}} (\Gamma G \mu)^2 \left[ \frac{x^{2\chi_m-4}}{2\chi_m-4} {}_2F_1(3-2\chi_m, 4-2\chi_m; 5-2\chi_m; -\frac{f\Gamma G\mu}{3H_m x}) \right]_1^{\sqrt{1+z_{\text{eq}}}}, \qquad (83) $$

where $C_m = 0.015$ and $\chi_m = 0.295$, and ${}_2F_1$ is the Gauss hypergeometric function.

The extra population of small loops plays an important role in determining the SGWB, particularly at high frequencies. Due to their small size ($<\Gamma G\mu t$) their contribution is cutoff at low frequencies and their lifetime is too short for them to survive from the radiation-era to the matter-era. The contribution of the extra population of small loops emitting GWs during the radiation era is given by the piece-wise formula
<figure><img src="image_14.png" /><figcaption>**Fig. 10** Analytical approximation to the SGWB generated by cosmic string networks with $G\mu = 10^{-10}$ and different values of $\alpha$ for Model 1. The solid lines represent the approximation to the SGWB in Eq. (78), while the dashed lines correspond to the SGWB obtained numerically. Image reproduced with permission from Sousa et al. (2020)</figcaption></figure>

$$ \Omega_{\text{GW}}^{\text{eps1}} = \frac{64\pi C_r \Omega_r (1/2 - 2\chi_r)}{3(1 - 2\chi_r)(2 - 2\chi_r)} G\mu_c^{2\chi_r - 1} \\ \times \begin{cases} 0 & \text{if } f < 4(1 + z_{\text{eq}})H_r(\Gamma G\mu)^{-1} \\ \left[4(1 + z_{\text{eq}})H_r/(\gamma_c f)\right]^{2\chi_r - 1} - (\Gamma G\mu/\gamma_c)^{2\chi_r - 1} & \text{if } f < 4(1 + z_{\text{eq}})H_r\gamma_c^{-1} \\ (2 - 2\chi_r) - 4(1 + z_{\text{eq}})H_r(1 - 2\chi_r)/(\gamma_c f) - (\Gamma G\mu/\gamma_c)^{2\chi_r - 1} & \text{if } f > 4(1 + z_{\text{eq}})H_r\gamma_c^{-1} \end{cases} \tag{84} $$

and the contribution of the extra population of small loops emitting GWs during the matter era is given by the piece-wise formula

$$ \Omega_{\text{GW}}^{m, \text{eps1}} = \frac{54\pi C_m H_m \Omega_m (1 - 2\chi_m)}{(3 - 2\chi_m)(2 - 2\chi_m)f} G\mu_c^{2\chi_m - 2} \left(\frac{3H_m}{\gamma_c f}\right) \\ \times \begin{cases} 0 & \text{if } f < 4H_m (\Gamma G\mu)^{-1} \\ \left[\frac{3H_m}{\gamma_c f}\right]^{2\chi_m - 3} \left\{1 - \left[3H_m/(\Gamma G\mu f)\right]^{3-2\chi_m}\right\} & \text{if } f < 4H_m \sqrt{1+z_{\text{eq}}}(\Gamma G\mu)^{-1} \\ \left[\frac{3H_m}{\gamma_c f}\right]^{2\chi_m - 3} \left[1 - (1+z_{\text{eq}})^{-(3-2\chi_m)/2}\right] & \text{if } f < 4H_m \gamma_c^{-1} \\ (3-2\chi_m)f\gamma_c/(3H_m) + (2\chi_m - 2) - \left(\frac{f\gamma_c}{3H_m\sqrt{1+z_{\text{eq}}}}\right)^{3-2\chi_m} & \text{if } f < 4H_m \sqrt{1+z_{\text{eq}}\gamma_c^{-1}} \\ \left[\frac{3H_m}{\gamma_c f}\right]^{-1} (3-2\chi_m)\left[1 - (1+z_{\text{eq}})^{-1/2}\right] & \text{if } f > 4H_m \sqrt{1+z_{\text{eq}}\gamma_c^{-1}} \end{cases} \tag{85} $$

Similarly to the previous section, these approximations assume that the GW are emitted by the fundamental mode, ie. by the oscillations of the loop. One can use Eq. (80) to construct predictions with an arbitrary number of modes (Fig. 11).

## 7.4 LISA detection prospects

Searches by current GW experiments (by LIGO/VIRGO and PTAs) on power spectra of the form $\Omega_{\text{GW}}(f) = Af^n$ have provided upper bounds on the amplitude $A$ for different fixed values of the spectral index $n$ (Abbott et al. 2019c; Lentati et al. 2015;
<figure><img src="image_15.png" /><figcaption>Fig. 11 Figure to demonstrate the validity of the templates for Model III. The dark dashed line is the result from numerical integration with only the fundamental mode and $G(z) = 1$</figcaption></figure>

Arzoumanian et al. 2018). Very recently, the NANOGrav and EPTA collaborations have even reported an excess that might be due to a SGWB signal (Arzoumanian et al. 2020; Chen et al. 2021b), although the origin is still uncertain. As no statistical significant evidence of the expected quadrupolar spatial correlations has been found, the detected signal can still be ascribed to some unknown systematics. More data is needed in order to discern whether this signal is a real SGWB or not, and if it is, whether it is due to BSM scenarios (Ellis and Lewicki 2021; De Luca et al. 2021b; Neronov et al. 2021) or to astrophysics (Sesana 2016).

The PLS curve (Thrane and Romano 2013) is a useful construction that graphically quantifies the ability of a detector to measure a SGWB with a spectrum characterised by a power-law in frequency (the extension to a generic frequency-dependence of the signal is discussed in Sect. 9). Specifically, we use the LISA PLS as introduced by Thrane and Romano (2013), but using the most updated LISA sensitivity curves based on the final configuration of LISA and new knowledge of its noise. (see LISA 2023 for all current LISA documentation.)

Claiming detection of a given SGWB from cosmic string loops (say for a given tension and other fixed string network parameters), can be roughly interpreted as the detection of the signal after 3 years of accumulated data (corresponding to 4 years of LISA operation), with a SNR $\ge$ 10. Since the shape can be more complicated than a simple power law, a more elaborated analysis following (Caprini et al. 2019; Flauger et al. 2021) is required to assess the SNR for a given detection, see also Karnesis et al. (2020). Here, for the moment, we simply quantify the parameter space compatible with a detection, without quantifying the SNR associated to such a detection. We do not reconstruct such parameter space with appropriate statistical techniques, though such work is already in the pipeline of coming work from the LISA collaboration.

In Figs. 12 and 13, we present numerical results for the SGWB generated by cosmic string loops, in a standard cosmological background. The LISA band is well
suited to set strong constraints on the string tension, thanks to the natural shift of the “bump” in the SGWB to larger frequencies as $G\mu$ decreases. This effect can be clearly seen in these figures, where we show a sequence of SGWB spectra over the LISA PLS, as we vary the string tension. We can find in this way the lowest $G\mu$ for which an intersection between the SGWB spectrum and the PLS still takes place. While the exact bound depends on our choice of model and $P_n$, in the regime LISA will probe, all three models—Model 2 and Model 3, as well as Model 1 (when calibrated to describe NG simulations by setting $\alpha = 0.1$ and $\mathcal{F} = 0.1$)—predict a string tension bound of $\mathcal{O}(10^{-17})$. We note that the trailing edge of the bump of the SGWB (that scales as $\Omega_{GW} \propto f^{3/2}$) will be the last part of the spectrum to pass through the LISA sensitivity band.

Note also that for string tensions $G\mu > 10^{-16}$, LISA will probe the high-frequency side of the SGWB bump, the particular shape of which depends on how the number of particle degrees of freedom change across the universe history. This is relevant because while the three models predict roughly equal lower bounds for the LISA window, Models 1 and 2 disagree with Model 3 at high frequencies. These discrepant regions will pass through the LISA band and hence, in the event of a detection, we could discern among the different models.

### 7.4.1 Anisotropies in the stochastic gravitational-wave background

GW sources with an inhomogeneous spatial distribution would lead to anisotropies in the SGWB, in addition to the anisotropies induced by the nature of spacetime along the line of propagation of GWs. The formalism (Jenkins and Sakellariadou 2018) to study anisotropies induced by the distribution of string loops was consequently applied for the string distribution of Models 2 and 3. It has been shown that the

<figure><img src="image_16.png" /><figcaption>**Fig. 12** Solid red curves show cosmic string SGWB curves for a range of $G\mu$ values. From the darkest most high up line to the lightest lowest one these read: $G\mu = 10^{-10}$, $G\mu = 10^{-13}$, $G\mu = 10^{-15}$ and $G\mu = 10^{-17}$. The $P_n$ used in computation of these spectra was inferred from simulations (Blanco-Pillado and Olum 2017), and the loop number density is from Model 2. The dashed orange curve shows the sensitivity of EPTA. The dark orange dash-dotted line shows the projected SKA sensitivity. The dotted black line shows the LISA PLS of SNR = 10</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Fig. 13** Identical to Fig. 12, however, $P_n \propto n^{-4/3}$ and using the loop number density from Model 3 (Lorenz et al. 2010)</figcaption></figure>

angular spectrum of the 2-point correlator is relatively insensitive to the particular choice of the string loop model: regardless of the model, the anisotropies are driven by local Poisson fluctuations in the number of loops, and the resulting angular power spectrum is spectrally white (i.e. $C_\ell = \text{constant}$ with respect to $\ell$) (Jenkins and Sakellariadou 2018).

In Fig. 14 we show the amplitude of the SGWB angular power spectrum for models 2 and 3 as a function of $G\mu$. We also include Model 1 with a small value of the initial loop size $\alpha$, as this has been studied in the literature as a source of significant anisotropies in the PTA frequency band (Kuroyanagi et al. 2017). We find that, regardless of the loop model and the string tension, the predicted $C_\ell$ spectrum is unfortunately far too small to be detected with LISA.

## 7.5 Extended and alternative considerations

In the previous sections we focused on GW signals from cosmic string loops in most conventional scenarios, i.e. NG strings evolving in the standard cosmology background with particle content of the SM. Nevertheless there are well motivated variations encoding new physics that can lead to subtle or dramatic deviations to the standard prediction of string SGWB spectra that we have demonstrated. In this section we review a few representative cases. We will start with the fact that long strings also radiate a SGWB that add up to the previously discussed SGWB from loops. We will also consider cosmic string scenarios with different loop sizes and forecast the parameter space available for exploration with LISA. We will then move on to the possible signals from networks of cosmic strings formed during a global symmetry breaking which as we will see are quite distinct as then particle production and not GW production becomes the dominant emission mode. Next we will turn to cosmic strings emerging from superstring theory and in particular the possibility of smaller reconnection probability they can feature. We finally consider the case of metastable cosmic strings.
<figure><img src="image_18.png" /><figcaption>**Fig. 14** Amplitude of the SGWB anisotropies for different cosmic string network models, as a function of the string tension. We use a representative LISA-band GW frequency of 1 mHz. Note that the spectra here are not normalised with respect to the monopole, so $\sqrt{C_l}$ is proportional to $\Omega_{GW}$. We note that the angular spectrum is constant with respect to $l$, since the anisotropies are driven by local Poisson fluctuations in the number of loops</figcaption></figure>

Another key aspect is that the general broad band nature of cosmic strings SGWB enables *cosmic archaeology* which would allow us to probe the expansion history prior to BBN as well as new particles beyond the SM that are difficult to access by other means. We will not cover this method here, instead referring the interested reader to Sect. 9 where methods of probing the expansion history are discussed and specifically to its Sect. 9.2.3 dealing with cosmic strings.

### 7.5.1 Long strings

We have focused so far on the GW emission by sub-horizon size string loops. Long strings, however, either infinite or super-horizon loops, emit also GWs. One contribution to this signal originates due to the relativistic motion of the long strings, as the network energy-momentum tensor adapts itself to maintain scaling during cosmic evolution. As a result, an emission of GWs sourced by the anisotropic stress of the network, takes place around the horizon scale at every moment (Krauss 1992; Jones-Smith et al. 2008; Fenu et al. 2009; Figueroa et al. 2013). This SGWB background is actually expected to be emitted by any scaling network of cosmic defects, independently of the topology and origin of the defects (Figueroa et al. 2013). It represents an *irreducible background* generated by any defect network (that has reached scaling). For NG cosmic string networks, such irreducible background represents however a sub-dominant signal compared to the background emitted from the loops discussed before. In the case of field-theory string networks (for which simulations to date indicate the absence of “stable” loops (Vincent et al. 1998;
Hindmarsh et al. 2009; Daverio et al. 2016; Hindmarsh et al. 2017b, 2021)), it is instead the dominant GW signal emitted by the network.

The energy density spectrum of this irreducible SGWB from long strings is predicted to be scale-invariant for the modes emitted during radiation domination (RD) (Figueroa et al. 2013). The amplitude of the background depends on the fine details of the so called unequal-time-correlator of the network’s energy-momentum tensor, which can be only accurately obtained from very large scale lattice simulations of defect networks. In the case of global defects, the scale-invariant GW power spectrum has been analytically estimated in Jones-Smith et al. (2008), Fenu et al. (2009). The amplitude of the spectrum plateau for global strings has been recently calibrated in lattice field theory simulations (Figueroa et al. 2020), obtaining

$$h_0^2 \Omega_{\text{GW}}^{(0)} \approx 3.13 \cdot 10^{-13} \left( \frac{G\mu}{10^{-6}} \right)^2. \qquad (86)$$

Using the latest Planck CMB constraints (Lopez-Eiguren et al. 2017), Figueroa et al. (2020) finds that the amplitude of the plateau satisfies

$$h^2 \Omega_{\text{GW}}^{(0)} < 9.7 \cdot 10^{-15}, \qquad (87)$$

with the number in the right hand side of the inequality corresponding to the amplitude when the CMB bound $10^{12}(G\mu)^2 \le 0.031$ is saturated. The amplitude in Eq. (87) is larger, for instance, than the maximum amplitude expected (as bounded by current CMB constraints (Ade et al. 2021)) for the quasi-scale invariant GW background in slow-roll inflation (Caprini and Figueroa 2018), $h_0^2 \Omega_{\text{GW}}^{(\text{inf})} \lesssim 10^{-16}$. It is however, still too small to be observed by LISA, and it is clearly subdominant when compared to the amplitude of the dominant GW signal from the long lived NG loops, which scales as $(G\mu)^{1/2}$.

One can also consider the contribution to the GW spectrum coming from the accumulation of small-scale structure on long strings. These kinks are the product of the multiple intercommutations that infinite strings suffer over the course of their cosmological evolution, and were noticed early on in numerical simulation of cosmic networks (Bennett and Bouchet 1988; Sakellariadou and Vilenkin 1990). The emission of GW from individual infinite strings modulated by kinks has been calculated in Sakellariadou (1990), Hindmarsh (1990). Using these results, one can also compute the spectrum produced by these kinks on a network assuming the simple model in which their characteristic scale is given by $\alpha t$. At high frequencies one can then estimate that the radiation-era plateau of this contribution should be (Vilenkin and Shellard 2000)

$$h^2 \Omega_{\text{GW}} \approx \frac{128\pi^2}{3\xi^2\alpha} h^2 \Omega_{\text{rad}} (G\mu)^2, \qquad (88)$$

which, for $\alpha \approx 0.1$ and $\xi_r = 0.271$, shows a rough agreement with the value obtained from field theory simulations. On the other hand, recently, Matsui and Kuroyanagi (2019) has calculated the GW spectrum produced by kink-kink collisions on long strings, and found that the amplitude is larger than in previous estimates. This is
because the characteristic scale $\alpha$ turns out to be much smaller than 0.1 according to their semi-analytic estimation of the kink number distribution.

### 7.5.2 Agnostic approach to loop size

In Sect. 7.4, we have analyzed the detection prospects for NG strings using either Model 1 (calibrated to NG simulations) or simulation-inferred Models 2 and 3. Here, we extend the analysis a bit further, by using Model 1 to study scenarios with different loop sizes. As we have seen in Sect. 7.3.3, although the typical shape of the SGWB spectrum generated by cosmic strings is roughly independent of $\alpha$, the amplitude of the radiation-era plateau and the characteristics of the peak—its shape, height and broadness—are highly dependent on the size of loops. As a matter of fact, the amplitude of the spectrum generally decreases with decreasing $\alpha$ and, therefore, LISA should, in general, be less sensitive to scenarios in which loops are created with a smaller size.

The $(\alpha, G\mu)$ parameter space available for exploration with LISA—which we plot in Fig. 15—is characterized in Auclair et al. (2020). This analysis finds that LISA should be able to probe about 16 orders of magnitude in loop size and, thus, it will have a good capability to detect string models that deviate from the standard NG scenario. Note, however, that the projected constraints on cosmic string tension are also less stringent as the loop size decreases. Nevertheless, LISA shall be able (conservatively) to probe cosmic string scenarios in which loop production is significant up to tensions

$$G\mu < 8 \times 10^{-12}, \qquad (89)$$

independently of loop size (Auclair et al. 2020). This corresponds to an improvement of almost 5 order of magnitude over current equivalent constraints (Sanidas et al.

<figure><img src="image_19.png" /><figcaption>Fig. 15 Projected constraints on $G\mu$ of the LISA mission for cosmic string scenarios characterised by different loop-size parameter $\alpha$ for $n_* = 1$ (dashed line) and $n_* = 10^5$, with $q = 4/3$ (dash-dotted line). The shaded area corresponds to the region of the $(\alpha, G\mu)$ parameter space that will be fully available for exploration with LISA. The dotted line corresponds to scenarios for which $\alpha = \Gamma G\mu$, so that the region above this line corresponds to cosmic string models in which loops are small, while the region below corresponds to the large loop regime. Image reproduced with permission from Auclair et al. (2020), copyright by IOP/SISSA</figcaption></figure>
2012) and of more than 2 orders of magnitude over the projected SKA constraints (Sanidas et al. 2013).

### 7.5.3 Global strings

Global cosmic strings are generically predicted in BSM frameworks with e.g. a post-inflationary global U(1) symmetry breaking which may be associated to axion-like DM. The GW emission from global strings have been less actively studied than that from NG strings, with only few papers having addressed an explicit computation of the SGWB from global/axion strings (see Jones-Smith et al. 2008; Fenu et al. 2009; Giblin et al. 2012; Figueroa et al. 2013; Chang and Cui 2020; Gouttenoire et al. 2020a; Figueroa et al. 2020; Gorghetto et al. 2021). Studying this topic is however particularly timely given the increasing interest in axion DM. In the following we review the semi-analytical approach taken in the recent work (Chang and Cui 2020).

The approach taken in Chang and Cui (2020) systematically follows a similar procedure as for NG strings which is reviewed in Sect. 7.2.1 but there are a few key differences for global strings. First, global strings have a time-dependent string tension (Copeland et al. 1990; Dabholkar and Quashnock 1990; Vilenkin and Shellard 2000)

$$ \mu(t) = 2\pi\eta^2\ln(L/\delta) \equiv 2\pi\eta^2 N, \qquad (90) $$

where $L \simeq H^{-1}\xi^{-1}$ is the string correlation length, $\xi$ is the number of long strings per horizon volume, $\delta \simeq 1/\eta$ is the string thickness, and $N \equiv \ln(L/\delta) \simeq \ln(\eta\xi^{-1}t)$ is time-dependent. Once reaching the scaling regime, the long string energy density evolves as

$$ \rho_{\infty} = \xi(t) \frac{\mu(t)}{t^2}, \qquad (91) $$

where $\xi(t)$ quickly approaches a constant for NG strings, yet needs to be determined for global strings. In addition, global string loops have an additional significant decay channel, through Goldstone emission. Taking these into account, the evolution equations of a global string network are as follows based on a Velocity-dependent One-Scale model (Martins 2019; Vilenkin and Shellard 2000):

$$ \left(2 - \frac{1}{N}\right) \frac{dL}{dt} = 2HL(1 + v^2) + cv + \sigma \frac{v}{N}, \quad \frac{dv}{dt} = (1 - v^2) \left[ \frac{k}{L} - 2Hv \right], \qquad (92) $$

where $k \sim 0.28$ is a momentum parameter, $v \sim 0.57$ is the averaged long string velocity, $c \sim 0.5$ is the loop chopping parameter and $\sigma \sim 5.83$ is the Golstone radiation parameter. These values are obtained (Chang and Cui 2020) from recent simulation results (Gorghetto et al. 2018; Kawasaki et al. 2018; Hindmarsh et al. 2020; Klaer and Moore 2017).

As presented in Eq. (90), reliable numerical simulations of global string networks require a huge scale separation in string length $L \propto 1/H(t)$ and core width $\delta$. Thus the global string simulation is much more challenging. Nevertheless, it has seen rapid
development in recent years (Gorghetto et al. 2018; Kawasaki et al. 2018; Hindmarsh et al. 2020; Klaer and Moore 2017; Buschmann et al. 2020; Figueroa et al. 2020), while uncertainties remain to be resolved with future higher resolution simulations. Due to the lack of dedicated simulation results for global string loop distribution at formation time, a reasonable benchmark inspired by Model 2 for the NG case is considered while the effects of alternative possibilities of loop distribution are discussed in Chang and Cui (2020).

Once formed, a loop oscillates and loses energy by the rate (Vilenkin and Vachaspati 1987; Battye and Shellard 1996, 1994; Vilenkin and Shellard 2000)

$$dE/dt = -\Gamma G\mu^2 - \Gamma_a \eta^2, \qquad (93)$$

where the right hand side represents GW and Goldstone radiation in order. Studies show that $\Gamma \simeq 50$ (Vilenkin 1981; Blanco-Pillado et al. 2011, 2014; Blanco-Pillado and Olum 2017), $\Gamma_a \simeq 65$ (Vilenkin and Shellard 2000). Consequently the length of a loop after its formation time $t_i$ would evolve as

$$\ell(t) \simeq \alpha t_i - \Gamma G \mu (t - t_i) - \kappa (t - t_i), \qquad (94)$$

where $\kappa \equiv \Gamma_a/(2\pi N)$. The SGWB from global strings can be computed following the similar procedure of computing the SGWB for NG strings while taking into account these distinctions. Figure 16 illustrates the results with varying symmetry breaking scale $\eta$, assuming standard cosmic history, with the NG string GW spectra shown in contrast. As can be seen, the global string amplitudes are more sensitive to scale $\eta$, i.e. $\Omega_{GW}^{\text{global}} \propto \eta^4$ and the NG string $\Omega_{GW}^{\text{NG}} \propto \eta$. The spectrum falls off at $f \sim \frac{1}{t_0} \sim 10^{-17}$ Hz due to the last GW emissions, today. Meanwhile, in the late matter domination period, it follows the frequency dependence $\Omega_{GW}^{\text{global}}(f) \propto f^{-1/3}$ for summation of higher oscillation normal modes $n \gg 10^5$ as demonstrated in very recent studies (Cui et al. 2020; Blasi et al. 2020). Then, at $f \sim 10^{-8}$ Hz, the logarithmic time dependence of $\mu$ made a gradually, logarithmically declining plateau towards high $f$ on the GW spectrum, instead of a large plateau in NG strings in RD period.

<figure><img src="image_20.png" /><figcaption>Fig. 16 SGWB spectrum from a global (solid) and gauge (dashed) string network in standard cosmology with loop length parameter $\alpha = 0.1$, and symmetry breaking scale $\eta = 10^{15}, 5 \times 10^{14}, 10^{14}$ GeV for red, green, blue, respectively. The grey region is the LISA sensitivity</figcaption></figure>
In very recent simulation studies (Saurabh et al. 2020), they have concluded that the global string loop lifetime is of the order the loop initial length $l \sim \alpha t$ which is in a good agreement with the prediction of Eq. (94) by energy conservation. However, the amplitude of the SGWB from large-scale simulations of global string networks from Figueroa et al. (2020), Figueroa et al. (2013), c.f. Eqs. (86) and (87), obtain an amplitude of the SGWB notably smaller than the above semi-analytical prediction. More work is therefore needed, with further investigation and improvements on both global string simulation and analytical analysis, in order to understand the origin of this discrepancy.

### 7.5.4 Cosmic superstrings

Several brane-inflationary scenarios predict the copious production of cosmic superstrings (see e.g. Copeland et al. 2004; Dvali et al. 2001; Jones et al. 2003): fundamental strings—or F-strings—and 1-dimensional Dirichlet branes—or D-strings—that grow to macroscopic sizes and play the cosmological role of cosmic strings. These cosmic superstrings, as a result of their quantum nature, do not always exchange partners when they collide: strings may simply pass through each other without intercommutation. In other words, the probability of intercommutation $P$ is —unlike that of ordinary strings (Shellard 1987; Verbiest and Achucarro 2011)— smaller than 1. In fact, it was shown that $10^{-3} \lesssim P \lesssim 1$ in collisions between F-strings and $10^{-1} \lesssim P \lesssim 1$ for D-string collisions (Jackson et al. 2005).

Cosmic superstring networks are then expected to lose energy less efficiently and, consequently, to be significantly denser than networks of ordinary cosmic strings and to generate a SGWB with a higher amplitude as a result. As a matter of fact, the loop-chopping parameter—which quantifies the energy that is lost in the form of loops—is expected to be such that $\tilde{c} \propto P^\gamma$, where $\gamma > 0$.<fnref n="21" /> For such weakly interacting networks, with $\tilde{c} \ll 1$, the amplitude of the SGWB is expected to roughly scale as (Avelino and Sousa 2012; Sousa and Avelino 2016)

$$ \Omega_{GW} \propto \tilde{c}^{-2} \propto P^{-2\gamma}. \qquad (95) $$

The constraints on $G\mu$ derived on the previous sections for $P=1$ are then conservative for cosmic superstrings: for $P<1$ the bounds should necessarily be tighter.

Note, however, that the length of loops produced by cosmic superstring networks is not known, since the precise number density of loops was not yet measured in simulations.<fnref n="22" /> Nevertheless, conservative “$\alpha$-independent” constraints on cosmic superstring tension—obtained using the radiation-era plateau of the SGWB generated by small loops—can be derived. In fact, LISA will be able to probe cosmic strings up

<footnote n="21">²¹ The exact value of the exponent $\gamma$ is still a matter of debate—$\gamma = 1/2$ was observed in NG simulations in Minkowski space (Sakellariadou 2005) and $\gamma = 1/3$ in both radiation- and matter-era simulations (Avgoustidis and Shellard 2006)—and, thus, here we discuss the effects of $P$ (mostly) qualitatively.</footnote>

<footnote n="22">²² There is some evidence that the reduction of the intercommuting probability is more efficient in suppressing the production of large loops than that of small loops (Avgoustidis and Shellard 2006), which seems to indicate that smaller $\alpha$ ($\sim \Gamma G\mu$) may be favoured for these networks.</footnote>
to a tension of $G\mu \sim 10^{-12}, 10^{-13}, 10^{-14}$, for $P = 10^{-1}, 10^{-2}, 10^{-3}$ respectively (Auclair et al. 2020).$^{23}$

Note also that there are relevant aspects of cosmic superstring dynamics that were not taken into account when deriving these constraints. In particular, when superstrings of different types collide, they are expected to bind together to create a (heavier) third type of string. This is expected to lead to networks with junctions and a hierarchy of tensions—whose dynamics differ from that of ordinary string networks (Copeland et al. 2004, 2006a, 2007, 2008; Avgoustidis and Shellard 2008; Rajantie et al. 2007; Sakellariadou and Stoica 2008; Avgoustidis and Copeland 2010; Avgoustidis et al. 2015) – and to have an impact on the shape and amplitude of the SGWB (Pourtsidou et al. 2011; Sousa and Avelino 2016). Moreover, there are several other important aspects regarding the GW emission by cosmic superstrings that need to be clarified—most notably the number and strength of the cusps (Elghozi et al. 2014) as well as the possible coupling of superstrings to other fields—before a detailed study of the parameter space available to LISA can be performed.

### 7.5.5 Metastable strings

Up to this point, the discussion in this section focused on stable cosmic strings, i.e. cosmic superstrings or strings whose stability is protected by the nontrivial vacuum topology in the underlying field theory. However, in certain scenarios, cosmic strings can become metastable, with important implications for the expected signal in GWs. In the following, we will discuss metastable strings in models with an enlarged gauge group at higher energies, as they often arise in the context of grand unified theories. For other models resulting in cosmic strings with a finite lifetime, see e.g. Kamada and Yamada (2015b), Kamada and Yamada (2015a), Bettoni et al. (2019).

Consider a $U(1)$ model giving rise to cosmic strings and its embedding in a gauge group $G$ with symmetry breaking pattern $G \to H' \times U(1) \to H$. The second symmetry breaking step will still produce cosmic strings as before. However, if $G$ and $H$ are such that the vacuum manifold of the broken subgroup $G/H$ is simply connected, these cosmic strings are metastable, i.e. classically stable, but quantum mechanically unstable against the nonperturbative Schwinger production of monopole-antimonopole pairs (Vilenkin 1982; Preskill and Vilenkin 1993; Martin and Vilenkin 1996, 1997). These monopoles are the topological defects associated with the first symmetry breaking step, $G \to H' \times U(1)$. They nucleate along the strings, such that the network breaks apart into individual string segments with a monopole on the one end and an antimonopole on the other end. These dumbbell-like objects (Nambu 1977) keep losing energy via the emission of GWs, until the defect network has fully disappeared. This late-time contribution to the GW spectrum requires further investigation. In the following, we will neglect it and focus on the emission of GWs from string loops prior to their decay.

Models based on the symmetry breaking pattern $G \to H' \times U(1) \to H$ have recently been investigated in Buchmüller et al. (2020a), Buchmüller et al. (2020b),

<footnote>$^{23}$ For larger loop sizes, the strength of the constraints may increase by up to 6 orders of magnitude (Auclair et al. 2020).</footnote>
Buchmüller (2021). While Buchmüller et al. (2020a), Buchmüller et al. (2020b) consider the breaking chain $SO(10) \to G_{SM} \times U(1)_{B-L} \to G_{SM}$, with $G_{SM}$ denoting the SM gauge group, Buchmüller (2021) studies the breaking of $G = SU(2)_R \times U(1)_{B-L}$ in an extended electroweak sector down to the hypercharge group $H = U(1)_Y$ via $H' = U(1)_R \subset SU(2)_R$. In both cases, $U(1)_{B-L}$ is the Abelian gauge symmetry associated with baryon-minus-lepton number, $B-L$, and both models thus predict the production of a network of metastable cosmic $B-L$ strings.<fnref n="24" />

The lifetime of the collapsing string network is controlled by the monopole-antimonopole nucleation rate per unit length (Leblond et al. 2009) (see also Monin and Voloshin 2008),

$$ \Gamma_d = \frac{\mu}{2\pi} e^{-\pi\kappa}, \quad \kappa = \frac{m^2}{\mu}, \qquad (96) $$

where $\mu$ denotes the cosmic string tension and $m$ is the mass of the monopoles nucleating on the string. The parameter $\kappa$ can, at least in principle, be calculated in terms of the parameters of the grand unified theory model; see Buchmüller (2021) for an example. The string network then decays around a redshift $z_d$, which can be estimated by equating the decay rate of an average-sized string with the expansion rate (Buchmüller et al. 2020a; Leblond et al. 2009),

$$ z_d \approx 8.4 \left( \Gamma G \mu \frac{\Gamma_d}{H_0^2} \right)^{1/4} . \qquad (97) $$

Here, $\Gamma \sim 50$ is defined in Eq. (61) and $H_0$ is the present-day Hubble parameter. In order to compute the SGWB signal from the collapsing string network, the lower boundary in the redshift integral in Eq. (63) needs to be replaced by $z_d$. Fixing $G\mu = 2 \times 10^{-7}$, this results in the GW spectra shown in the left panel of Fig. 17. These spectra saturate the CMB bound on $G\mu$ and are close to the sensitivity reach of existing ground-based interferometers; see the right panel of Fig. 17. At high frequencies, the SGWB signal exhibits the usual plateau, while going to lower frequencies, it decays in proportion to $f^{3/2}$. For $\sqrt{\kappa} \lesssim 6$, the string network is only short-lived, such that the signal dies off at frequencies above the LISA band, whereas for $\sqrt{\kappa} \gtrsim 8$, LISA will not be able to distinguish the signal from a standard string signal. In this case LISA will be able to access $G\mu$ values that would otherwise be excluded by PTA measurements. A LISA detection of a standard string signal at $G\mu$ in excess of current PTA bounds would therefore point to metastable strings with $\sqrt{\kappa} \lesssim 8$. In fact, values of $\sqrt{\kappa} \simeq 8$ are probed by existing PTA observations (Buchmüller et al. 2020b). Even larger values of $\sqrt{\kappa}$ can leave an imprint in CMB spectral distortions (Kite et al. 2021). The most interesting value from the perspective of LISA is $\sqrt{\kappa} \sim 7$, which can be easily achieved in realistic models (Buchmüller 2021). For $\sqrt{\kappa} \sim 7$, LISA will see the turnover in the spectrum caused by the collapse of the

<footnote n="24">²⁴ Both models also predict the production of a monopole gas during the first symmetry breaking step. These monopoles are, however, diluted by a stage of inflation preceding the production of the metastable string network and are hence irrelevant.</footnote>
<figure><img src="image_21.png" /><figcaption>Fig. 17 Left panel: SGWB signals from a metastable string network for $\tilde{G}\mu = 2 \times 10^{-7}$ and $\alpha = 0.1$, for different values of $\sqrt{k}$; see Eq. (96). The solid lines represent the numerical result, the dashed lines show the analytical approximation in terms of a flat plateau at high frequencies and an $f^{3/2}$ rise at low frequencies. Right panel: Existing and future constraints on the metastable string parameter space. The hashed vertical band marks the prediction of the $U(1)_{B-L}$ model in Buchmüller et al. (2020a), where both plots appeared for the first time. Copyright by the author(s)</figcaption></figure>

string network and hence uncover valuable information on new physics close to the energy scale of grand unification.

# 8 Inflation

Section coordinators: E. Dimastrogiovanni. Contributors: D. Comelli, E. Dimastrogiovanni, M. Fasiello, D.G. Figueroa, J. Fumagalli, L. Iacconi, A. Malhotra, S. Matarrese, A. Mazumdar, L. Pilo, S. Renaux-Petel, A. Ricciardone, R. Rollo, G. Tasinato, V. Vennin, D. Wands, L. Witkowski.

## 8.1 Introduction

Can LISA detect GWs produced at the time of the big bang? This section aims to investigate such a possibility focusing on cosmic inflation.

Inflation is the leading paradigm for describing the very early phases of cosmic expansion. Furthermore, it offers a mechanism for GW production that can be probed with GW experiments. The inflationary era is a short phase of accelerated expansion believed to have occurred within the first instants of the history of our universe and leading to the hot big bang phase, which characterises the standard cosmological evolution. This framework was initially proposed and developed to solve basic problems within the standard big bang cosmology (Guth 1981; Starobinsky 1980; Linde 1982; Albrecht and Steinhardt 1982). After the original proposal, it was soon realised that inflation comes equipped with a mechanism of particle production, an inevitable consequence of quantum mechanics applied to an accelerating cosmological spacetime. The same mechanism is responsible for generating both the primordial scalar anisotropies (Mukhanov and Chibisov 1981; Hawking 1982; Starobinsky 1982; Guth and Pi 1982; Bardeen et al. 1983) that source the evolution
of cosmic structures at large scales, and a SGWB of primordial GWs (Grishchuk 1974; Starobinsky 1979; Rubakov et al. 1982; Fabbri and Pollock 1983).

The basic setup of GW production during inflation is as follows. Consider the dynamics of spin-2, transverse-traceless fluctuations $h_{ij}(t, \vec{x})$ around a flat FLRW line element

$$ds^2 = -dt^2 + a^2(t) [\delta_{ij} + h_{ij}(t, \vec{x})] dx^i dx^j, \qquad (98)$$

with $a(t)$ being the scale factor. During inflation one typically assumes that the Hubble parameter $H = \dot{a}/a$ is nearly constant, leading to a quasi-exponential expansion. The dynamics of spin-2 fluctuations is obtained by expanding the Einstein-Hilbert action at quadratic order in $h_{ij}$ around the homogeneous FLRW solution above. Cosmological perturbations can be quantized according to the basic rules of quantum field theory (see e.g. Birrell and Davies 1984). Quantum mechanics converts the large gradients characterising the spacetime geometry, associated with the rapid cosmological expansion, into the production of spin-2 quanta whose amplitude, at wavelengths larger than the Hubble horizon $H^{-1}$, freezes to a nearly scale-invariant value

$$P_h(k) \approx \frac{2}{\pi^2} \left(\frac{H_{\text{inf}}}{M_{\text{Pl}}}\right)^2 \left(\frac{k}{k_*}\right)^{n_T}, \qquad (99)$$

with $k_*$ a pivot scale, $H_{\text{inf}}$ is the inflationary Hubble scale when $k_*$ exited the Hubble radius during inflation and $P_h$ is the dimensionless tensor power spectrum (here $M_{\text{Pl}}$ is the reduced Planck mass). In an exact de Sitter background, the tensor spectrum would be exactly scale invariant, $n_T = 0$. In slow-roll inflation, however, a small (“slow-roll suppressed”) red tilt is developed, $n_T \approx -2\epsilon$, where $\epsilon \equiv \frac{M_{\text{Pl}}^2}{2} \left(\frac{1}{V} \frac{dV}{d\phi}\right)^2$ is a slow roll parameter (and $V$ is the potential of the inflaton field $\phi$).

After inflation ends, the Hubble scale starts to increase, reaching the size of the wavelength of the fluctuations produced during inflation. The latter re-enter the horizon with a very large occupation number, behaving as classical stochastic variables, and forming a primordial SGWB (see e.g. Caprini and Figueroa 2018; Maggiore 2018 for reviews). A detection of the (as of today undetected) inflationary SGWB would provide key information on the physics of the early universe. Moreover, it would be the first direct experimental evidence of the quantization (see Riotto 2003 for a review and Polarski and Starobinsky 1996 for a discussion on the transition to a semi-classical behaviour) of spin-2 gravitational interactions: the production of GWs is based on quantum mechanical notions applied to cosmology.

The simplest class of inflationary models relies on a single scalar field, the inflaton, characterised by a standard kinetic term and driving the acceleration. In this configuration, the inflaton homogeneous and time-dependent profile rolls slowly down an almost flat potential throughout the inflationary phase (Liddle and Lyth 2000; Baumann 2011). Such scenarios are associated to rather specific predictions in terms of the inflationary SGWB properties. They consist in an almost scale-invariant spectrum, whose amplitude is too small for detection by direct GW experiments such as LISA and whose slope is slightly red-tilted (the amplitude decreases as frequency
increases). At present the most promising path to detection of an inflationary SGWB from standard single-field slow-roll models is through its footprints on the B-mode polarisation of the CMB radiation (see e.g. Kamionkowski and Kovetz 2016 for a review).

There are nevertheless plenty of good reasons to go beyond vanilla models of inflation. Such a step is motivated for example by the need to embed inflationary physics in the broader context of particle physics and especially within quantum theories of gravity. Upon exploring a richer inflationary dynamics, one must also revisit the corresponding predictions of primordial GW properties.

Particularly relevant for this work is the fact that in scenarios beyond single-field slow-roll, the amplitude of the primordial SGWB can be enhanced, and this may well occur in the LISA frequency band. As a result, LISA has the potential and the opportunity to probe our understanding of the early universe and, in particular, of the inflationary phase. In this section we expand on such a possibility. We refer the reader to Bartolo et al. (2016b) for previous work on the topic by the LISA Cosmology Working Group.

We shall consider here the following scenarios:

* When embedding inflation in particle physics and string theory setups (Lyth and Riotto 1999; Baumann and McAllister 2015), several additional fields may enter the game thereby changing the basic predictions of inflation. As we shall see in Sect. 8.2, the presence of extra degrees of freedom brings about the possibility of enhanced cosmological correlators and associated signatures. Of particular interest for us will be the case where such dynamics affects the tensor sector and leads to an amplification of the primordial GW signal. Additional properties of the SGWB that may reveal a multi-field mechanism include chirality, non-Gaussianity, and specific frequency profiles that can help distinguish the primordial signal from its astrophysical counterpart.

* A SGWB with sufficiently large amplitude to allow detection by LISA can be induced at second order in perturbations by a boosted scalar power spectrum. The mechanism leading to such enhancement typically gives rise also to oscillatory features. In Sect. 8.3 we discuss the characteristic oscillatory frequency profiles of the SGWB that can be probed in the LISA frequency ranges. In the process, we discuss novel methods to identify the different signatures associated to a number of inflationary models.

* The inflationary mechanism can be realised in scenarios with alternative spacetime symmetry breaking patterns, in which the background of the fields driving inflation depends also on space coordinates. As we discuss in Sect. 8.4, the dynamics of the corresponding tensor sector are different than in standard scenarios, and can lead to a primordial SGWB signal within reach of LISA.

* Explicit UV completions of inflation motivated by non-local versions of Starobinsky's $R + \alpha R^2$ model are discussed in Sect. 8.5. There we also show how considering such approaches reflects on the properties of the primordial GW spectrum.
* Inflation is an early period of accelerated expansion that must come to end so as to give way to cosmological evolution within the framework of standard big bang cosmology. In the transition phase the inflaton field couples with or decays into SM particles. This process may lead to GW production and to a SGWB spectrum whose properties are investigated in Sect. 8.6.

* All the scenarios mentioned above point to the fact that the primordial SGWB from inflation can be produced by qualitatively distinct mechanisms, all of which are markedly different from astrophysical processes. Section 8.7 aims to put together and investigate the features of the inflationary SGWB—chirality, frequency profile, anisotropies—that can be probed by LISA.

## 8.2 Gravitational waves sourced by additional fields during inflation

### 8.2.1 Axion-gauge field inflation

Ubiquitous in particle physics, the existence of axion-like particles driving inflation has received considerable attention (see Pajer and Peloso 2013 for a comprehensive review) starting with the well-known natural inflation proposal (Freese et al. 1990; Adams et al. 1993). The appeal of such a setup relies in part on their approximately shift-symmetric potential protecting the inflaton mass from large corrections. As a result of non-perturbative contributions from gauge field configurations (instantons), the potential acquires the typical cosine profile

$$V(\chi) = \Lambda^4 \left[ 1 + \cos\left(\frac{\chi}{f_{\chi}}\right) \right]. \qquad (100)$$

The axion-like field $\chi$ is driven by a potential at a scale $\Lambda$ and with decay constant $f_{\chi}$. In the absence of other fields, a nearly scale-invariant scalar spectral index, as required by CMB observations, leads to a large axion decay parameter $f_{\chi} \sim M_P$. On the other hand, a number of reasons favour inflationary realisations with a sub-Planckian $f_{\chi}$: (i) quantum gravity effects are expected to break the shift symmetry, as well as any global symmetry, at the Planck scale through the formation of a (virtual) BH; (ii) in typical string theory constructions one finds $f_{\chi} < M_P$ (Banks et al. 2003). It follows that a field content allowing for a smaller decay constant whilst preserving the naturalness of the potential, is of interest in this context. Coupling the axion-inflaton with gauge fields is a simple and intriguing possibility. Perhaps the most studied example is that of a Chern–Simons-type coupling

$$\mathcal{L} \supset \frac{\lambda}{4f_{\chi}} \chi F \tilde{F}, \qquad (101)$$

which preserves the shift symmetry and “dissipates” some of the inflaton kinetic energy into the gauge sector, in so flattening its effective potential without the need to resort to a large $f_{\chi}$. Several specific axion-gauge fields realisations are found in the literature, employing both Abelian (Anber and Sorbo 2010) and non-Abelian (Adshead and Wyman 2012) gauge modes. Their phenomenology is a rich one, but
we shall focus here mostly on the gravity (i.e. tensor) sector. The Abelian and non-Abelian scenarios share some key features:

• Inspection of the equations of motion for the gauge sector in Eqs. (102) and (103), reveals that the effect of the coupling to the axion-like inflaton depends also on the popularisation of the gauge field, and it is controlled by the parameter $\xi \equiv \lambda\dot{\chi}/(2f_{\chi}H)$. More specifically, denoting by $A_{\pm}$ the mode functions of the two circular polarisations of the gauge field, one has (Garretson et al. 1992; Anber and Sorbo 2010; Maleknejad and Sheikh-Jabbari 2011; Adshead and Wyman 2012)

$$ \text{Abelian case :} \qquad A_{\pm}'' + \left( k^2 \pm \frac{2k}{\tau} \xi \right) A_{\pm} = 0; \qquad (102) $$

$$ \text{non-Abelian case:} \quad t_{\pm}'' + \left[ k^2 + \frac{2}{\tau^2} (m_Q \xi \pm k \tau (m_Q + \xi)) \right] t_{\pm} \approx 0; \qquad (103) $$

where $\tau$ is conformal time. The quantity $t$ in Eq. (103) is the wave-function associated with the tensor degrees of freedom $t_{ia}$, stemming from choosing SU(2) as the gauge group (Maleknejad and Sheikh-Jabbari 2013); the parameter $m_Q$ is defined as $m_Q \equiv gQ/H$ where $g$ is the gauge coupling, $H$ the Hubble rate, and $Q$ the background of one of the (three) gauge sector scalar modes. We note the slow roll relation $\xi \simeq m_Q + 1/m_Q$, so that the two parameters that control the production in the Abelian and non-Abelian case coincide in the large $\xi$ limit.

It is intuitively clear then that the effect of the Chern–Simons coupling, driven by $\xi$, is stronger as one approaches the end of inflation where a large kinetic term $\dot{\chi}$ breaks the slow-roll condition.

• Given that the amplitude of the relevant gauge degrees of freedom is enhanced by the coupling, and that these modes source GWs, the tensor power spectrum is typically blue in these scenarios. This makes them of immediate interest for GW detectors, such as LISA (see e.g. Figs. 18 and 19, respectively for the Abelian and for the non-Abelian cases).

The fact that minimal single-field slow-roll scenarios predict instead a slightly redtilted GW power spectrum, well below the sensitivity of LISA and possibly of the proposed BBO, is very much relevant here. A detection of a primordial GW signal at small scales will strongly point to a multi-field inflation mechanism such as the ones discussed in this and the next subsection (Fig. 20).

• Polarisation-dependent equations of motion and solutions for gauge fields feed directly into a polarisation dependent GW spectrum. If the sourced contribution is comparable or larger than the ever-present vacuum tensor fluctuations the GW spectrum will be a chiral one. This is particularly interesting a signature because (i) it is testable with LISA by means of the kinematically induced dipole (Domcke et al. 2020b) and (ii) is rather uniquely associated with a Chern–Simons-type coupling, at least in terms of chirality of primordial origin.

Further characterisation of these inflationary scenarios necessarily includes the study of their scalar sector, especially vis a vis CMB constraints as well as scalar, tensor and mixed higher-point functions. The possibility of large tensor non-Gaussianity (i.e. the GW bispectrum) (Anber and Sorbo 2012; Namba et al. 2016; Agrawal et al. 2018) is
<figure><img src="image_22.png" /><figcaption>**Fig. 18** Axion-Abelian gauge field coupling can produce a large amplitude and varying width GW background depending on the evolution of fields during inflation. One can have a bump/Gaussian peak at scales where the axion rolls fast. The figure shows three examples that are relevant for current and forthcoming/future interferometers, namely PTA (frequencies around $10^{-10}$–$10^{-7}$ Hz), LISA (frequencies $10^{-5}$–$10^{-1}$ Hz) and Advanced LIGO-Advanced Virgo-KAGRA (around 1-100 Hz). In all those examples, the upper limit of the GW background is set by constraints from enhanced curvature perturbations and PBH bounds. Interestingly, density perturbations peaked at PTA and LISA scales can also generate PBHs which can be phenomenologically relevant and possibly constitute some fraction of DM. Image adapted from García-Bellido et al. (2017), Ünal (2019)</figcaption></figure>

worth mentioning here, with the reminder that its shape function is typically very similar to the so-called equilateral template. As a result, it is best tested at CMB scales (namely by studying the impact on the CMB of modes with wavelength not much smaller than the present horizon), rather than at higher frequencies where propagation effects strongly suppress the signal (Bartolo et al. 2019b). Moreover, models with axions and gauge fields can also produce a SGWB peaked at scales relevant for various GW detectors, Domcke et al. (2016), García-Bellido et al. (2016) including aLIGO/AVirgo/Kagra, LISA and PTAs. They can also produce peaked density perturbations that can seed an abundance of PBHs compatible with that of DM (Bugaev and Klimai 2014; García-Bellido et al. 2016; Domcke et al. 2017; García-Bellido et al. 2017), as we briefly touch upon in Sect. 8.2.3. We shall also point to the intriguing possibility of similar scenarios being responsible for the matter–antimatter asymmetry of the universe (Maleknejad 2014; Caldwell and Devulder 2018).

For the sake of completeness we should also mention some caveats one may extract from recent literature on the subject as well as several emerging new research paths. Coupling the axion with the gauge sector allows for a sub-Planckian dynamics by essentially feeding the axion kinetic energy into the gauge sector. This effect is regulated by the coupling $\lambda$. The ensuing flattening of the inflaton effective potential makes for a viable scalar spectral index $n_s$ at large scales. At the same time, the (enhanced) gauge sector also sources GWs. It turns out that the parameter space granting the appropriate $n_s$ may correspond to an overproduction of GWs (Dimastrogiovanni and Peloso 2013).

In addition, the amplified gauge fields also source scalar density perturbations, and it challenging to generate observable GW without overproducing scalar modes (Barnaby and
<figure><img src="image_23.png" /><figcaption>**Fig. 19** The GW energy density $\Omega_{GW}$ as a function of frequency sourced by SU(2) axion-gauge field setups (solid blue line and yellow dashed line). A signal outside the reach of certain probes may still be accessed via LISA at small scales (and e.g. LiteBIRD at very large scales). The vertical bars are the expected $1\sigma$ error bars estimated with and without the presence of astrophysical foregrounds (light and dark shared areas, respectively). Different colours are used for LiteBIRD (green), SKA (orange), LISA (blue), and ET (purple). Logarithmic binning in wavenumber is employed with $\Delta \ln k = 1.2$. For the sake of comparison, the tensor spectrum predicted in some illustrative single-field slow-roll scenario (dashed lines) are displayed, including their BICEP2/Keck/Planck upper bound $r = 0.06$. The quantity $k_p$ is the pivot scale and the quantity $\sigma$ in the upper left box is a parameter of the SU(2) model, typically of order 1–10. Image reproduced with permission from Campeti et al. (2021), copyright by IOP/SISSA</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Fig. 20** Feynman diagrams for the GW power spectrum. Left panel: The gauge fields $A_\mu$ source GW nonlinearly in the Abelian case. Right panel: The SU(2) scenario, where gauge tensor degrees of freedom $t_{ij}$ source GWs already at tree level</figcaption></figure>

Peloso 2011; Namba et al. 2016; Papageorgiou et al. 2018, 2019). Several routes have been explored to preserve the intriguing GW phenomenology of these models whilst also satisfying CMB and PBH bounds. One of the most economical choices is to introduce a spectator scalar field (Namba et al. 2016; Dimastrogiovanni et al. 2017; Obata and Soda 2016). This weakens the production of density perturbations and the corresponding scalar sector constraints on the axion-like particle and delivers viable scenarios both in the Abelian and non-Abelian cases. The back-reaction dynamics (Peloso et al. 2016; Maleknejad and Komatsu 2019) as well as the requirements stemming from perturbativity bounds (Ferreira et al. 2016; Peloso et al. 2016; Dimastrogiovanni et al. 2018; Papageorgiou et al. 2018, 2019) have been recently studied for these setups.
The possibility of Schwinger pair production has been explored in the same context as well as the coupling with fermionic fields (Domcke and Mukaida 2018; Mirzagholi et al. 2020b). In the same spirit, it is natural to investigate couplings with SM particles (Ferreira and Notari 2017; Domcke et al. 2020a). Also of note are several successful embeddings of axion gauge fields models within supergravity and string theory constructions (Dall'Agata 2018; McDonough and Alexander 2018; Holland et al. 2020).

### 8.2.2 Non-minimally coupled (spinning) fields

As a useful approach that is complementary to the one in the previous section, we will briefly expand here on a recently introduced EFT of non-minimally coupled (spinning) fields (Bordin et al. 2018). The presence of multiple fields during inflation is completely natural from the top-down perspective (Baumann and McAllister 2015). The benefit of an effective approach to such dynamics lies in the fact it provides a general, unified, description of the inflationary field content and the ensuing cosmological signatures. It is well-known that the presence of spinning fields makes for a richer phenomenology,<fnref n="25" /> one that is certainly worth exploring.

The presence of increasingly demanding unitary constraints as one goes up the spin ladder, typically implies that higher ($s \ge 2$) spin fields decay within of a few e-folds during inflation, to the detriment of their imprints on early universe observables. This is remedied by coupling the extra spinning content directly to the inflaton: the background breaks de Sitter isometries and weakens the related unitarity constraints, thus allowing for a lighter (more long-lived) particle content. As a result, extra spin-2 fields can for example directly source the GW spectrum providing the leading contribution to the signal (Bordin et al. 2018).

It is instructive to report a few of the leading operators in the Lagrangian coupling an extra spin-2 field $\sigma$ with the standard tensor modes $\gamma_{ij}$ and the Goldstone boson $\pi$, related to the scalar curvature $\zeta$ via $\zeta \sim -H\pi$:

$$ 
\begin{aligned}
S \supseteq & S_{\text{free}}^{(2)} + S_{\text{int}}^{(2)} + S_{\text{int}}^{(3)} = \int dt d^3 x a^3 (L_{\text{free}}^{(2)} + L_{\text{int}}^{(2)} + L_{\text{int}}^{(3)}) \\
& L_{\text{free}}^{(2)} = \frac{1}{4} ((\dot{\sigma}^{ij})^2 - c_2^2 a^{-2} (\partial_i \sigma^{jk})^2 - \frac{3}{2} (c_0^2 - c_2^2) a^{-2} (\partial_i \sigma^{ij})^2 \\
& \qquad \qquad \qquad \qquad - m^2 (\sigma^{ij})^2) \\
& L_{\text{int}}^{(2)} = - \frac{\rho}{\sqrt{2\epsilon H}} a^{-2} \partial_i \bar{\partial}_j \pi_c \sigma^{ij} + \frac{1}{2} \rho \gamma^c_{ij} \sigma^{ij} \\
& L_{\text{int}}^{(3)} = - \frac{\rho}{2\epsilon_1 H^2 M_{\text{Pl}}} a^{-2} (\partial_i \pi_c \partial_j \pi_c \dot{\sigma}^{ij} + 2H \partial_i \pi_c \partial_j \pi_c \sigma^{ij}) \\
& \qquad \qquad \qquad \qquad - \mu (\sigma^{ij})^3 + \dots,
\end{aligned}
\quad (104) $$

where $c_i$ is the sound speed for the helicity-i mode of the spin-2 particle $\sigma$. The first line of Eq. (104) contains the free quadratic Lagrangian for $\sigma$ whilst the second and third lines respectively describe the quadratic and cubic mixing of $\sigma$ with standard

<footnote n="25">²⁵ For example, scalar and mixed non-Gaussianities display a characteristic extra angular behaviour in the squeezed limit (Arkani-Hamed and Maldacena 2015). Even more relevant in this context is the fact that extra tensor modes, such as those of spin-2 particles, can directly source GWs.</footnote>
scalar and tensor modes. In addition to vacuum fluctuations, the GW power spectrum receives a contribution proportional to the quadratic coupling $\rho$, which may well be the leading one. Upon allowing a time dependent $\rho$ or a time dependent helicity-2 sound speed $c_2$, it is possible to obtain a blue GW spectrum (Iacconi et al. 2020a, b) and, more in general, one with a non-trivial scale dependence (see Fig. 21). There is a substantial portion of parameter space in this EFT that delivers a GW power spectrum detectable by LISA.

Further, even within bounds ensuring perturbativity control, lack of gradient instabilities, and compliance with CMB bounds, this setup supports large squeezed non-Gaussianities, including the purely tensorial three-point function. It follows that the EFT field content may also be tested by probing GW anisotropies (see Sect. 8.7.3) as well as via cross-correlating the SGWB with the CMB (see Sect. 8.7.4 below).

### 8.2.3 Models that produce peaks in the power spectrum of scalar perturbations

Primordial density waves in the early radiation-dominated universe interact gravitationally with GWs and therefore generate a SGWB (Mollerach et al. 2004; Ananda et al. 2007; Baumann et al. 2007). As discussed in Sect. 10.3, LISA is sensitive to GWs with frequency $f$ generated from primordial density perturbations with power spectrum $P_{\zeta}(f) \approx 10^{-4}$ which re-entered the Hubble-horizon at a temperature $T \approx 10^5(f/mHz)$ GeV. These density waves can be generated from quantum fluctuations that are swept up by the accelerated expansion during a preceding period of inflation and exit the Hubble-scale $N$ e-folds before the end of inflation. By assuming instantaneous reheating (a finite period of reheating can reduce this number), one obtains

<figure><img src="image_25.png" /><figcaption>**Fig. 21** The GW energy density of the EFT when including the contribution of an extra spin-2 field whose helicity-2 mode has a time-dependent sound speed. The two colours correspond to different values of $\rho/H$, respectively (3 $\times 10^{-3}$, 4 $\times 10^{-4}$), supporting a signal within reach for both SKA and LISA, or for LISA only. In generating both lines, the time dependence of the sound speed has been kept constant to $s_2 \equiv \dot{c}_2/(Hc_2) \approx -0.2$. Image reproduced with permission from Iacconi et al. (2020b)</figcaption></figure>
$$N \approx 22 + \ln\left(\frac{V_{*}^{1/4}}{10^{15} \text{ GeV}}\right) - \ln\left(\frac{f}{\text{mHz}}\right), \qquad (105)$$

where $N = \int H dt$ is the integrated expansion, and $V_*$ is the energy density during inflation.

It is often assumed that the power spectrum of primordial density perturbations from inflation can be described by a power law, $P_{\zeta} \propto f^{n-1}$, whose amplitude and spectral index are fixed by observations on CMB scales, where $P_{\zeta} \approx 10^{-9}$ and $n < 1$. This is not general (Chluba et al. 2015). However, within slow roll, the deviation from a simple power law (running of the spectral index) is limited (second order in slow-roll parameters) and thus a significant enhancement in the primordial power spectrum on the scales probed by LISA requires the breakdown of the usual slow-roll approximations.

Quantum fluctuations in the fields driving the accelerated expansion naturally generate primordial density fluctuations in the subsequent radiation era, and there may be additional fluctuations from fields that, although subdominant during inflation (so-called spectator fields), affect the subsequent evolution between the end of inflation and the epoch when modes re-enter the Hubble-horizon. The dimensionless density contrast, $\zeta$, can be identified with the perturbation in the local integrated expansion, $\zeta = \delta N$. In the large-scale limit (neglecting spatial gradients) this can be related to the field fluctuations at first order $\zeta = (\partial N / \partial \bar{\Phi}) \cdot \delta \bar{\Phi}$, where $N(\bar{\Phi})$ corresponds to the background expansion as a function of the local field values, $\bar{\Phi}$, during inflation.

For a single inflaton field (or more generally for adiabatic field fluctuations in a multi-field setting (Liddle et al. 1998; Gordon et al. 2000)) $\partial N / \partial \phi = H / \dot{\phi}$ while the amplitude of quantum field fluctuations at Hubble exit are given by $P_{\delta\phi} \simeq H / 2\pi c_s$, where $c_s$ is the sound speed. Given that the Hubble rate at Hubble-exit ($k = aH$) necessarily decreases during inflation, there are essentially two ways to boost the primordial density perturbations in single field models. Either one can decrease the speed of the inflaton, $\dot{\phi}$ faster than the slow-roll evolution, known as ultra-slow roll (Inoue and Yokoyama 2002; Kinney 2005), or one can rapidly decrease the sound speed, $c_s$, e.g. due to non-decoupling of massive fields (Achucarro et al. 2011). Ultra-slow roll describes a regime where the potential gradient becomes negligible and the field becomes friction dominated, $\dot{\phi} \propto e^{-3N}$, which can occur in models with a near-inflection point in the potential. Typically the ultra-slow roll phase is transient, leading to a localised step (Starobinsky 1992; Leach et al. 2001) or broad peak (García-Bellido and Ruiz Morales 2017) in the power spectrum. In some cases the departure from slow-roll may interrupt inflation, before inflation restarts in an ultra-slow roll phase (Roberts et al. 1995; Leach and Liddle 2001; Ragavendra et al. 2021).

If the energy scale of inflation is as low as $V_*^{1/4} \sim 10^6$ GeV then LISA could be sensitive to physical processes that end inflation and reheat the universe. The coherent oscillation of the inflaton (and, possibly, other field) that characterised this stage may lead to resonant amplification of scalar metric perturbations near the Hubble scale at the end of inflation, a phenomenon known as metric preheating
(Nambu and Taruya 1997; Finelli and Brandenberger 1999; Jedamzik et al. 2010; Easther et al. 2011). Density perturbations produced close to the Hubble scale at this stage may become large in models where inflation ends due to an instability, triggering a phase transition. This may be second-order, as in hybrid inflation models (Linde 1994; Copeland et al. 1994), where a tachyonic instability in a second “waterfall” field leads to the rapid growth of quantum fluctuations, or in the geometrical destabilisation of inflation, due to a geometrical instability in a negatively-curved field-space (Renaux-Petel and Turzyński 2016). A FOPT at the end of inflation requires the rapid nucleation of bubbles of a broken symmetry phase, whose collisions both reheat the universe and can generate GWs (Caprini et al. 2016).

More generally fluctuations in fields orthogonal to the instantaneous trajectory in field-space, and hence independent of the adiabatic field fluctuations, can enhance the primordial scalar power spectrum, $P_\zeta$, after Hubble exit (Wands et al. 2002). Any non-geodesic “turn” in field space during inflation converts isocurvature field perturbations at Hubble-exit into an additional source of curvature perturbations on super-Hubble scales (Gordon et al. 2000; Groot Nibbelink and van Tent 2002; Palma et al. 2020; Fumagalli et al. 2023). In a different mechanism, spectator axion fields might also exhibit a brief phase of fast roll during inflation, generating a localised peak in the power of density perturbations (Namba et al. 2016). Alternatively, the effect of spectator fields may be completely subdominant during inflation but boost the primordial power spectrum after the end of inflation as in the curvaton (Enqvist and Sloth 2002; Lyth and Wands 2002; Moroi and Takahashi 2001) or other modulated reheating (Dvali et al. 2004; Kofman 2003) scenarios. Quite generally we expect otherwise light scalar fields present during inflation to acquire effective masses of order the Hubble scale during inflation (Dine et al. 1995; Baumann and Green 2012), suppressing their fluctuations on large scales, leading to a steep blue tilt, dominating the power spectrum on small scales, as in the axion-like curvaton model (Kawasaki et al. 2013).

## 8.3 Small-scale primordial features

Embeddings of inflation in high energy theory motivate the exploration of inflationary mechanisms beyond the single-field slow-roll framework. For example, UV completions of inflation in string compactifications typically introduce many new degrees of freedom that ultimately contribute to the inflationary dynamics. Thus, from the UV point of view, single-field slow roll models have the semblance of toy-models that capture the essence of inflation without being fully realistic.

Avenues for going beyond single-field slow-roll are plentiful and, as a result, a large number of models have been constructed. However, instead of proceeding model by model, there is also a more systematic way for going beyond the simplest version of inflation. The idea is to characterise departures from single-field slow-roll in terms of their effect on the scalar power spectrum.

These signatures of a departure from single-field slow-roll have been termed “features” after corresponding characteristic properties of the scalar power spectrum.
See Chen (2010), Chluba et al. (2015), Slosar et al. (2019) for reviews on this topic. Using the language of features, departures from single-field slow-roll inflation can be assigned to one of the two following broad classes<fnref n="26" />:

1. A so-called “sharp feature” is characterised by an oscillation in the scalar power spectrum that is periodic in the wavenumber $k$. This arises whenever there is some sudden transition during inflation, like e.g. a step in the inflation potential or a sharp turn in the trajectory.

2. A so-called “resonant feature” denotes an oscillation in the scalar power spectrum in $log(k)$. It arises when some components of the background oscillate with a frequency larger than the Hubble scale, inducing a resonance with the oscillations of the quantum modes of the density perturbations, a typical example being axion monodromy inflation (Silverstein and Westphal 2008; Flauger et al. 2010).

Over the relevant range of scales, the scalar power spectrum corresponding to each class can be written as

$$ \text{Sharp: } P_{\zeta}(k) = P_0(k) \left[ 1 + A_{\text{lin}} \cos (\omega_{\text{lin}} k + \vartheta_{\text{lin}}) \right], \quad (106) $$

$$ \text{Resonant: } P_{\zeta}(k) = P_0(k) \left[ 1 + A_{\text{log}} \cos (\omega_{\text{log}} \log(k/k_{\star}) + \vartheta_{\text{log}}) \right], \quad (107) $$

with $k_{\star}$ being some arbitrary reference scale introduced for dimensional reasons. That is, a sharp or resonant feature is described by an oscillation with amplitude $A_{\text{log/lin}}$ about an envelope $P_0(k)$. The precise form on the envelope will depend on the model, but can be taken as sufficiently smooth over the period of oscillations.<fnref n="27" /> For sharp features, the frequency $\omega_{\text{lin}} \sim 1/k_{\star}$, where $k_{\star}$ corresponds to the scale that crosses the Hubble radius at the time of the feature. As for resonant ones, one has $\omega_{\text{log}} \sim M/H$, where M is the frequency of the background oscillations.

At large scales ($k \lesssim 1 \text{ Mpc}^{-1}$) features are severely constrained by CMB and LSS data, which mandate a nearly scale-invariant power spectrum with amplitude $P_{\zeta} \sim 10^{-9}$. In contrast, at small scales ($k \gg 1 \text{ Mpc}^{-1}$), CMB and LSS constraints do not apply and the scalar power spectrum can depart significantly from scale-invariance. This opens the possibility that the feature constitutes the dominant contribution to the scalar power spectrum at small scales.

For such contribution to be sufficiently large to make the SGWB detectable by current or forthcoming interferometers, the scalar power spectrum at the scale of the feature needs to be significantly enhanced compared to its value at CMB scales. Specifically, $P_{\zeta}$ should increase from $\sim 10^{-9}$ to $\sim 10^{-4}$ by moving $k$ from the CMB to the LISA momentum scale. GW observatories such as LISA are thus sensitive to

<footnote n="26">²⁶ Oscillations of heavy fields act as “primordial standard clocks”, leading to feature signals resembling that of a sharp feature at larger scales and of a resonant feature at smaller scales.</footnote>

<footnote n="27">²⁷ In realistic models the amplitude and frequency of the oscillation as well as phase offset can run with the scale $k$, but if this running is sufficiently “slow”, the templates (106) and (107) will still be applicable over a suitable range of scales.</footnote>
feature models where the scalar power spectrum exhibits a peak at small scales, where the CMB and LSS constraints do not apply, as previously explained in Sect. 8.2.3. This is not an unrealistic expectation: such localised enhancements of the power spectrum occur frequently in inflation models with features. In fact, the mechanism of amplification of scalar fluctuations is often also responsible for producing the oscillatory feature and vice versa, as observed e.g. in Ballesteros et al. (2019), Palma et al. (2020), Fumagalli et al. (2023), Braglia et al. (2020), Fumagalli et al. (2021a), Tasinato (2021), Braglia et al. (2021). GW observatories such as LISA can thus test models of inflation at scales inaccessible to CMB and LSS surveys.

What is interesting for GW astronomy is that features in the scalar power spectrum lead to corresponding features in the frequency profile of the corresponding scalar-induced GWs (Fumagalli et al. 2021a; Braglia et al. 2021; Fumagalli et al. 2021b). Here we focus on the contribution to the SGWB sourced by scalar fluctuations when they re-enter the horizon during RD after inflation, see Sect. 10.3.<fnref n="28" /> Remarkably, this inherently nonlinear effect does not lead to oscillations being washed out. One rather finds that a sharp feature leads to a periodic modulation in the spectral shape of the SGWB contribution, while a resonant feature produces a corresponding log-periodic modulation. As a result, over some range of scales the spectral shape of $\Omega_{GW}(k)$ can be matched by the following templates (Fumagalli et al. 2021a; Braglia et al. 2021; Fumagalli et al. 2021b):

$$ \text{Sharp: } \Omega_{\text{GW}}(k) = \bar{\Omega}_{\text{GW}}(k) \left[ 1 + A_{\text{lin}} \cos \left( \omega_{\text{lin}}^{\text{GW}} k + \phi_{\text{lin}} \right) \right], \quad (108) $$

$$ \text{Resonant: } \Omega_{\text{GW}}(k) = \bar{\Omega}_{\text{GW}}(k) \left[ 1 + A_{\text{log},1} \cos \left( \omega_{\text{log}} \log(k/k_{\star}) + \phi_{\text{log},1} \right) \right. \\ \left. \qquad + A_{\text{log},2} \cos \left( 2\omega_{\text{log}} \log(k/k_{\star}) + \phi_{\text{log},2} \right) \right], \quad (109) $$

with $\omega_{\text{lin}}^{\text{GW}} = \sqrt{3}\omega_{\text{lin}}$. Here, $\bar{\Omega}_{\text{GW}}(k)$ refers to the GW fraction with the oscillatory component averaged out.

The above is best illustrated by an example. As has been observed in Palma et al. (2020), Fumagalli et al. (2023), an amplification of scalar fluctuations can be achieved in multi-field inflation when the inflationary trajectory exhibits a strong turn. If the turn is also sufficiently sharp, i.e. it is executed in a short time interval, this leads to a sharp feature. In simple cases the scalar power spectrum can then be computed analytically as (Palma et al. 2020; Fumagalli et al. 2021a)

$$ \frac{P_{\zeta}(k)}{P_0} = \frac{e^{2\sqrt{(2-\kappa)\kappa\eta_{\perp}\delta}}}{2(2-\kappa)\kappa} \times \sin^2 \left( e^{-\delta/2} \kappa\eta_{\perp} + \arctan \left( \frac{\kappa}{\sqrt{(2-\kappa)\kappa}} \right) \right), \quad (110) $$

which is valid for $\kappa \le 2$ and where $\kappa \equiv k/k_{\star} = k/(k_{\text{f}}\eta_{\perp})$. Here $\eta_{\perp}$ quantifies the strength of the turn in terms of the departure from a geodesic, $\delta$ is the duration of the turn in units of e-folds and $k_{\text{f}}$ is the scale crossing the Hubble radius at the time of the

<footnote n="28">²⁸ In addition, there will also be a scalar-induced contribution to the SGWB sourced during inflation. The latter is slow-roll suppressed compared to the post-inflationary contribution and will be ignored in the following, even though it can dominate in selected models, see e.g. Zhou et al. (2020).</footnote>
sharp turn (see Fumagalli et al. 2021a for generalisations). As can easily be checked, Eq. (110) is indeed of the form of Eq. (106) characteristic of a sharp feature.

In Fig. 22a we plot $P_{\zeta}(k)$ in Eq. (110) for two example models and in Fig. 22b the corresponding GW energy density $\Omega_{GW}(k)$. We can then make the following observations, which will not only hold for the model at hand, but will be generic for a sharp feature and can be also generalised to the case of a resonant feature:

* The overall spectral shape $\bar{\Omega}_{GW}(k)$ is determined by the envelope $P_0(k)$. If this is sufficiently narrowly-peaked, the spectral shape of $\bar{\Omega}_{GW}(k)$ consists of a broad lower peak at lower frequencies and a narrow principal peak at higher frequencies. This is what is observed in Fig. 22b. The principal peak occurs at $k = 2/\sqrt{3} k_{\star}$, where $k_{\star}$ is the frequency where the envelope $P_0$ has its maximum, and this can be understood to arise from resonant amplification (Kohri and Terada 2018).

* The oscillatory running of the SGWB induced by the feature is visible as a modulation of $\Omega_{GW}(k)$ on the principal peak, see again Fig. 22b. In the vicinity of the maximum the modulation is well-matched by the template in Eq. (108) for a sharp feature, as can be seen in Fig. 22c where we plot the ratio $\Omega_{GW}(k)/\bar{\Omega}_{GW}(k)$: Especially near the centre of the peak this exhibits a sinusoidal oscillation with near-constant amplitude.

* The modulation in $\Omega_{GW}(k)$ can be understood as a superposition of resonance peaks sourced by the individual maxima of $P_{\zeta}(k)$. Labelling the maxima in $P_{\zeta}(k)$ by $k_i$, the resonance analysis predicts a series of peaks in $\Omega_{GW}(k)$ at (Cai et al. 2019c):

$$k_{\max,ij} = \frac{1}{\sqrt{3}}(k_i + k_j), \quad \text{with } k_{\max,ij} > |k_i - k_j|. \qquad (111)$$

Applying this to the sharp and resonant feature case one can show that this predicts the periodic structure recorded in the templates Eqs. (108) and (109) (Fumagalli et al. 2021a, b).

* The amplitude of oscillation in $\Omega_{GW}(k)$ is typically attenuated compared to the amplitude of oscillation in $P_{\zeta}(k)$. For example, a sharp feature with $A_{\text{lin}} = 1$ leads to a modulated GW spectrum with $A_{\text{lin}} \sim \mathcal{O}(10\%)$ (Fumagalli et al. 2021a; Braglia et al. 2021).

* Consider a feature in the scalar power spectrum with an associated peak at $k_{\star} = k_{\text{CMB}}e^{N_{\star}}$, with $N_{\star}$ the number of e-folds after horizon-crossing of the CMB modes. This bump in $P_{\zeta}(k)$ will produce a peak in $\Omega_{GW}$ whose frequency $f_{\text{peak}}$ is related to $N_{\star}$ as

$$N_{\star} \approx \ln\left(\frac{f_{\text{peak}}}{\text{Hz}}\right) + 37. \qquad (112)$$

For the peak to fall into the frequency range of maximal sensitivity of LISA, $f = 10^{-3}-10^{-2}$ Hz, the enhancement and the feature in $P_{\zeta}(k)$ has to occur
<figure><img src="image_26.png" /><figcaption>**Fig. 22** Scalar power spectrum $P_ζ$ (**a**) and the scalar-induced GW energy density $\Omega_{GW}$ (**b**) for two inflation models (see text for details) with a sharp feature and enhanced fluctuations from a strong turn in the inflationary trajectory, together with the PLS for SNR threshold $\text{SNR}_{\text{th}} = 1$ and total observation time $T = 3$ years assumed in the analysis of (Fumagalli et al. 2021a). The black dashed line in (**a**) shows the envelope of $P_ζ$ and in (**b**) the corresponding GW energy density $\bar{\Omega}_{GW}$. The $O(1)$-oscillations in $P_ζ$ are processed into $O(10\%)$ modulations on the principal peak of $\Omega_{GW}$. Over the frequency range of the principal peak, the modulations in $\Omega_{GW}$ can be modelled as cosine-oscillations about $\bar{\Omega}_{GW}$, as can be seen in (**c**). To translate $k/k_★$ into $f/\text{Hz}$, we considered $N_★ = 31.5$, i.e. the peak in $P_ζ$ occurs at wavenumber $k = k_★$, which is 31.5 e-folds larger than the CMB value. Images adapted from Fumagalli et al. (2021a)</figcaption></figure>

$N_★ \sim 30–32$ e-folds after the generation of the CMB modes, making LISA sensitive to a wide range of models producing features during the later stages of inflation.

We thus find that features in the scalar power spectrum, encoding departures from single-field slow-roll inflation, lead to corresponding features in the scalar-induced contribution to the SGWB. If the feature is associated with a sufficient enhancement of the relevant scalar fluctuations, the resulting contribution to the SGWB is in principle detectable by GW observatories like LISA. As features encode departures from single-field slow-roll without reference to explicit models, this leads to the exciting prospect that SGWB measurements can be used to learn about inflation in a model-independent fashion.
## 8.4 Effective field theory of broken space-reparametrisations

LISA can shed light also on the symmetries characterising the early universe, and in particular the inflationary period. Standard models of inflation are based on the assumption that the inflaton background field is only time dependent, and this is well described within the EFT of inflation framework (Creminelli and Vernizzi 2017). However within the standard EFT framework, it is hard to produce GWs which have an amplitude large enough to be detected by interferometers. On the other hand, if space-diffeomorphisms are broken during inflation, and so the background field can have a space-dependent vacuum expectation value, then GWs are amplified at small scales becoming a potential target for GW detectors, like LISA. In particular, the breaking of space-diffeomorphism allows the graviton to be massive, since there are no symmetries preventing tensor fluctuations from acquiring a mass during inflation. If this is the case, the graviton is characterised by an action which is the most general one for tensor fluctuations:

$$S_T^{(2)} = \frac{M_P^2}{8} \int d^4 x a^2(\eta) [h'_{ij} h^{ij} - (k^2 + a^2(\eta) m_g^2) h_{ij} h^{ij}], \quad (113)$$

where $a(\eta)$ is the scale factor and $m_g$ is the graviton mass, which is the distinctive feature of a setup that breaks space-diffeomorphism. Whereas such action can be derived within an EFT framework, there are models which realise such a possibility, like Solid (Endlich et al. 2014) and Supersolid (Koh et al. 2013; Cannone et al. 2015a, b; Bartolo et al. 2016a) inflation. They are described by Lagrangians which contain three (or four) scalar fields, which respect some internal symmetries ensuring the homogeneity and isotropy of the background. Such models can lead to a blue spectrum for tensor modes, as well as to other distinctive properties in the tensor spectrum (Endlich et al. 2014; Ricciardone and Tasinato 2017; Bartolo et al. 2016a; Cannone et al. 2015b; Akhshik et al. 2014; Akhshik 2015). Such an action leads, in the small graviton mass limit ($|m_h/H| \ll 1$), to a tensor primordial power spectrum like

$$P_h = \frac{2 H^2}{\pi^2 M_P^2} \left(\frac{k}{k_*}\right)^{n_T}, \quad (114)$$

with the tensor spectral index equal to

$$n_T = \frac{2 m_h^2}{3 H^2}. \quad (115)$$

Notice that a blue spectrum, $n_T > 0$, requires a positive $m_h^2$. This is the case of interest, since it enhances the tensor spectrum at small scales, and can lead to a signal detectable by LISA. A model-independent analysis of such kind of models, and the possibility to detect the GW signal in the range of frequencies and energy densities probed by LISA has been carried out in Bartolo et al. (2016b).
$$ \varphi^0 = \bar{\varphi}(t), \qquad \varphi^a = x^a, \quad a = 1, 2, 3; \qquad (116) $$

this time 4D diffeomorphisms are completely broken, leaving only a global symmetry as a leftover, where an SO(2) rotation is performed both to the fields and to the spatial coordinates. To allow for a FLRW background solution, the action for the scalar fields dynamics needs to be symmetric under internal rotations and shift transformations

$$ \varphi^a \rightarrow \mathcal{R}_b^a \varphi^b, \qquad \mathcal{R}^i \mathcal{R} = \mathbf{1}, \qquad \varphi^A \rightarrow \varphi^A + c^A, \qquad A = 0, 1, 2, 3. \qquad (117) $$

In general, the fluctuations around Eq. (117) can be interpreted as the phonons of a supersolid (Son 2005; Celoria et al. 2017), while the special case where only the “spatial” Stueckelberg fields $\varphi^a$ (with $a = 1, 2, 3$) are present corresponds to a solid (Endlich et al. 2013). Similar models have been proposed as massive deformation of gravity (Ballesteros et al. 2016; Celoria et al. 2018). In the scalar sector there are two propagating degrees of freedom that mix non-trivially both at early and late times and, after exiting the horizon, give rise to non-trivial cross-correlations with distinctive features for primordial non-Gaussianity. One key feature is that the production of GWs during inflation can be enhanced by the cubic interaction of the graviton with phonons (Celoria et al. 2021); the spectral index of the secondary produced GWs is blue-tilted with the amplitude within the LISA sensitivity.

## 8.5 UV complete models of $R + \alpha R^2$ model

Starobinsky's $R + \alpha R^2$ model of inflation matches the latest CMB data extremely well (Akrami et al. 2020; Starobinsky 1983). Since massless gravity contains only derivative interactions, it invites higher and infinite covariant derivative contributions as well. One particularly attractive model of higher derivative theory of gravity which contains infinite covariant derivatives and generalizes the quadratic curvature theory of gravity has been presented in Biswas et al. (2006), Modesto (2012), Biswas et al. (2012a). It has several cosmological implications, from providing the initial conditions from a non-singular cosmology to modifying the power spectrum for the primordial GWs (Koshelev et al. 2018, 2020). One particular subset of the action has been studied deeply, which generalizes Starobinsky's model and recovers it in the IR limit:

$$ S = \int d^4 x \sqrt{-g} \left( \frac{M_p^2}{2} R + \frac{\lambda}{2} \left[ R \mathcal{F}_1(\Box_s) R + W_{\mu\nu\rho\sigma} \mathcal{F}_2(\Box_s) W^{\mu\nu\rho\sigma} \right] \right). \qquad (118) $$

Here $\Box_s = \Box/M_s^2$ with $M_s$ being the scale of non-locality, and $\lambda$ is a dimensionless parameter useful to control the effect of higher curvature contributions. In the limit $M_s \to \infty$, one recovers the local Starobinsky's model of inflation. The ghost-free condition around an inflationary, i.e. approximate de Sitter background, constrains the gravitational form factors $\mathcal{F}_1$ and $\mathcal{F}_2$ (Biswas et al. 2016, 2017). Noticeably, the scalar power spectrum does not get any modification compared to the local Starobinsky's model of $R + \alpha R^2$ inflation (Starobinsky 1983). However, the
gravitational power spectrum, the tilt and the tensor to scalar ratio $r$ all get modified in an interesting manner (Koshelev et al. 2018, 2020):

$$ 
\begin{aligned}
\mathcal{P}_T &= \frac{1}{12\pi^2 \mathcal{F}_1(m^2/M_s^2)} (1 - 3\epsilon)e^{-2\gamma_T(-\bar{R}/2M_s^2)}|_{k=aH}, \\
n_T &\equiv \frac{d \ln \mathcal{P}_T}{d \ln k}|_{k=aH} \approx -\frac{3}{2N^2} - \left(\frac{2}{N} + \frac{3}{2N^2}\right) \frac{\bar{R}}{2M_s^2} \gamma_T(-\frac{\bar{R}}{2M_s^2})|_{k=aH}, \quad (119) \\
r &= \frac{12}{N^2} e^{-2\gamma_T(-\bar{R}/2M_s^2)}|_{k=aH}.
\end{aligned}
$$

The crucial difference here in comparison with the local Starobinsky model is that the tensor power spectrum is scaled by an exponential factor of $\gamma_T$ evaluated at the pole of the tensor mode $\bar{R}/(6M_s^2)$, where $\bar{R} = 3.7 \times 10^{-8}(55/N_*)^3 M_p^2$ denotes the Ricci scalar evaluated at the pivot scale $k_* = aH$, and $N_*$ denotes the number of e-foldings of inflation corresponding to the pivot scale. Accordingly, also the tensor tilt gets modified. The ghost-free condition demands that $\gamma_T$ is an entire function (namely, that it can be represented as a power series that converges everywhere in the complex plane) (Biswas et al. 2016, 2017). In the local $R^2$ model one obtains $r = 12/N^2 = 3(1 - n_s)^2$ as it follows from the original computation of scalar and tensor power spectra generated during inflation (Starobinsky 1983). From the CMB data (Ade et al. 2021) we can infer that the constraint $r < 0.036$ implies $\gamma_T > -1.05$ at $N_* = 55$, while there are no constraints on the tilt in the tensor power spectrum, due to lack of data. LISA could be able to test it, although there are indications that at LISA frequencies the nonlocal model is indistinguishable from local Starobinsky inflation (Koshelev et al. 2018; Calcagni and Kuroyanagi 2021) (section 9.3). Additional mechanisms (Bartolo et al. 2016b) or a modification in the theory (e.g., in the form factor) could reopen this opportunity.

## 8.6 Preheating

Preheating is characterised by non-perturbative particle production mechanisms (Traschen and Brandenberger 1990; Kofman et al. 1994; Shtanov et al. 1995; Kaiser 1996; Khlebnikov and Tkachev 1996; Prokopec and Roos 1997; Kaiser 1997; Kofman et al. 1997; Greene et al. 1997; Kaiser 1998), which typically take place after inflation in many models of particle physics (see Allahverdi et al. 2010; Amin et al. 2014; Lozanov 2019 for reviews). Following the end of inflation, interactions between the inflaton and some other field species—the preheat field(s)—can induce an exponential growth of the modes of the preheat field(s) within certain bands of momenta. A paradigmatic example of this is parametric resonance (Kofman et al. 1994, 1997; Greene et al. 1997; Figueroa and Torrenti 2017b), though there are other mechanisms (see below). The field gradients created during this stage can generate a sizeable anisotropic stress to source GWs, with the specific details of the resulting GW spectrum depending strongly on the considered scenario (Khlebnikov and Tkachev 1997; Easther and Lim 2006; Easther et al. 2007; García-Bellido and Figueroa 2007; García-Bellido et al. 2008; Dufaux et al. 2007, 2009; Kusenko and
Mazumdar 2008; Kusenko et al. 2009; Fenu et al. 2009; Figueroa et al. 2011). In general, preheating mechanisms are very efficient in transferring a significant fraction of the total energy available from the inflationary sector into the preheat field(s), and as a result, a large amount of GWs are radiated in the process. Furthermore, if the inflaton-preheat field coupling is tuned to certain values, the resulting SGWBs may develop large anisotropies at cosmological scales (Bethke et al. 2013, 2014).

For illustrative purposes, let us consider a standard scenario where the inflaton oscillates around the minimum of its potential after the end of inflation. For example we can consider a power-law potential $V(\phi) = \frac{1}{p}\lambda\mu^{4-p}\phi^p$, with $\lambda$ a dimensionless coefficient, $\mu$ some mass scale, and $p$ an integer index $p \ge 2$. Denoting as $t_\star$ the end of inflation, the inflaton oscillates for $t \gtrsim t_\star$ with a time-dependent frequency $\Omega_{\text{osc}} \equiv \omega_\star(t/t_\star)^{1-2/p}$, $\omega_\star \equiv \sqrt{\lambda}\mu^{(2-p/2)}\phi_\star^{(p/2-1)}$, where $\phi_\star \equiv \phi(t_\star)$ (Figueroa and Torrenti 2017b). To be specific, let us also consider a quadratic interaction $g^2\phi^2\chi^2$ between the inflaton $\phi$ and a the preheat field $\chi$, with $g$ a dimensionless coupling constant. If the resonance parameter $q_\star \equiv g^2\phi_\star^2/\omega_\star^2$ is much larger than unity $q_\star \gg 1$, the preheat field is excited through a process of broad resonance during the inflaton oscillations, with the amplitude of the resonant modes growing exponentially inside a Bose-sphere of radius $k \lesssim k_\star \sim q_*^{1/4} \omega_\star$. This radiates GWs efficiently within a similar band of momenta. At the end of the process the resulting energy density spectrum of GWs exhibits a peak with amplitude and location (redshifted to frequencies today) given by (Figueroa and Torrenti 2017a)

$$f \approx 8 \cdot 10^9 \left( \frac{\omega_{\star}}{\rho_{\star}^{1/4}} \right) \epsilon_{\star}^{\frac{1}{4}} q_{\star}^{\frac{1}{4}+\eta} \text{ Hz}, \qquad (120a)$$

$$\Omega_{\text{GW}}^{0}(f) \approx \mathcal{O}(10^{-9}) \times \epsilon_{\star} C \frac{\omega_{\star}^{6}}{\rho_{\star} M_{\text{P}}^{2}} q_{\star}^{-\frac{1}{2}+\delta}, \qquad (120\text{b})$$

where $\rho_\star$ is the total energy density at $t = t_\star$, $\eta$ and $\delta$ are parameters that account for non-linearities of the system, and $C$ is a constant that characterises the strength of the resonance. The factor $\epsilon_\star \equiv (a_\star/a_{\text{RD}})^{1-3w}$ parametrises the expansion history between the end of inflation and the onset of RD, assuming an averaged EoS $w \ne 1/3$ during this period (if $w = 1/3$ then $\epsilon_* = 1$). The values for $C$, $\eta$, and $\delta$, can only be determined with classical lattice simulations in a model by model basis.

Leaving aside non-linear effects for simplicity, we observe that the frequency and amplitude of the peak, roughly scale as $f \sim q_*^{1/4}$ and $\Omega_{\text{GW},0} \sim q_*^{-1/2}$. This means that in order to shift the peak to observable frequencies, we need to decrease $q_\star$; but in doing so, we further decrease the amplitude of the signal. Using this linear approximation we see that decreasing the peak frequency say by an order of magnitude, implies the reduction of the amplitude of the background by two orders of magnitude. Furthermore, we see that these backgrounds can only be expected down to a minimum frequency, as efficient GW production requires $q_* > 1$ to sustain broad resonance during the field instability, and hence $q_* = 1$ marks a minimum frequency we can think of. Even though such scaling behaviours are modified when
considering corrections due to non-linearities, the logic persists (Figueroa and
Torrenti 2017a), and as a result these backgrounds are forced to be peaked at high-
frequencies way above those accessible by LISA. To be concrete, considering for
instance a quartic potential around the minimum, we obtain $f \simeq (10^7 - 10^8)$ Hz and
$\Omega_{GW,0} \simeq (10^{-13} - 10^{-11})$, when assuming the range $q_{\star} \in (1, 10^4)$.<fnref n="29" /> For a quadratic
potential around the minimum of the potential, we find instead (assuming $\epsilon_{\star} = 1$) a
frequency $f \simeq (10^8 - 10^9)$ Hz and amplitude $\Omega_{GW,0} \simeq (10^{-12} - 10^{-11})$, when
considering a resonance parameter $q_{\star} \in (10^4, 10^6)$. The conclusion is clear, these
backgrounds are completely out of reach of current and (so far) planned direct
detection GW experiments.

If the scalar field interactions induce a tachyonic effective mass in the preheat
species, GWs can also be strongly produced, e.g. during hybrid preheating (García-
Bellido and Figueroa 2007; García-Bellido et al. 2008; Dufaux et al. 2009). In this
case, contrary to parametric resonance, the present day frequency and amplitude of
the generated SGWB might be tuned to peak at small frequencies, while retaining a
large amplitude. This depends on the model parameters, but in general a strong fine-
tuning is required for this to happen. A similar circumstance arises in the case of the
GWs produced from oscillons produced during hilltop preheating (Antusch et al.
2017, 2018; Amin et al. 2018), which upon similar fine tuning of the potential
parameters, can also lead to observable GW backgrounds at low frequencies.

If the field species involved in preheating are of a different nature than just scalar
fields, then new channels of GW production open up. For example, GWs can be
produced during the out-of-equilibrium excitation of fermions after inflation, both for
spin-1/2 (Enqvist et al. 2012; Figueroa and Meriniemi 2013; Figueroa 2014) and
spin-3/2 (Benakli et al. 2019) fields. In this case a SGWB with large amplitude is also
forced to peak at high frequencies. GWs can also be generated when the preheat
fields are (Abelian and non-Abelian) gauge fields. For example gauge fields could be
coupled to a charged scalar field via standard gauge covariant derivatives like in
Dufaux et al. (2010), Figueroa et al. (2016), Tranberg et al. (2018), or to a pseudo-
scalar field through a derivative axial coupling as in Adshead et al. (2018), Adshead
et al. (2020a), Adshead et al. (2020b). Preheating can be remarkably efficient in the
second case (Adshead et al. 2015; Cuissa and Figueroa 2019), with the energy
density produced in GWs possibly reaching up to ~1% of the total energy in the
system for the strongest coupling strengths (Adshead et al. 2020a, b). The peak
frequency is however also very large in these scenarios as well.

In conclusion, preheating mechanisms are capable of creating very large SGWBs,
but these are naturally peaked at very high frequencies, which are typically beyond
the LISA window. Some particular models can sustain a sizeable amplitude at LISA
frequencies, but only at the expense of a strong fine-tuning of their parameters.

<footnote n="29">²⁹ The GW spectrum in the quartic potential case also features additional peaks, see Figueroa and Torrenti (2017a) for more details.</footnote>
## 8.7 Summary of distinctive gravitational-wave observables from inflation

In what follows we describe the relevant observables that can be used to break possible degeneracies between the SGWB from inflation, and the one due to other cosmological sources, as well as to astrophysical ones.

### 8.7.1 Chirality

Although parity violation has so far been observed only in weak interactions, it is important to investigate whether the same phenomenon may occur also in the very early universe, during inflation. In the primordial context, we shall refer to chirality whenever the two polarisations of GWs have different amplitudes, see Sect. 5.

For convenience, one may introduce the parameter $\chi \equiv |P_{\gamma}^{L} - P_{\gamma}^{R}| / \sum_{\lambda} P_{\lambda}^{\lambda}$, whose range extends up to $\chi = 1$. Several mechanisms for generating chiral GWs are found in the literature. The common trait to such realisations is the presence of a Chern–Simons-type interaction. In the case of the EM field strength such term reads

$$ \mathcal{L} \supset \phi F_{\mu\nu} \tilde{F}^{\mu\nu}, \qquad (121) $$

where $\tilde{F}^{\mu\nu}$ is the dual EM field. The reader will recognise the similarities between Eq. (101) and (121). The comparison underscores the fact that the field corresponding to $F$ does not need to be the Standard Model photon, nor an Abelian gauge field. In modified gravity theories, Chern–Simons gravity is obtained by promoting the three dimensional gravitational Chern–Simons term to 4D (Jackiw and Pi 2003; Alexander and Yunes 2009):

$$ \mathcal{L} \supset f(\phi) R_{\mu\nu\rho}^{\sigma} \tilde{R}_{\sigma}^{\mu\nu\rho}, \qquad (122) $$

where $\tilde{R}$ is the dual Riemann tensor. The matter field $\phi$ in the last two equations is typically (but not necessarily) identified as the inflaton in models of the early phase acceleration. The parity violation ensuing from Chern–Simons-type coupling in the early universe can be tested across a vast range of scales, from CMB (Lue et al. 1999; Thorne et al. 2018) to interferometers (Crowder et al. 2013; Smith and Caldwell 2017; Domcke et al. 2020b), including LISA. Both types of interactions in Eqs. (121) and (122) have been extensively studied (Anber and Sorbo 2010; Satoh 2010; Crisostomi et al. 2018; Bartolo and Orlando 2017; Bartolo et al. 2019d; Qiao et al. 2020; Bartolo et al. 2021; Bordin and Cabass 2020)), sometimes considered together (Lue et al. 1999; Mirzagholi et al. 2020a), and investigated as emerging from quantum gravity (in the sense of Hořava–Lifshitz theories, Horava 2009) (Takahashi and Soda 2009).

### 8.7.2 Frequency profile

The single-field slow-roll inflationary paradigm predicts a slightly red-tilted GW spectrum. The tilt of the tensor power spectrum is indeed proportional to the slow-roll parameter, $n_T \simeq -2\epsilon \ll 1$. Such a GW signal is well below the sensitivity of
LISA and could even elude detection from more sensitive instruments, including the proposed BBO. This is the key notion behind the claim that a detection of primordial GWs by LISA would be a smoking gun for an inflationary mechanism that goes beyond the single-field slow-roll scenario. With some interesting exceptions (Mylova et al. 2018) one may go further and state that a primordial signal at LISA scales (and sensitivity) is strongly suggestive of multi-field (or multi-clock) dynamics.

A typical case in point is that of GW sourced by extra fields during inflation. Whenever the sourced contribution goes beyond vacuum fluctuations, the richer multi-field dynamics can give rise to a non-trivial frequency dependence, from a blue GW spectrum to a “bump-like” structure and to oscillatory features mimicking the ones of the scalar power spectrum. For all these examples there exist explicit realisations, as detailed in Sects. 8.2 and 8.3.

In studying the SGWB frequency profile one ought to be aware of the bounds at CMB scales, those set at intermediate scales by PTAs and, perhaps most importantly for LISA, those set at the relatively close frequencies accessed by LIGO/Virgo.

### 8.7.3 Anisotropies

The angular resolution of the LISA detector might enable a detection of another peculiar feature useful in the SGWB characterisation process: the anisotropy (direction dependence) in the energy density. Such anisotropies contain information about the generation process of GWs and their propagation across cosmic inhomogeneities. Using a Boltzmann equation approach, the contribution coming from the generation mechanism retains a frequency dependence, which is peculiar of the SGWB, and, due to the non-thermal nature of the graviton distribution function at their decoupling time (Bartolo et al. 2019a, 2020b). The contribution arising from the propagation of GWs through large-scale cosmological scalar (and tensor) background perturbations, happens to be larger compared to the same effect for CMB photons. So, in a similar way to CMB photons, the SGWB from e.g. inflation is affected by the Sachs-Wolfe and the Integrated Sachs-Wolfe effects. The former is a gravitational redshift, due to the difference of the gravitational potential at the moment of the SGWB production and today. The latter is due to the variation of the gravitational potential along the line of sight from the SGWB production to its the detection, integrated in time (Bartolo et al. 2019a, 2020b).

It turns out that LISA is sensitive to these effects in the angular power spectrum if the SGWB isotropic energy density is $\Omega_{GW} \ge 10^{-12}$ (Bartolo et al. 2019a, 2020b; Valbusa Dall'Armi et al. 2021). For a representative inflationary model, namely axion inflation (Cook and Sorbo 2012; García-Bellido et al. 2016), it was shown that the predicted SGWB can be within the reach of LISA and exhibits anisotropies with a large frequency dependence, a possible distinctive target for the SGWB detection and characterisation in LISA.

In the remaining part of this subsection we will discuss GW anisotropies that result from squeezed primordial non-Gaussianity. This is yet another case in point for the use of GW anisotropies as a probe of the production mechanism for GWs. Let us briefly describe the genesis of these non-Gaussianity-sourced anisotropies. A
squeezed primordial bispectrum encodes a coupling between two short-wavelength modes and one long-wavelength mode. The effect of the very same coupling is also manifest at the level of the power spectrum of the short-wavelength modes, in the form of a modulation by the long-wavelength mode (Jeong and Kamionkowski 2012; Dai et al. 2013; Brahma et al. 2014; Dimastrogiovanni et al. 2014, 2016). The magnitude and form of this modulation depends on the specific type of interactions (as dictated by the inflationary Lagrangian) and on the nature (e.g. scalar vs. tensor) of the long-wavelength mode. Let us consider the case of a tensor bispectrum:

$$ \langle h_{\mathbf{k}_1} h_{\mathbf{k}_2} h_{\mathbf{q}} \rangle = (2\pi)^3 \delta^{(3)}(\mathbf{k}_1 + \mathbf{k}_2 + \mathbf{q}) B^{\text{ttt}}(\mathbf{k}_1, \mathbf{k}_2, \mathbf{q}). \quad (123) $$

In the squeezed limit, $q \ll k_1 \approx k_2$, and for models preserving statistical isotropy and parity, the long wavelength tensor mode imprints a quadrupolar modulation in the primordial power spectrum of the short-wavelength modes (Dimastrogiovanni et al. 2020):

$$ P_h^{\text{mod}}(\mathbf{k}, \mathbf{x}) = P_h(k) [1 + Q_{\ell m}(\mathbf{k}, \mathbf{x}) \hat{n}_{\ell} \hat{n}_{m}], \quad (124) $$

where $\mathbf{k} = k\hat{n}$ and we have defined

$$ Q_{\ell m}(\mathbf{k}, \mathbf{x}) = \int \frac{d^3 q}{(2\pi)^3} e^{i\mathbf{q} \cdot \mathbf{x}} \sum_{\lambda_3} h_{\ell m}^{\lambda_3}(\mathbf{q}) F_{\text{NL}}^{\text{ttt}}(\mathbf{k}, \mathbf{q}). \quad (125) $$

Here $F_{\text{NL}}^{\text{ttt}}(\mathbf{k}, \mathbf{q})$ is the amplitude of the tensor bispectrum in the squeezed limit, normalised by the product of the power spectra, $P_h(k) \cdot P_h(q)$. The anisotropy in Eq. (124) determines a contribution to the energy density contrast (see Eq. (35) for the definition of $\delta_{\text{GW}}$) given by

$$ \delta_{\text{GW}}(k, \hat{n}) = Q_{\ell m}(\mathbf{k}, \mathbf{d}) \hat{n}_{\ell} \hat{n}_{m}, \quad (126) $$

where $\mathbf{d} = -\hat{n}d$, $d \equiv \eta_0 - \eta_{\text{in}}$ being the separation in conformal time between horizon re-entry of the k-mode and the present time.

In a similar way, for inflationary models predicting a long-short mode coupling between scalars and tensors, a long-wavelength scalar fluctuation will modulate the power spectrum of GWs on small scales:

$$ P_h^{\text{mod}}(\mathbf{k}, \mathbf{x}) = P_h(k) \left[ 1 + \int \frac{d^3 q}{(2\pi)^3} e^{i\mathbf{q} \cdot \mathbf{x}} \zeta(\mathbf{q}) F_{\text{NL}}^{\text{stt}}(\mathbf{k}, \mathbf{q}) \right]. \quad (127) $$

This generates an anisotropic component for the GW energy density measured at present time, which leads to

$$ \delta_{\text{GW}}(k, \hat{n}) = \int \frac{d^3 q}{(2\pi)^3} e^{-i d \hat{n} \cdot \mathbf{q}} \zeta(\mathbf{q}) F_{\text{NL}}^{\text{stt}}(\mathbf{k}, \mathbf{q}), \quad (128) $$

with $F_{\text{NL}}^{\text{stt}}(\mathbf{k}, \mathbf{q})$ being the amplitude of the scalar-tensor-tensor bispectrum in the squeezed limit, normalised by $P_{\zeta}(q) \cdot P_h(k)$.
Primordial non-Gaussianity is one of the most informative probes of field
interactions during inflation. Already at the level of the bispectrum (the lowest-order
non-Gaussian correlator), a variety of constraints can be placed on inflationary
models based on its momentum dependence (a.k.a. shape) and overall amplitude. In
particular, bispectra that receive substantial contribution in the squeezed momentum
configuration are a smoking gun for inflationary dynamics that go beyond the single-
field slow-roll realisation. This is due to the fact that, whenever the inflationary
dynamics is guided by a single effective degree of freedom, the leading order
contribution to the squeezed bispectrum can be removed by a gauge transformation
(Maldacena 2003; Tanaka and Urakawa 2011; Pajer et al. 2013). Single-clock models
therefore predict a suppressed bispectrum. On the other hand, typical classes of
models that can lead to sizeable soft limits for cosmological correlators include
multiple fields, excited initial states, or an inflationary background that breaks space
diffeomorphism invariance.

While at large (e.g. CMB) scales primordial bispectra can be directly constrained,
this is not the case at small scales due to loss of coherence from propagation in the
perturbed universe (see e.g. discussion in Sect. 10.3). Anisotropies are therefore a
key observable for constraining tensor and mixed non-Gaussianity at interferometer
scales.

As an example, in the left panel of Fig. 23 we display the auto-correlation
$\langle \delta_{\text{GW},\ell_1 m_1} \delta_{\text{GW},\ell_2 m_2} \rangle = \delta_{\ell_1 \ell_2} \delta_{m_1 m_2} C_{\ell_1}^{\text{GW}}$ arising from anisotropies of the kind described
in Eq. (128) (red line). For the sake of comparison, in the same plot we show the
auto-correlation for anisotropies imprinted as gravitons travel in the perturbed
background after horizon re-entry (blue line) (Alba and Maldacena 2016; Contaldi
2017; Bartolo et al. 2019a, 2020b).

The recent works (Contaldi et al. 2020; Bartolo et al. 2022) study the capability of
LISA to access the SGWB anisotropies. See Sect. 12 for more details.

8.7.4 Cross-correlations

Anisotropies such as those just discussed in Sect. 8.7.3 are produced by the effect of
long-wavelength (scalar or tensor) perturbations and are therefore correlated with
anisotropies in the CMB (Adshead et al. 2021; Malhotra et al. 2021). These cross-
correlations provide an additional observable for constraining primordial non-
Gaussianity (see right panel of Fig. 23 for an application). As an example, we report
in Fig. 24 the error in determining $F_{\text{NL}}^{\text{stt}}$ using a measurement of these cross-
correlations. The calculation of the error follows that of Malhotra et al. (2021) and it
is adapted here, under that assumption that Taiji and LISA should happen to fly
together (Ruan et al. 2020b; Orlando et al. 2021). We see that for the combination
$\{r = 0.05, n_T = 0.30\}$ one would be able to detect $F_{\text{NL}}^{\text{stt}}$ of order $10^4$ or larger.

For any given inflationary model generating non-trivial scalar-tensor-tensor or
tensor-tensor-tensor bispectra in the squeezed limit, the overall level of the GW
anisotropies is determined by the specific form of the angular dependence of the
bispectra and by the magnitude of the $F_{\text{NL}}$ parameters. Typically, there exist
minimum values of $F_{\text{NL}}$ below which the GW anisotropy map is dominated by noise.
<figure><img src="image_27.png" /><figcaption>Fig. 23 Left panel: Auto-correlation of GW anisotropies induced by a scalar-tensor-tensor squeezed bispectrum with $F_{\text{NL}}^{\text{stt}} \sim \mathcal{O}(10^3)$ (red line), and of GW anisotropies arising from propagation through the perturbed background (blue line). Right panel: Cross-correlation of GW anisotropies and CMB temperature anisotropies in the same two cases</figcaption></figure>

Cross-correlations can be particularly helpful in these cases as they have the potential to constrain smaller levels of non-Gaussianity (Malhotra et al. 2021).

# 9 Tests of non-standard pre-Big-Bang nucleosynthesis cosmology via the SGWB

Section coordinator: G. Calcagni. Contributors: I. Ben-Dayan, G. Calcagni, C-F. Chang, Y. Cui, D.G. Figueroa, S. Kuroyanagi, M. Lewicki, A. Mazumdar, G. Servant, P. Simakachorn.

## 9.1 Introduction

In the standard model of cosmology, the universe begins with a inflationary epoch that first induces an exponential growth of the universe and then reheats it with a very hot plasma. In the standard picture, it is the energy density of such a plasma that dominates the expansion of the universe until the matter and DE domination final stages. It is however worth remembering that this picture is based on several

<figure><img src="image_28.png" /><figcaption>Fig. 24 The expected $1\sigma$ error in determining $F_{\text{NL}}^{\text{stt}}$ with cross-correlations of CMB temperature anisotropies and GW anisotropies, as a function of $\ell_{\text{max}}$</figcaption></figure>
theoretical assumptions and has only been tested up to temperatures $T \approx 5$ MeV thanks to data from CMB (Aghanim et al. 2020) and successful predictions of BBN (Cyburt et al. 2016). The period between that temperature and inflation still holds many unknowns (Allahverdi et al. 2020). In this section, we will discuss how detection of a SGWB by LISA would allow us to probe the evolution of the universe during that time. We will begin with the possibility of probing the expansion rate including the number of relativistic degrees of freedom thanks to the SGWB sources discussed in the previous sections and move on to possible departures from GR.

According to the standard inflationary paradigm in GR, scalar and tensor perturbations are generated during an early era of accelerated expansion induced by a slow-rolling phase. In the vanilla inflationary scenario, the primordial tensor spectrum is red-tilted at CMB scales, i.e. its amplitude decreases with frequency. If inflation is followed by the long-lasting RD period, then also the produced SGWB, late-time remnant of the primordial tensor spectrum, turns out to be red-tilted and its amplitude at the mHz—Hz frequencies remains too small to be detected by any present or planned near-future GW detector. However, both deviations from the vanilla inflationary scenario (see Sect. 8) and the usual RD epoch of the standard model of cosmology can enhance the SGWB to amplitudes large enough to reach the LISA sensitivity. In the latter case, for instance, if the expansion rate of the universe between the end of inflation and the onset of RD, is dictated by an effective EoS parameter $w \neq 1/3$, the inflationary SGWB is naturally tilted at high frequencies, for the modes that re-entered the Hubble radius during such non-standard epoch. If the EoS is stiff, i.e. $w > 1/3$, then the tilt is blue, thus enhancing the power of the GW spectrum at large frequencies, and hence making this branch of the spectrum potentially observable with LISA. If, on the contrary, the EoS is $w < 1/3$, then the tilt at high frequencies is red, and hence the spectrum goes below the standard plateau amplitude due to modes crossing during RD.

Any modification of the expansion rate occurring after a FOPT also leaves its imprint on the SGWB produced by that phase transition. The modified redshift results in an overall shift of the entire signal but, more interestingly, characteristic features might also be produced in the signal. This is possible for modes with lengths larger than the horizon size at the time of the transition. If these modes enter the horizon during a time of modified expansion, the typical $f^3$ slope at low frequencies will be modified allowing us to decode the EoS at that time. However, as this feature is always present significantly below the main amplitude peak of the FOPT SGWB (see Sect. 6), measuring the expansion using this method will always be much more challenging than the initial detection of the spectrum.

Cosmic strings are another possible GW source which would also allow us to probe the expansion history. As a network of strings would continuously produce GWs throughout most of cosmological history (see Sect. 7), also the expansion history is encoded in the resulting spectral shape. This, in fact, includes not only drastic modifications such as periods of non-standard expansion, but also subtler modifications including the number of degrees of freedom contributing to the radiation dominating at early times. As we will see, an accurate observation of a
SGWB from cosmic strings will always give us some information on the expansion rate of the very early universe.

Finally, we will also explore the possibility to test quantum gravity cosmological models via SGWB observations at LISA. After excluding a plethora of models whose SGWB does not reach the instrument's sensitivity curve, we single out two pre-Big-Bang scenarios as candidates that could generate a signal relevant for LISA science.

## 9.2 Non-standard expansion histories

### 9.2.1 Inflation

Equation (99) provides the amplitude of the vacuum tensor modes produced at super-horizon scales during inflation. During the expansion history following inflation, the tensor modes re-enter successively back inside the Hubble radius, turning into a proper classical (yet stochastic) background of GWs. Once the modes become sub-Hubble, they start oscillating with a decaying amplitude $h_{ij} \propto 1/a$ (this is independent of the expansion rate (Caprini and Figueroa 2018)), propagating as relativistic degrees of freedom. Since modes with different wavelengths re-enter the Hubble radius at different moments of cosmic evolution, different modes may propagate through different periods of expansion. As a result, modes with very different wavelengths may sustain very different amplitudes with respect to each other, depending on the rate of expansion of the universe at their time of Hubble-crossing. To characterise the spectrum of the SGWB today, we can write the energy density spectrum normalised to the critical density $\rho_{\text{crit}} = 3M_{\text{Pl}}^2H^2$, as (Caprini and Figueroa 2018)

$$ \Omega_{\text{GW}}(t, k) \equiv \frac{1}{\rho_{\text{crit}}} \frac{d\rho_{\text{GW}}(t, k)}{d \ln k} = \frac{k^2}{12a^2(t)H^2(t)} P_h(t, k), \quad (129) $$

where the power spectrum $P_h(t, k)$ at a generic time $t$ is related to the power $P_h(k)$ generated super-horizon during inflation by a transfer functions,

$$ P_h(t, k) \equiv T_h(t, k)P_h(k), \quad T_h(t, k) \approx \frac{1}{2} \left[ \frac{a_k}{a(t)} \right]^2. \quad (130) $$

The transfer function $T_h(t, k)$ depends on the expansion history between the moment $t = t_k$ of horizon re-entry of a given mode $k$, and a later moment $t > t_k$. Here $t_k$ is implicitly defined from the condition $a_k H_k \equiv k$, with $a_k \equiv a(t_k)$, and $H_k \equiv H(t_k)$. The factor $1/2$ is due to the time-average of the rapidly-oscillating sub-horizon wave.

The transfer function today can be evaluated as

$$ T_h(k) \approx \frac{1}{2} \left( \frac{a_k}{a_0} \right)^2 \approx \frac{1}{2} G_k \Omega_{\text{rad}}^{(0)} \left( \frac{a_0 H_0}{a_k H_k} \right)^2, \quad G_k \equiv \left( \frac{g_{*,k}}{g_{*,0}} \right) \left( \frac{g_{s,0}}{g_{s,k}} \right)^{4/3}, \quad (131) $$

where $g_s$ and $g_*$ are the relativistic number density of species contributing to the total entropy and energy densities, respectively. Using $\Omega_{\text{rad}}^{(0)} \approx 9 \cdot 10^{-5}$, $g_{s,0} \approx 3.91$,