$g_{*,0} = 2$, and $g_{s,k} \simeq g_{*,k} \simeq 106.75$, so that $G_k \simeq 0.65$, we obtain that tensor modes crossing during the RD epoch lead to an energy density spectrum with amplitude

$$ \left. \Omega_{\text{GW}}^{(0)} \right|_{\text{inf}} \approx G_k \frac{\Omega_{\text{rad}}^{(0)}}{12\pi^2} \left(\frac{H_{\text{inf}}}{m_{\text{Pl}}}\right)^2 \left(\frac{H_{\text{inf}}}{H_{\text{max}}}\right)^2 \left(\frac{f}{f_{\text{Pl}}}\right)^{n_T} \approx 5 \times 10^{-16} \left(\frac{H_{\text{inf}}}{H_{\text{max}}}\right)^2 \left(\frac{f}{f_{\text{Pl}}}\right)^{n_T}, \qquad (132) $$

with $n_T$ given by the original inflationary tilt. In other words, the RD energy spectrum retains the (quasi-) scale invariant spectral shape of the original inflationary tensors. We note that in evaluating Eq. (132) we have considered for concreteness that $g_{s,k}, G_{*,k}$ equal the SM degrees of freedom before electroweak symmetry breaking, and hence they are independent of $k$. In reality the number of SM relativistic degrees of freedom change with the scale, but for simplicity in Eq. (132) we considered an identical suppression for all the modes as $G_k \simeq 0.65$, so that we can provide a single number for the amplitude of this plateau.

If between the end of inflation and the onset of RD there is a period of expansion characterised by an EoS different than that of radiation $w \neq 1/3$, the inflationary GW energy density spectrum develops a tilt within the range of scales corresponding to the modes crossing the Hubble radius during such period. As a result, the (quasi-) scale invariance of the original tensor spectrum is lost. This feature is actually quite interesting from an observational point of view, as we might detect or constrain in this way the post-inflationary expansion history, and hence the properties of the fields driving the expansion between inflation and RD (Giovannini 1998, 1999; Boyle and Steinhardt 2008; Watanabe and Komatsu 2006; Boyle and Buonanno 2008; Kuroyanagi et al. 2009, 2011, 2015; Figueroa and Tanin 2019ab; Gouttenoire et al. 2021a).

In order to understand better the above discussion, let us consider that between the end of inflation and the onset of RD, the expansion rate is dictated by an effective EoS parameter $w \neq 1/3$. In scenarios where the inflaton potential is a monomial, $V(\phi) \propto \phi^p$, the inflaton oscillates around the minimum of such a potential after inflation, so that an effective EoS, averaged over inflaton oscillations, emerges as $w \simeq (p-2)/(p+2)$ (Turner 1983; Lozanov and Amin 2017). For $p<4$ the EoS is in the range $0<w<1/3$. For $p=2$, the oscillations of a massive free field lead to an energy density redshifting on average (over one oscillation cycle) as $\langle \rho_{\phi} \rangle_{\text{osc}} \propto 1/a^3$, analogous to that of non-relativistic particle species. Thus, an EoS of $w \simeq 0$ emerges in that case. Nothing prevents however the possibility of a stiff dominated stage with EoS $w = w_s$ for $1/3 < w_s < 1$. Such a case can be actually achieved quite naturally if the kinetic energy of the inflaton dominates after inflation, either through inflaton oscillations (Turner 1983) under a steep potential (e.g. $V(\phi) \propto \phi^p$ with $p > 4$), or simply by an abrupt drop of the inflationary potential at at moment that triggers itself the end of inflation.

Propagating the inflationary tensor modes through the epoch starting immediately after the end of inflation, leads to a GW energy density spectrum today, expressed as a function of present-day frequencies $f = k/(2\pi a_0)$, like (Figueroa and Tanin 2019a)
$$ \Omega_{\text{GW}}^{(0)}(f) = \Omega_{\text{GW}}^{(0)}|_{\text{inf}} \times W(f/f_{\text{RD}}) \times C_w \left(\frac{f}{f_{\text{RD}}}\right)^{2\left(\frac{3w-1}{3w+1}\right)} \approx \Omega_{\text{GW}}^{(0)}|_{\text{inf}} \\ \times \begin{cases} 1, & f \ll f_{\text{RD}} \\ C_w \left(\frac{f}{f_{\text{RD}}}\right)^{2\left(\frac{3w-1}{3w+1}\right)}, & f \gg f_{\text{RD}} \end{cases} \qquad (133) $$

with $C_w$ a numerical prefactor in the range $1<C_w<2^{5/2}/\pi \approx 1.80$ for $1/3<w<1$, and where $f_{\text{RD}} \equiv k_{\text{RD}}/(2\pi a_0)$ is the frequency corresponding to the horizon scale at the onset of RD, $k_{\text{RD}} = a_{\text{RD}}H_{\text{RD}}$, and $W(x)$ is a window function connecting the two sides of the spectrum: the high-frequency branch on the one hand, corresponding to modes crossing during the non-standard period ($f \gg f_{\text{RD}}$), and the low frequency branch corresponding to modes crossing during RD ($f \ll f_{\text{RD}}$). The tilt in Eq. (133) indicates that models with a stiff epoch between the end of inflation and the onset of RD are actually very appealing observationally, as in these scenarios a blue tilt is developed in the GW energy density spectrum at large frequencies, with $0<n_T<1$ (Giovannini 1998, 1999; Riazuelo and Uzan 2000; Sahni et al. 2002; Tashiro et al. 2004; Boyle and Buonanno 2008; Giovannini 2008, 2009; Caprini and Figueroa 2018; Bernal and Hajkarim 2019; Figueroa and Tanin 2019a). This potentially opens up the possibility of detection by upcoming direct detection GW experiments, including LISA. The presence of a stiff period is actually well-motivated theoretically in scenarios like Quintessential inflation (Peebles and Vilenkin 1999; Peloso and Rosati 1999; Huey and Lidsey 2001; Majumdar 2001; Dimopoulos and Valle 2002; Wetterich 2014, 2015; Hossain et al. 2014; Rubio and Wetterich 2017), gravitational reheating (Ford 1987; Spokoiny 1993) (with the caveats explained in Figueroa and Tanin (2019b)), the Higgs-reheating scenario (Figueroa and Byrnes 2017) and generalisations (Dimopoulos and Markkanen 2018; Opferkuch et al. 2019). The GW spectrum in these scenarios is controlled by $w, f_{\text{RD}},$ and $H_{\text{inf}}$, and the parameter space compatible with a detection by various experiments has been recently analysed in Gouttenoire et al. (2021a), Bernal and Hajkarim (2019), Figueroa and Tanin (2019a). In Fig. 25 we can see the form of the spectrum in the left panel (including scale-dependent changes in $G_k$ due to changes in the number of relativistic species). In the right panel of the same figure we plot the parameter space that LISA can probe. The tilt in Eq. (133) indicates that models with a stiff epoch between the end of inflation and the onset of RD are actually very appealing observationally, as in these scenarios a blue tilt is developed in the GW energy density spectrum at large frequencies, with $0<n_T<1$ (Giovannini 1998, 1999; Riazuelo and Uzan 2000; Sahni et al. 2002; Tashiro et al. 2004; Boyle and Buonanno 2008; Giovannini 2008, 2009; Caprini and Figueroa 2018; Bernal and Hajkarim 2019; Figueroa and Tanin 2019a). This potentially opens up the possibility of detection by upcoming direct detection GW experiments, including LISA. The presence of a stiff period is actually well-motivated theoretically in scenarios like Quintessential inflation (Peebles and Vilenkin 1999; Peloso and Rosati 1999; Huey and Lidsey 2001; Majumdar 2001; Dimopoulos and Valle 2002; Wetterich 2014, 2015; Hossain et al. 2014; Rubio and Wetterich 2017), gravitational reheating (Ford 1987; Spokoiny 1993) (with the caveats explained in Figueroa and Tanin
<figure><img src="image_29.png" /><figcaption>Fig. 25 Left panel: GW energy spectra including changes in the number of relativistic degrees of freedom, depending on whether the transition from stiff domination to RD takes place (red dashed line) or after (blue solid line) than the QCD phase transition. For comparison we show the corresponding spectra (grey solid line) without correcting for changes in the number of degrees of freedom. Image reproduced with permission from Figueroa and Tanin (2019a), copyright by IOP/SISSA. Right panel: Parameter space region that LISA can probe after removing the region (above the dashed horizontal line) incompatible with the BBN bound on GW backgrounds. Here $w_s$ refers to the EoS for a stiff background</figcaption></figure>

(2019b)), the Higgs-reheating scenario (Figueroa and Byrnes 2017) and generalisations (Dimopoulos and Markkanen 2018; Opferkuch et al. 2019). The GW spectrum in these scenarios is controlled by $w$, $f_{\text{RD}}$, and $H_{\text{inf}}$, and the parameter space compatible with a detection by various experiments has been recently analysed in Bernal and Hajkarim (2019), Figueroa and Tanin (2019a), Gouttenoire et al. (2021a). In Fig. 25 we can see the form of the spectrum in the left panel (including scale-dependent changes in $G_k$ due to changes in the number of relativistic species). In the right panel of the same figure we plot the parameter space that LISA can probe.

Consistency with upper bounds on SGWBs, like the BBN and CMB constraints, rules out a significant fraction of the observable parameter space. (See discussion at the end of Sect. 5.1.) This renders for instance this signal as unobservable by Advanced LIGO, independently of the parameter space {$w, f_{\text{RD}}, H_{\text{inf}}$} (Figueroa and Tanin 2019a). The GW background remains detectable in LISA, but only in a small island of parameter space, see right panel of Fig. 25 corresponding contrived scenarios having a low EoS $0.46 \lesssim w \lesssim 0.56$, a high inflationary scale $H_{\text{inf}} \gtrsim 10^{13}$ GeV at the same, and a low transition frequency $10^{-11} \text{ Hz} \lesssim f_{\text{RD}} \lesssim 3.6 \times 10^{-9} \text{ Hz}$ (or equivalently a low reheating temperature $1 \text{ MeV} \lesssim T_{\text{RD}} \lesssim 150 \text{ MeV}$). There exists other bounds beyond those from BBN and CMB on the duration of the stiff era (Gouttenoire et al. 2021a).

A short kination era ($w = 1$) happening much after reheating and inside the radiation era can be realized naturally by rotating-axion dynamics in the early universe (Gouttenoire et al. 2021b; Co et al. 2022; Gouttenoire et al. 2021a). Such kination typically lasts for 2–10 e-foldings and occurs after a matter era such that the inflationary SGWB gets boosted into a peak shape that is highly observable by future planned GW observatories while evading the BBN and CMB constraints in Fig. 25. LISA can probe the striking signature for kination eras happening at 100 GeV–
100 TeV energy scales (Gouttenoire et al. 2021a, b). Interestingly, such GW peaked signature would be a unique probe of axion physics.<fnref n="30" />

If there is an intermediate phase with EoS in the range $0 \le w < 1/3$, the transfer function of the inflationary GW spectrum in Eq. (133) develops instead a red-tilted high-frequency branch, corresponding to the modes propagating through that phase. If we could detect the (quasi-)scale invariant GW plateau part of the spectrum, the transition to the high-frequency branch due to the non-standard era, could be used to determine the reheating temperature $T_{\text{RD}}$ of the universe. This is because the end of the non-standard era also corresponds to the onset of RD, and hence the pivot frequency $f_{\text{RD}}$ separating the two branches in the energy density spectrum informs directly about the energy scale at the onset of RD (Nakayama et al. 2008; Kuroyanagi et al. 2009, 2011; D’Eramo and Schmitz 2019). Unfortunately, given the suppression of the amplitude of the plateau, c.f. Eq. (132), the determination of $T_{\text{RD}}$ by this method cannot be probed with LISA, and only a futuristic experiment like BBO or DECIGO (Crowder and Cornish 2005; Seto et al. 2001; Sato et al. 2017) might have a chance.

Another interesting opportunity is to look for a signal of particles beyond the SM. When particles becomes non-relativistic as the temperature drops below the mass and decouple from thermal equilibrium, they no longer contribute to the radiation energy density. This can be seen as a change in the values of $g_{*,k}$ and $g_{s,k}$, resulting in a temporary speed-up of the Hubble expansion rate. SM particles induce changes in $g_{*,k}$ and $g_{s,k}$ in the frequency range of $\sim 10^{-12}$ to $10^{-5}$ Hz (Watanabe and Komatsu 2006), which is beyond the range of LISA sensitivity curve. However, any change due to BSM particles with mass near 100 TeV leave an imprint on the SGWB spectrum at the LISA sensitivity. This effect is not detectable in the case of standard slow-roll inflation, since the amplitude of the SGWB is far below the LISA sensitivity as in Eq. (132), while a strongly blue-tilted primordial spectrum would provide an opportunity to test BSM physics by LISA (Caldwell et al. 2019).

### 9.2.2 Phase transitions

We discussed the production of GWs in FOPTs in detail in Sect. 6. We will now discuss the impact a period of non-standard evolution would have on these backgrounds. The first most obvious change will come in through modified red-shifting of the spectra from the time of their production until today (Allahverdi et al. 2020; Gouttenoire et al. 2021a). Simply parametrising the non standard evolution through a barotropic parameter $w$ we find for the abundance

<footnote n="30">³⁰ Such a peaked GW signature from a kination era induced by a spinning complex scalar field in the context of self-interacting DM was found and studied in Li et al. (2014), Rindler-Daller and Shapiro (2014), Li et al. (2017), Li and Shapiro (2021) where the initial rotation is postulated and the preceding matter era is a consequence of inflaton oscillations rather than axion dynamics.</footnote>
$$ \Omega_{\text{GW},0} = \left(\frac{a_*}{a_0}\right)^4 \left(\frac{H_*}{H_0}\right)^2 \Omega_{\text{GW},*} = 1.67 \times 10^{-5} h^{-2} \left[ \frac{100}{g_s(T_{\text{RD}})} \right]^{\frac{4}{3}} \left[ \frac{g_*(T_{\text{RD}})}{100} \right] \left(\frac{H_*}{H_{\text{RD}}}\right)^{\frac{2(3w-1)}{3w+3}} \Omega_{\text{GW},*}, \quad (134) $$

and for the frequency

$$ \begin{aligned} f_0 &= \frac{a_*}{a_0} f_* = \frac{a_{\text{RD}} H_{\text{RD}}}{a_0} \frac{a_* H_*}{a_{\text{RD}} H_{\text{RD}}} \frac{f_*}{H_*} \\ &= 1.65 \times 10^{-5} \, \text{Hz} \left( \frac{T_{\text{RD}}}{100 \, \text{GeV}} \right) \left( \frac{100}{g_s(T_{\text{RD}})} \right)^{\frac{1}{3}} \left( \frac{g_*(T_{\text{RD}})}{100} \right)^{\frac{1}{2}} \left( \frac{f_*}{H_*} \right) \left( \frac{H_*}{H_{\text{RD}}} \right)^{\frac{3w+1}{3w+3}}, \end{aligned} $$

where the asterisk (*) denotes quantities at the time of the transition, while “RD” refers to the moment when the universe becomes radiation dominated. In a decelerating universe $H_*/H_{\text{RD}} > 1$ and the direction of these modifications depends on the barotropic parameter. In radiation domination, $w = 1/3$, the amplitude does not change as GWs redshift at the same rate as the background. For an expansion dominated by energy density redshifting slower than radiation $w < 1/3$ the amplitude will keep decreasing while in the opposite case of an extra component redshifting faster $w > 1/3$ the abundance would increase looking at only the effect of redshift.

Another crucial modification comes from the conditions at the time of the transition. While typically the GW production is not crucially modified (Barenboim and Park 2016; Guo et al. 2021), the transition is linked to the radiation component and if that is subdominant the GW abundance also needs to be corrected. In terms of the strength of the transition this leads to a simple modification by

$$ \alpha \propto \rho_R^{-1} \rightarrow (\rho_{\text{tot}} - \rho_V)^{-1}, \quad (135) $$

where $\rho_R$ is radiation density while $\rho_V$ is the difference between the initial and final state energy densities. Finally $\rho_{\text{tot}}$ is the total energy density which in domination of any additional component is of course the dominant contribution.

This modification makes observation of scenarios where the transition occurs during the domination of an extra component rather unlikely. Even in the optimistic case of kination (with $w = 1$) where the redshift of the signal would increase the amplitude from Eq. (135) we get $(\rho_V/\rho_{\text{tot}})^2 \approx (H_*/H_{\text{RD}})^{-4\frac{3w-1}{3w+3}}$ which results in an overall suppression. Due to this, to get any hope of observation we would need the modified expansion to end almost immediately after the transition (Barenboim and Park 2016) despite the fact that the two are governed by completely separate mechanisms. Thus, the most promising models from the observational point of view are the ones in which the transition occurs in standard circumstance, but it is strong enough to itself modify the expansion history, as in this case Eq. (135) does not apply. One example here would be a very strong transition in which the field dominates the total energy density and after the transition as the field oscillates around its minimum it is causing a matter dominated period which will last until the field decays (Ellis et al. 2020c).
All the modifications of the spectrum we discussed up to now, essentially come down to shifts of the entire spectrum that could also be mimicked by simple changes in the parameters of the transition. We now switch to a modification changing the spectral shape which could provide a smoking gun signal for a modified expansion period.

This modification has to do with the fact that at scales larger than the horizon size at the time of the transition the features of our source become irrelevant. The source at this large scale is effectively just white noise and GWs entering the horizon at later time during RD always predict a spectrum proportional to $f^3$ (Caprini et al. 2009a; Ünal 2019; Cai et al. 2020) at frequencies corresponding to superhorizon scales $f<f_{H_*} = a_*H(a_*)/2\pi$. This behaviour, however, depends on the expansion rate at the time when GWs enter the horizon (Hook et al. 2021). In our relevant example of matter domination, it would create an $f^1$ plateau in the spectrum for $f_{H_{RD}}<f<f_{H_*}$. The main issue with observation of this feature comes from the fact that the peak of the signal has to do with the characteristic scale of the transition, which is typically much smaller than horizon size such that $f_*/f_{H_*} \propto \beta/H$, which is typically much bigger than one. As a result, this feature will typically appear significantly below the peak and, given that the abundance decreases quickly off peak, it will also have a much smaller abundance. Concerning LISA scales, in Fig. 26, we show several examples of this modification on bubble collision spectra calculated according to Lewicki and Vaskonen (2021) using parameter examples not difficult to realise in classically conformal models (Ellis et al. 2020c). From top to bottom the lines show the decay width of the scalar field $\Gamma_\phi/H = 1$, $\Gamma_\phi/H = 10^{-2}$ and $\Gamma_\phi/H = 10^{-4}$ which results in longer matter-dominated periods and reheating temperatures of $T_{RD} = 10^4$ GeV, $T_{RD} = 10^3$ GeV and $T_{RD} = 10^2$ GeV, respectively. We also stress that departure from standard cosmology that we have just explained, is not the only option leading to relevant chances of the SGWB frequency profile. For instance, the occurrence of an intermediate matter and kination era inside the radiation epoch as

<figure><img src="image_30.png" /><figcaption>**Fig. 26** The lines show the GW spectrum produced in FOPT by bubble collisions for the indicated transition parameters. From top to bottom, the decay width of the scalar field decreases leading to a lengthening period of effectively matter-dominated expansion as the field oscillates around its minimum before finally decaying. The resulting reheating temperatures read $T_{RD} = 10^4$ GeV, $T_{RD} = 10^3$ GeV and $T_{RD} = 10^2$ GeV, respectively</figcaption></figure>
motivated by e.g. early axion dynamics also impacts the low-frequency slope of the SGWB signal (Gouttenoire et al. 2021b).

### 9.2.3 Cosmic strings

We now turn to the case of GW spectra produced by cosmic strings that we first discussed in Sect. 7. We focus on analytical modelling described in Sect. 7.2.1 with numerical factors set to agree with Model 2 described in Sect. 7.2.2.

These spectra are very convenient laboratories to study the expansion history of the early universe, since in RD they simply produce a flat plateau. Any features beyond that can be linked to modifications of the expansion rate from minor modifications caused by variations in the number of degrees of freedom to simply a different power-law caused by domination of energy density red-shifting differently than radiation (Cui et al. 2018, 2019; Gouttenoire et al. 2020a, b). The feature will appear in the spectrum above at the characteristic frequency $f_{RD}$ corresponding to the temperature $T_{RD}$ at which the expansion begins to follow the standard radiation-dominated picture (with SM number of degrees of freedom) (Cui et al. 2019):

$$f_{RD} = (8.67 \times 10^{-3} \text{ Hz}) \left(\frac{T_{RD}}{\text{GeV}}\right) \left(\frac{10^{-11}}{G\mu}\right)^{1/2} \left[\frac{g_{*}^{\text{SM}}(T_{RD})}{g_{*}^{\text{SM}}(T_{0})}\right]^{\frac{8}{6}} \left[\frac{g_{s}^{\text{SM}}(T_{0})}{g_{s}^{\text{SM}}(T_{RD})}\right]^{\frac{7}{6}}, \quad (136)$$

where $G\mu$ is the tension characterising the string network. We stress that, for local strings, the relevant energy scale is not the size of the horizon, in contrast with phase transitions for instance. This is because loops do not suddenly decay after production, in contrast with other cosmological sources of GW. So, for LISA, the cosmic-strings GW spectrum is sensitive to MeV scale, rather than the 100 GeV–TeV scale (Gouttenoire et al. 2020a). In Fig. 27, we show the part of the ($G\mu, T_{RD}$) parameter space in which LISA can probe the spectra at $f_{RD}$ testing the standard cosmological expansion rate.

For global strings, the turning-point frequency turns out to be (Gouttenoire et al. 2020a)

$$f_{RD} = \left(\frac{T_{RD}}{\text{GeV}}\right) \left(\frac{0.1}{\alpha}\right) \left[\frac{g_{*}^{\text{SM}}(T_{RD})}{g_{*}^{\text{SM}}(T_{0})}\right]^{\frac{1}{4}} \times \begin{cases} 8.9 \times 10^{-7} \text{ Hz for } 10\%, \\ 7.0 \times 10^{-8} \text{ Hz for } 1\%, \end{cases} \quad (137)$$

where $\alpha$ is the loop size with respect to the Hubble horizon and the percentage

<figure><img src="image_31.png" /><figcaption>**Fig. 27** Range of temperature $T_{RD}$ to which LISA can probe the cosmological expansion rate using a spectrum from a cosmic string network with the indicated string tension $G\mu$. The grey region indicates temperatures where modifications of the expansion rate would already be in tension with BBN</figcaption></figure>
indicates the relative deviation $|(\Omega_{\text{st}} - \Omega_{\text{nonst}})/\Omega_{\text{st}}|$ of the non-standard spectrum with respect to the standard one. Thus, in contrast to local strings, the turning-point frequency is independent of the string tension (cf. Eqs. (136) and (137)).

Coming back to the local-string case, assuming that at earlier times the expansion is governed by a simple EoS with $\rho \propto a^{3(1+w)}$, we can approximate the slope of the spectrum at $f > f_{\text{RD}}$ as

$$ \Omega_{\text{GW}}(f > f_{\text{RD}}) \propto \begin{cases} f^{2\frac{3w-1}{3w+1}} & \text{for } w \ge 5/21, \\ f^{-1/3} & \text{for } w < 5/21. \end{cases} \qquad (138) $$

We show examples of modified spectra for spectrum from a cosmic string network $G\mu = 10^{-10}$ including early matter domination (with $w = 0$) and kination (with $w = 1$) lasting up until $T_{\text{RD}} = 5$ MeV and $T_{\text{RD}} = 5$ GeV in Fig. 28. A change in the number of degrees of freedom would instead create a smooth step in the spectrum with the total change in abundance given by (Cui et al. 2019),

$$ \Omega_{\text{GW}}(f \gg f_{\text{RD}}) \approx \Omega_{\text{GW}}^{\text{SM}}(f) \left( \frac{g_{*}^{\text{SM}}}{g_{*}^{\text{SM}} + \Delta g_{*}} \right)^{\frac{1}{3}}, \qquad (139) $$

where the SM index denotes quantities computed assuming the number of degrees of freedom as in the SM while $\Delta g_*$ is the number of new degrees of freedom decoupling at $T_{\text{RD}}$. According to this simple formula probing the abundance of the plateau with $\mathcal{O}(1\%)$ accuracy would allow us to discover inclusion of even several new degrees of freedom with larger amounts requiring less accuracy. The range of temperatures we can survey is again given by Fig. 27 although ascertaining the exact accuracy concerning our possible estimation of the number of degrees of freedom in this range would require further scrutiny. It was also shown that measurements of the turning-point frequency in the SGWB from cosmic strings due to a temporary matter era induced by a massive unstable particle can enable to probe unstable particles with

<figure><img src="image_32.png" /><figcaption>Fig. 28 GW spectrum produced by a cosmic string network with $G\mu = 10^{-10}$ together with spectra produced if that network evolved in an early matter domination period (with $w = 0$) or kination (with $w = 1$) lasting up until $T_{\text{RD}} = 5$ MeV and $T_{\text{RD}} = 5$ GeV. We note that the blue tilted branches in the figure are only shown for representative purposes, as they violate BBN and CMB constraints [recall discussion at the end of Sect. 5.1], so for a realistic effect compatible with those bounds, a much smaller stiff EoS is required, somehow larger than (but very close to) $w \approx 1/3$</figcaption></figure>
lifetime in the range $10^{-8}$ to 0.1 s, extending significantly the well-known BBN bound (Gouttenoire et al. 2020b).

In addition to the potential of probing a new EoS and new relativistic degrees of freedom with a SGWB from cosmic strings, (pre-)inflationary dynamics may leave distinctive imprints in GW signals. Although a cosmic string network formed during or before inflation is diluted by inflation, it may come back into horizon at later times and lead to detectable GW signals. In this case, the SGWB typically gets suppressed especially at higher frequencies (Guedes et al. 2018). Nevertheless, well-motivated BSM scenario of an intermediate short inflationary period inside the radiation era can lead to detectable effects. With LISA, intermediate inflation scales in the range $0.1-10^{12}$ GeV can be probed (Gouttenoire et al. 2020a). Besides, a GW burst signal could be significant and be the leading signal (Cui et al. 2020), in contrast to the conventional cases where the SGWB is the more sensitive means to probe cosmic strings. The dedicated study of the correlated burst and SGWB signals could reveal information about when the strings come back into the horizon as well as related inflationary dynamics.

Non-standard pre-BBN cosmology can significantly distort the vanilla shape of the SGWB from cosmic strings. It can turn it into a peaked shape (Gouttenoire et al. 2020a). Note that a peaked shape can also be generated in standard cosmology in the presence of a particle production cut-off. However, this works only for small string tension and lies outside of LISA. So, a peaked GW spectrum from cosmic strings is a smoking-gun signature of a non-standard era (less stiff than radiation, i.e. matter or inflationary era). Note that this peak is much broader than the peak from the FOPT SGWB signal and can be distinguished.

## 9.3 SGWB in quantum gravity

In Sect. 4.1 we saw some modifications to the propagation of GWs of astrophysical origin that could happen in quantum gravity. Here we will concentrate on the potential quantum-gravity effects on primordial SGWBs. (see Addazi et al. 2022 for an overview of the experimental and theoretical constraints on quantum gravity.)

Early-universe models embedded in or inspired by theories of quantum gravity have the potential to leave an observable imprint in the SGWB. In fact, these cosmological scenarios usually predict modifications in the shape of the primordial scalar and tensor spectra, which are generated either from the quantum fluctuations of an inflationary field or by an alternative mechanism. If these modifications include a blue tilt at high frequencies, they could overcome the CMB bounds and still give rise to a primordial SGWB $\Omega_{GW}(f)$ reaching the sensitivity curves of the present or future interferometers.

We can classify early universe models related to quantum gravity according to the tilt of the tensor spectrum, which we call $P_t(f)$ to distinguish the full spectrum from the exact power-law parametrization $\Delta_{h,\inf}^2$ in Eq. (99).

This classification is based on the overall trend of the full spectrum. We should note that in some models a positive tensor index $n_T > 0$ at CMB scales does not
imply a blue tilt at higher frequencies. The models are the following (see Calcagni and Kuroyanagi 2021 for a full list of references):

*   *Red-tilted tensor spectrum:* the large class of string-cosmology models with flux compactification (not to be confused with the cosmic-strings models discussed above); one semi-classical solution in Wheeler-DeWitt canonical quantum cosmology; loop quantum cosmology, where gravity is quantized canonically using the Ashtekar–Barbero connection. There are three main approaches to describe perturbations in loop quantum cosmology: with a dressed metric, in hybrid quantization and with effective constraints, the latter in turn being characterised by the presence of inverse-volume corrections, holonomy corrections or both. The first two and the third in the presence of inverse-volume corrections predict a red-tilted tensor spectrum.

*   *Blue-tilted tensor spectrum:* another semi-classical solution in Wheeler-DeWitt canonical quantum cosmology; loop quantum cosmology in the effective-constraints approach in the presence of holonomy corrections; non-local Starobinsky inflation, a model embedded in non-local quantum gravity where early-universe acceleration is driven by curvature; string-gas cosmology (again, unrelated to cosmic-strings models), where primordial spectra are generated by a thermal bath of strings most of which are wrapped around compact extra dimensions; new ekpyrotic scenario, where the spectra are generated by the collision of branes; Brandenberger–Ho non-commutative inflation, an inflationary stage realised in a geometry which is fuzzy at microscopic scales; multi-fractional inflation, an inflationary era realised in a geometry whose dimension changes with the probed scale; pre-Big-Bang scenario, where the dualities of string theory suggest the existence of a phase prior to the big bang.

Models with red-tilted tensor spectrum do not generate an observable SGWB detectable at ground- or spaced-based interferometers. Regarding models with a blue-tilted spectrum, quantum modifications in Wheeler-DeWitt canonical cosmology are too small to be observable by any of the present or planned interferometers (Calcagni and Kuroyanagi 2021). Loop quantum cosmology in the effective-constraints approach in the presence of holonomy corrections is already ruled out observationally (Bolliet et al. 2016). The tensor spectrum of non-local Starobinsky inflation (see Sect. 8.5) tends asymptotically to the spectrum of local Starobinsky inflation, hence it is unobservable even if the amplitude at interferometer scales is larger than at CMB scales (Calcagni and Kuroyanagi 2021). String-gas cosmology, the new ekpyrotic scenario, Brandenberger–Ho non-commutative inflation and multi-fractional inflation can all reach the sensitivity of DECIGO but, unfortunately, not of LISA (Calcagni and Kuroyanagi 2021).

A model worth further investigation is the *pre-Big-Bang scenario* (Gasperini and Veneziano 1993), where CMB constraints are respected and, at the same time, the amplitude of the SGWB at high frequencies can increase to touch the sensitivity curves of GW experiments, possibly including LISA (Gasperini 2016). In Gasperini (2016), the sensitivity threshold of a six-link, 5 M-km-arm-length configuration was
assumed (Caprini et al. 2016). An updated analysis could determine the parameter space of the model opening an observational window.

While the above model has a singular hypersurface at which perturbations are matched, another non-standard pre-BBN expansion history is realised by bouncing models where the universe slowly contracts (Ben-Dayan 2016; Ben-Dayan and Kupferman 2019; Artymowski et al. 2021), and then expands without a big-bang singularity or inflation (e.g. (Battefeld and Peter 2015; Lehners 2010) for reviews). In conformal time, before the bounce the scale factor goes as $a(\eta) \sim (-\eta)^b$, where $0<b \ll 1$ for slow contraction. Considering vacuum fluctuations on top of this background metric and using the standard power-law parametrization $P_h \sim f^{n_T}$, the expected SGWB spectral tilt from vacuum fluctuations is $n_{T,\text{vacuum}} \simeq 3 - |2b - 1| \simeq 2 + 2b$, resulting in a blue spectrum (Starobinsky 1979; Boyle et al. 2004; Artymowski et al. 2021) possibly slashing through the LISA sensitivity band. Hence, the SGWB signal is directly probing the geometry of the universe $a(\eta)$ making LISA the natural arena to test alternatives to inflation. As is well known, bouncing models with only one scalar field typically predict vacuum fluctuations generating spectra with $n_{T,\text{vacuum}} = n_{s,\text{vacuum}} - 1 \sim 2$. One still has to make sure that the slow-contracting bouncing models actually generate a viable scalar spectrum in accord with CMB observations, for instance by a curvaton mechanism or some other alternative. Specifically, coupling to gauge fields has been considered with an interaction term $\mathcal{L} \propto I^2(\phi)(F^2 - \gamma F \tilde{F})$, where $I$ is a function of the scalar. As seen in Sect. 8.7.1, gauge fields can provide an additional source of chiral GWs, whose tilt we denote by $n_{T,\text{sourced}}$ and is such that $n_{T,\text{sourced}} = n_{s,\text{sourced}} - 1 \approx -0.04$ in accord with CMB observations, however $r \sim 1/9$ which is above current bounds. Hence, the scalar spectrum needs to be generated by some curvaton or entropic mechanism from a second scalar field. Contrary to the inflationary case in Sect. 8.2.1, the spectral tilt $n_{T,\text{sourced}}$ is constant in the bouncing case (Ben-Dayan 2016; Ben-Dayan and Kupferman 2019). Once CMB observations are determined by a curvaton mechanism, then depending on the parameters of the model the sourced tensor spectrum may be detected by LISA with $0.15 < n_{T,\text{sourced}} < 0.31$ or $0.85 < n_{T,\text{sourced}} < 1.1$. If one wishes to match CMB observations solely with the sourced spectrum, then one open question is carrying out the calculation of the spectra across the bounce and check that the model consistently generates sourced spectra observed at CMB scales and vacuum spectra observed by LISA. This is the opposite to what happens in the inflationary models discussed in Sect. 8.2, where the tensor spectrum observed on CMB scales must be the vacuum one while LISA can only observe the sourced spectrum.

Another model worth exploring in the LISA context might be the non-perturbative gravity and bouncing universe of Biswas et al. (2006), Biswas et al. (2010), Biswas et al. (2012b), where it is possible to explore a pre-Big-Bang phase by modelling it by a string-gas dominated Hagedorn phase (Biswas et al. 2007, 2014). In this non-perturbative extension of Einstein’s gravity, there exists a bouncing phase given by the scale factor $a(t) = \cosh \lambda t$ in the presence of radiative matter and a non-zero cosmological constant. Such a solution would also permit a stiff fluid in a Hagedorn phase where the universe is primarily dominated by string winding modes (Biswas
et al. 2007, 2014). Their eventual decay into radiation, and exit from the bounce would yield the standard big-bang cosmology. Primordial perturbations are created during the thermal phase with thermal statistical initial conditions which lead, specifically, to a blue tilt in the tensor power spectrum, while keeping the amplitude of the matter spectrum within the PLANCK observation. The blue tilt in the power spectrum, along with the matter power spectrum, constrains the scale entering in the gravitational sector to be around $10^{-4} M_{\text{Pl}}$ (Biswas et al. 2014).

In parallel with the observability of these and any other model related to quantum gravity with a blue-tilted SGWB (e.g. Dapor and Liegener 2020), it will be important to assess their theoretical robustness, which we have not discussed here. In general, models directly derived from quantum gravity should yield more rigid predictions and be better falsifiable than those only inspired by quantum-gravity phenomena. A faithful and realistic description of a high-energy, high-curvature generation mechanism would move this sector of GW astronomy beyond the level of ad hoc model building towards a mini-program contributing to LISA science, according to the following algorithm:

1. Selection or construction of early-universe models embedded into or inspired by quantum gravity, under the criterion of giving a blue-tilted spectrum $P_t(f)$ of primordial tensor fluctuations.

2. Control of the underlying theoretical steps leading from the fundamental theory to the cosmological model: assumptions, approximations, parameter space, fine tunings.

3. Generation of the SGWB spectrum $\Omega_{\text{GW}}(f)$ from the primordial spectrum $P_t(f)$ via transfer functions.

4. Comparison of the theoretical SGWB with the LISA sensitivity curve and constraints on the parameter space of the model. Non-detection of a SGWB can be used to rule out theories or to constrain their parameter space. If it was detected, one could further investigate characteristics of the SGWB such as anisotropies or the local spectral tilt at LISA frequencies to extract information on the underlying physics.

In order to achieve the last goal, we will capitalise on the LISA SGWB search and use some of the pipelines described in Sect. 12.

# 10 Primordial black holes

Section coordinators: S. Clesse, J. Garcia-Bellido. Contributors: S. Clesse, V. De Luca, J.M. Ezquiaga, G. Franciolini, J. Garcia-Bellido, R. Kumar Jain, S. Kuroyanagi, I. Musco, T. Papanikolaou, M. Peloso, S. Renaux-Petel, A. Riotto, E. Ruiz Morales, M. Scalisi, O. Sergijenko, C. Unal, C. Joana, V. Vennin, D. Wands.

## 10.1 Introduction

The idea that BHs may have formed in the early universe comes back to the late 1960s with the precursor work of Zel'dovich and Novikov (1967) and to the 1970s
with the works of Hawking (1971), Carr and Hawking (1974), Carr (1975) and of Chapline (1975). Already in Carr and Hawking (1974) and Chapline (1975) it was mentioned that such PBHs could contribute to the suspected DM in the universe or to the seeds of MBHs. The first formation scenarios in the context of inflation were proposed in the 1990s (Dolgov and Silk 1993; Carr et al. 1994; García-Bellido et al. 1996) but these usually led to (evaporating) PBHs of very small mass. In the late 1990s, stellar-mass PBHs have been seriously considered as a DM candidate, following the possible detection (e.g. in the MACHO survey) of several microlensing events towards the Magellanic clouds (Aubourg et al. 1993; Alcock et al. 1997). However, the EROS (Tisserand et al. 2007) and OGLE (Wyrzykowski et al. 2010, 2011a, b; Calchi Novati et al. 2013) surveys later set more stringent limits on the PBH abundance, and at the same time, very stringent constraints from CMB observations were claimed in Ricotti et al. (2008). Figure 29 summarizes the current constraints on the PBH abundance in the idealized limit that all PBHs have the same mass.

Since 2016, the real game-changer that has rekindled the idea that PBHs may exist and constitute from a significant fraction to the totality of the DM (Bird et al. 2016; Clesse and García-Bellido 2017b; Sasaki et al. 2016) has been the first GW detection from a BH merger by Advanced LIGO/Virgo (Abbott et al. 2016b). Nowadays, the importance of the different PBH binary formation channels, the possible abundance of PBHs, and their viable mass function, are subject to an intense activity and are under discussion (for recent reviews, see e.g. Carr and Kuhnel (2020), Carr et al. (2021c)). Furthermore, since PBHs are formed by the collapse of large density perturbations, PBHs are accompanied by a SGWB sourced by these perturbations at second-order. It has been calculated that if BHs detected by LIGO/Virgo have primordial origin, there is an inevitable accompanying SGWB peaking around PTA frequencies (Ando et al. 2018; García-Bellido et al. 2017).<fnref n="31" /> In late 2020, NANOGrav has claimed the possible detection of a SGWB at nHz frequencies (Arzoumanian et al. 2020) which, among several options, can be explained as the result of the density perturbations at the origin of stellar-mass PBH formation (De Luca et al. 2021b; Kohri and Terada 2021) (although many other explanations exist. In this context, LISA will search for the GW signatures of PBHs (García-Bellido et al. 2017; Cai et al. 2019a; Bartolo et al. 2019b; Ünal 2019) and will be complementary to ground-based GW detectors (Maggiore et al. 2020; Reitze et al. 2019) and EM probes (Kashlinsky et al. 2019), in order to prove or exclude the existence of PBHs, to evaluate the possible contribution to the DM, to the seeds of MBHs at high redshift, and to distinguish PBHs from SOBHs, on a wide range of mass scales. Any firm detection would open a new window on the physics at play in the very early universe and a possible way to solve various long-standing astrophysical and cosmological puzzles (Clesse and García-Bellido 2018; Carr et al. 2021b).

<footnote n="31">³¹ This SGWB is ineludible in the sense that it does not require any further assumption other than GR and large density perturbations. It is a standard SGWB formed by anisotropic stress which is quadratic order in scalar perturbations.</footnote>
<figure><img src="image_33.png" /><figcaption>**Fig. 29** Most stringent limits on the DM fraction made of PBHs, $f_{\text{PBH}}$, coming from the Hawking evaporation producing extragalactic gamma-ray (EG$\gamma$), $e^{\pm}$ observations by Voyager 1 (Ve$^{\pm}$), positron annihilations in the Galactic Center (GCe$^{+}$), gamma-ray observations by INTEGRAL (INT), microlensing searches by Subaru HSC (HSC), MACHO/EROS (E), OGLE (O) and Icarus (I), from CMB limits (CMB), GW observations by LIGO/Virgo (LVC), wide binaries in the galactic halo (WB), the ultra-faint dwarf galaxies Eridanus II (EII) and Segue 1 (S1), X-rays towards the galactic center (XrB) and Lyman-$\alpha$ limits (Lz). For microlensing and CMB limits, the different lines indicate some degree of uncertainties, respectively due to PBH clustering and disk vs. spherical accretion. Microlensing limits only apply to the fraction of PBHs uniformly distributed in galactic halos and are less stringent if a non-negligible fraction $f_{\text{clust}}$ of PBHs are in clusters. We show the limits for $f_{\text{clust}} = 0, 0.4$ and $0.8$. For LVC, the rate suppression of early binaries still allows $f_{\text{PBH}} \gtrsim 0.1$. All these limits apply to monochromatic models and can be model dependent. Recasting them to realistic PBH models with arbitrary mass functions requires a careful analysis. Image adapted from De Luca et al. (2021b)</figcaption></figure>

This section is organized as follows: after reviewing the principal formation scenarios (Sect. 10.2), we will consider the principal sources of GWs related to PBHs and their detectability with LISA, which are the SGWB generated at second-order by the non-linear cosmological density fluctuations at the origin of PBH formation (Sect. 10.3), the PBH binaries and hyperbolic encounter (Sect. 10.4) which may unveil a possible primordial origin of MBH seeds if they are detected at high redshifts (Sect. 10.5).

## 10.2 Formation scenarios

Hereafter we first provide a rapid overview of the principal mechanisms that can lead to large curvature fluctuations and PBH formation, which can be related to the early-universe phenomena discussed in Sects. 6, 7 and 8. We then review the general theory of PBH formation from large curvature fluctuations. Finally we discuss some recent developments related to PBH formation, such as non-linear and non-Gaussian effects, thermal history, that are all relevant for the estimation of the curvature
threshold to lead to PBH formation, as well as their possible abundance and mass distribution.

### 10.2.1 Origin of curvature fluctuations

There exists a broad variety of scenarios leading to large curvature fluctuations and PBH collapse, which can be classified as below. Some scenarios are detailed in the corresponding sections of this white paper.

1. **Single-field inflation:** Slow-roll single field inflation models can produce a power spectrum amplitude of curvature fluctuations on smaller scales larger than the ones probed by CMB and large scale structures observations. The simplest example is to have two subsequent inflationary phases, the second one lasting less than 50 e-folds (roughly). Potentials with inflection points (García-Bellido and Ruiz Morales 2017; Germani and Prokopec 2017), like critical Higgs inflation (Ezquiaga et al. 2018), can also lead to a transient enhancement of the primordial power spectrum. More generically, in the slow-roll approximation, any potential leading to a transient reduction of the speed of the inflaton will lead to large curvature fluctuations, eventually leading to PBH formation. Another possibility is to invoke a variation of the sound speed during inflation, see e.g. Özsoy et al. (2018).

2. **Multi-field inflation:** Large curvature fluctuations may also arise in multi-field models, e.g. during the waterfall phase of *hybrid inflation* (Clesse and García-Bellido 2015) or due to *turning trajectories* in the inflationary landscape (Fumagalli et al. 2023). Eventually, the power spectrum will not only exhibit a broad or sharp peak on scales that are relevant for PBHs, but for a sufficiently sharp turn, this is accompanied by oscillatory features that may lead to specific signatures in the PBH population and to oscillatory patterns in the scalar-induced SGWB (Fumagalli et al. 2021a; Braglia et al. 2021; Palma et al. 2020). It is also possible that curvature fluctuations are generated by the tunneling of the inflaton towards a local minimum of the field space (Garriga et al. 2016; Deng and Vilenkin 2017; Kusenko et al. 2020) and the subsequent *bubble collapse*. Another intriguing possibility has been proposed by employing axions and their interaction with the *gauge fields* to enhance primordial density perturbations and producing PBHs (Bugaev and Klimai 2014; Linde et al. 2013; Domcke et al. 2017; García-Bellido et al. 2016, 2017; Özsoy and Lalak 2021; Özsoy 2021). In Natural Inflation (Freese et al. 1990; Adams et al. 1993), the inflaton is a pseudo-scalar particle protected from quantum corrections for super-Planckian excursions via the shift symmetry. Also, theories with UV completion predict a large number of pseudo-scalar particles that could be the inflaton or spectator fields. In this scenario, curvature fluctuations are sourced by enhanced vector modes and can have a non-Gaussian distribution.

3. **Quantum diffusion:** The backreaction of quantum fluctuations during inflation makes the dynamics of the fields stochastic, and allow them to explore wider regions of the inflationary potential. This makes the tails of the distribution functions of primordial density fluctuations much heavier (Vennin and
Starobinsky 2015; Pattison et al. 2017; Ezquiaga et al. 2020a; Figueroa et al. 2021), since they decay exponentially instead of in a Gaussian way, which boosts the production of PBHs afterwards.

4. **Curvaton and stochastic spectator fields:** PBHs could have formed due the existence of a curvaton field, e.g. with a simple modification of the original curvaton scenario (Lyth and Wands 2002). The primordial curvature perturbations on CMB scales are produced by the inflaton, which acts very similarly to the standard single-field scenario, while the curvaton field becomes responsible for perturbations on smaller scales, at the origin of PBH formation. In another scenario (Carr et al. 2021a), a stochastic spectator field experiences quantum fluctuations during inflation making it exploring a wide range of the potential, including its slow-roll part but without having any impact on the inflationary dynamics. In regions where the field acquires a value allowing slow-roll, after inflation but when these are still super-horizon, additional expansion is produced locally, which generates curvature fluctuations that later collapse into PBHs. This model has the advantage that the primordial curvature power spectrum remains at level observed at CMB scales almost everywhere, except in PBH-forming regions.

5. **Preheating:** If inflation is followed by a preheating, when the inflaton oscillates coherently at its ground state and decays to other degrees of freedom, resonant amplification of the quantum field fluctuations take place (Kofman et al. 1994, 1997). These are accompanied with a resonant amplification of curvature fluctuations (Finelli and Brandenberger 1999; Bassett et al. 1999; Jedamzik and Sigl 2000; Bassett and Viniegra 2000), which may collapse and form PBHs (Martin et al. 2020; Auclair and Vennin 2021). The PBHs that form are typically very light and only if reheating completes at very low energy (below the electroweak scale) the formed PBHs would have a relevant mass for the LISA frequency range (for GW radiation from PBH binaries). It is however also possible that these curvature fluctuations source a detectable SGWB at second order (Papanikolaou et al. 2021).

6. **Phase transitions:** The formation of PBHs may have been facilitated in FOPTs (Jedamzik and Niemeyer 1999), in non-equilibrium second order phase transitions (Rubin et al. 2000) and in specific realizations of the QCD transition (Davoudiasl 2019).

7. **Early matter era:** The curvature threshold for the PBH collapse depends on the EoS of the universe. It goes to zero in a matter dominated era, and so it is possible that PBH have formed from standard $\mathcal{O}(10^{-5})$ inflationary curvature fluctuation if the early universes has undergone a transient matter-dominated era. Several mechanisms have been proposed to be at the origin of such early matter era, see e.g. Nayak and Singh (2012), Ballesteros et al. (2020), Green et al. (1997).

8. **Cosmic strings, domain walls:** Topological defects may have lead to the production of PBHs. For instance, cosmic string cusps can collapse gravitationally into PBHs with a mass function that could extend up to stellar masses (Jenkins and Sakellariadou 2020). The collapse of domain walls, e.g. produced
by tunnelling during inflation (Deng et al. 2017; Liu et al. 2020) or in QCD axion models (Ferrer et al. 2019; Ge 2020) is another class of PBH formation scenarios.

9. **Primordial magnetic fields:** It is suspected that the required seeds for extragalactic magnetic fields have an origin in the early universe. These primordial magnetic fields induce an anisotropic stress that can act as a source of large super-Hubble curvature fluctuations, leading to PBH formation with a broad possible range of masses (Saga et al. 2020).

### 10.2.2 Standard theory of primordial black holes formation from Gaussian curvature fluctuations

If PBHs have been formed due to the collapse of large primordial curvature perturbations, the fraction of the density of the universe made of PBHs at the time of their formation, $\beta_{\text{form}}$ (usually defined per unit of logarithmic mass interval), is determined by the probability that the amplitude of a primordial curvature fluctuation, measured in terms of the mass excess of the density contrast $\delta$ is above a certain threshold $\delta_c$ when the perturbation re-enter the cosmological horizon, corresponding to the Hubble length $1/H$, where $H$ denotes the Hubble parameter. The corresponding PBH mass $M$ is therefore linked to the Hubble horizon mass of the collapsing fluctuation, and so also to its size and its corresponding wavenumber $k_M = aH$ (and $a$ is the scale factor),

$$M = \gamma M_{\text{H}} = \left( \frac{3H^2}{8\pi G} \right) \frac{c^3}{H^3}, \qquad (140)$$

where $\gamma$ is the ratio between the final PBH mass and the collapsing Horizon mass, which depends on the formation model details (a realistic range is 0.1 ÷ 1). This should take into account that, when the perturbation amplitude $\delta$ is close enough to the threshold $\delta_c$ (i.e. $(\delta - \delta_c) \lesssim 10^{-2}$), the mass follows the scaling law of critical collapse:

$$M = K(\delta - \delta_c)^{\Gamma} M_{\text{H}}, \qquad (141)$$

where $\Gamma$ depends only on the EoS ($\Gamma \simeq 0.36$ in the radiation dominated universe) and $K$ depends on the shape of the perturbation (typically 1 ÷ 10).

If the curvature fluctuations originate from inflation, the PBH mass can be related to the time of Hubble exit of the corresponding scale during inflation, expressed in terms of the number of e-folds before the end of inflation $N_*$ (see e.g. García-Bellido et al. 1996), which gives, in the case of instantaneous reheating,

$$M \sim 4\pi \frac{M_{\text{Pl}}^2}{H_{\text{inf}}} e^{2N_*}, \qquad (142)$$

where $M_{\text{Pl}}$ is the reduced Planck mass and $H_{\text{inf}}$ is the (almost constant) Hubble rate during inflation.

In most inflationary scenarios, the distribution of primordial curvature fluctuations is Gaussian, or almost Gaussian. Some recent scenarios have proposed more complex
distributions that can be useful to alleviate or reduce the need of a large power spectrum amplitude. These are shortly discussed later and assuming that the curvature fluctuations are described by a Gaussian distribution the fraction of PBHs at formation is then usually given by

$$ \beta_{\text{form}}(M) = 2 \int_{\delta_c}^{\infty} \frac{1}{\sqrt{2\pi}\sigma} e^{-\frac{\delta_c^2}{2\sigma^2}} d\delta = 1 - \operatorname{erf}\left(\frac{\delta_c}{\sqrt{2}\sigma}\right) = \operatorname{erfc}\left(\frac{\delta_c}{\sqrt{2}\sigma}\right). \quad (143) $$

This expression however does not take into account that when $\delta$ is larger then a certain value $\delta_{\text{max}}$ (originally estimated $\delta_{\text{max}} \sim 1$ by Bernard Carr in 1975 (Carr 1975)) the perturbation forms a separate closed universe, topologically disconnected. A more accurate version can nevertheless be obtained re-normalizing the previous expression as

$$ \beta_{\text{form}}(M) = \frac{\int_{\delta_c}^{\delta_{\text{max}}} \frac{1}{\sqrt{2\pi}\sigma} e^{-\frac{\delta_c^2}{2\sigma^2}} d\delta}{\int_{0}^{\delta_{\text{max}}} \frac{1}{\sqrt{2\pi}\sigma} e^{-\frac{\delta_c^2}{2\sigma^2}} d\delta} = 1 - \frac{\operatorname{erf}\left(\frac{\delta_c}{\sqrt{2}\sigma}\right)}{\operatorname{erf}\left(\frac{\delta_{\text{max}}}{\sqrt{2}\sigma}\right)} = \frac{\operatorname{erfc}\left(\frac{\delta_c}{\sqrt{2}\sigma}\right) - \operatorname{erfc}\left(\frac{\delta_{\text{max}}}{\sqrt{2}\sigma}\right)}{1 - \operatorname{erfc}\left(\frac{\delta_{\text{max}}}{\sqrt{2}\sigma}\right)}. \quad (144) $$

The variance of the field of density perturbations $\sigma$, according to the Gaussian distribution of $\delta$ used before, is given by

$$ \sigma^2 = \langle \delta^2 \rangle = \int_0^\infty \frac{dk}{k} \mathcal{P}_\delta(k, r) = \frac{16}{81} \int_0^\infty \frac{dk}{k} (kr)^4 \tilde{W}^2(k, r) T^2(k, r) \mathcal{P}_\zeta(k), \quad (145) $$

where $\mathcal{P}_\delta(k, r)$ and $\mathcal{P}_\zeta(k)$ are the density and the curvature power spectrum, while $\tilde{W}(k, r) = 3(\sin kr - kr \cos kr)/(kr)^3$ is the Fourier transform of the top-hat smoothing function and $T(k, r) = \tilde{W}(k, r/\sqrt{3})$ is the linear transfer function, both computed at the horizon crossing scale. All this shows that a larger power spectrum could increase significantly the fraction of PBHs because the abundance of PBHs is exponentially sensitive to the value of the amplitude.

One can then calculate the contribution of PBHs to the density of the universe today, in terms of the fraction of DM that they represent, given by

$$ f_{\text{PBH}}(M) \sim 2.4 \times \beta_{\text{form}}(M) \sqrt{\frac{M_{\text{eq}}}{M}}, \quad (146) $$

where $M_{\text{eq}} \simeq 2.8 \times 10^{17} M_{\odot}$ is the horizon mass at matter-radiation equality, and the numerical factor corresponds to $2 \times (1 + \Omega_b/\Omega_{\text{CDM}})$. In this standard scenario, where PBHs form from Gaussian curvature perturbations, one does not expect spatial clustering at formation larger than the one predicted by the Poisson distribution (Ali-Haïmoud 2018; Desjacques and Riotto 2018; Ballesteros et al. 2018; Moradinezhad Dizgah et al. 2019).
### 10.2.3 Beyond Gaussianity

The above formalism applies to Gaussian perturbations but it can be generalized to non-Gaussian distributions, by replacing the integrand in the left-hand side of Eq. (144) with the appropriate distribution function.

At the perturbative level, even a small amount of local non-Gaussianity (Lyth 2012; Byrnes et al. 2012), or the inevitable non-linear (and hence non-Gaussian) relation between the primordial curvature perturbations and density perturbations can have important effects (Young and Byrnes 2015; Franciolini et al. 2018; De Luca et al. 2019; Young et al. 2019). In the regime of large non-Gaussianity (Linde et al. 2013; García-Bellido et al. 2016; Ünal et al. 2021), the amplitude of the power spectrum can produce same abundance of PBHs even if it is many orders of magnitude smaller with respect to Gaussian case. The existence of non-Gaussianity also influences the amount of SGWB accompanying the PBHs (Nakama et al. 2017; García-Bellido et al. 2017; Cai et al. 2019a; Ünal 2019). The fact that the power spectrum needs to deviate from the quasi scale invariance observed on CMB scales to reach fluctuations allowing PBHs to form also implies that it is natural to expect the non-Gaussianity to be strongly scale dependent too. Moreover, the merger rate of PBHs is strongly impacted by non-Gaussianity (Young and Byrnes 2020; Atal et al. 2020), because non-Gaussianity couples the long and short wavelength primordial density perturbations, leading to initial spatial clustering of the PBHs.

Since PBHs require large fluctuations to form, a perturbative description may not be sufficient. Quantum diffusion of the inflationary fields leads to large deviations from Gaussian statistics on the tails of the distribution functions of primordial density fluctuations (Pattison et al. 2017; Biagetti et al. 2018; Ezquiaga and García-Bellido 2018; Figueroa et al. 2021, 2022), which acquire an exponential (rather than Gaussian) profile (Vennin and Starobinsky 2015; Pattison et al. 2017; Ezquiaga et al. 2020a; Figueroa et al. 2021, 2022). It implies that PBHs can be formed with a much smaller amplitude of the power spectrum than what would be inferred using Gaussian statistics. This cannot be properly taken into account with the usual, perturbative parametrisation of non-Gaussian statistics (such as those based on computing the few first moments of the distribution and the non-linearity parameters $f_{NL}$, $g_{NL}$, etc.), which can only account for polynomial modulations of Gaussian tails, and needs to be described with a non-perturbative approach such as the stochastic-$\delta N$ formalism (Fujita et al. 2013; Vennin and Starobinsky 2015).

### 10.2.4 The threshold for primordial black holes and the impact of thermal history

A crucial parameter of the formalism presented above is the critical value $\delta_c$, distinguishing between cosmological perturbations collapsing into PBHs ($\delta > \delta_c$) and those ones bouncing back into the surrounding medium ($\delta < \delta_c$). This is a fundamental parameter because the resulting PBH abundance is exponentially sensitive to its value. The analysis of the gravitational collapse of curvature perturbations to form PBHs and the appropriate threshold condition has been an active line of research in the past years (Kopp et al. 2011; Harada et al. 2013; Young et al. 2014). It has been estimated using analytical methods (Harada et al. 2013) but
the best approach is to use fully relativistic simulations of PBH formation in spherical symmetry (Musco et al. 2005). Important results having emerged from recent studies are that its exact value is impacted by non-linear and non-Gaussian effects, and that it depends on the radial profile of the overdensity (Musco 2019; Young et al. 2019; De Luca et al. 2019; Kehagias et al. 2019) as well as on the shape of the primordial power spectrum (Yoo et al. 2018; Germani and Musco 2019).

Very recently, a new semi-analytical method tested against simulations in numerical relativity has been proposed in Musco et al. (2021), computing $\delta_c$ from the shape of the power spectrum, applied to a few particular cases (power-law spectra, log-normal or Gaussian peak, ...). In the radiation dominated universe the typical range of the threshold lies within $0.4 < \delta_c < 0.6$, with the larger values corresponding to a more peaked shape of the peak of the power spectrum.

On super horizon scales the non-linear amplitude of the curvature profile $\zeta$ is important for the formation of PBHs, and the energy density contrast $\delta\rho/\rho_b$, when the universe is radiation dominated, is expressed in terms of $\zeta$ as

$$ \frac{\delta\rho}{\rho_b} = \frac{\rho(r,t) - \rho_b(t)}{\rho_b(t)} = -\frac{1}{a^2 H^2} \frac{8}{9} e^{-5\zeta(r)/2} \nabla^2 e^{\zeta(r)/2}. \quad (147) $$

It can be shown that the amplitude $\delta$ is a quadratic function of the curvature profile (see for example (Musco 2019) for more details)

$$ \delta = -\frac{2}{3} \zeta'(r_m) [2 + r_m \zeta'(r_m)] = \delta_G \left(1 - \frac{3}{8} \delta_G\right), \quad (148) $$

where $\delta_G \equiv -\frac{4}{3} r_m \zeta'(r_m)$ is the linear component of the amplitude $\delta$. The value of $r_m$ is defined by the location of the peak of the compaction function $C \equiv 2\Delta M/R$ (where $\Delta M = M - M_b$ is the excess of mass with respect the background) and is given by $\zeta(r_m) + r_m \zeta'(r_m) = 0$. Given the value of $\delta_c$ the threshold value of the linear component $\delta_G$ is included within $0.5 \lesssim \delta_{c,G} \lesssim 0.9$, and from Eq. (148) we see that $\delta_{max} = 8/3$, corresponding to the maximum value of $\delta_G$, above which $\delta$ becomes negative, and does not describe a cosmological perturbation of our universe.

The threshold $\delta_c$ is also sensitive to the EoS at the time of formation. For example, the QCD phase transition makes the EoS to drop, increasing the production of PBHs of mass $O(M_\odot)$ (Jedamzik 1997; Byrnes et al. 2018). The reheating at the end of inflation should have filled the universe with radiation. In the absence of extensions beyond the SM, the universe remains dominated by relativistic particles with an energy density decreasing as the fourth power of the temperature as the universe expands. The number of relativistic degrees of freedom remains constant ($g_* = 106.75$) until around 200 GeV, when the temperature of the universe falls to the mass thresholds of SM particles.

The first particle to become non-relativistic is the top quark at $T \simeq m_t = 172$ GeV, followed by the Higgs boson at 125 GeV, and the Z and W bosons at 92 and 81 GeV, respectively. These particles become non-relativistic at nearly the same time and this induces a significant drop in the number of relativistic degrees of freedom down to $g_* = 86.75$. There are further changes at the b and c quark and $\tau$-lepton thresholds but these are too small to appear in Fig. 30. Thereafter $g_*$ remains approximately
<figure><img src="image_34.png" /><figcaption>**Fig. 30** Left panel: Evolution of the relativistic degrees of freedom $g_*$ as a function of the temperature. The grey vertical lines correspond to the masses of the electron, pion, proton/neutron, W, Z bosons and top quark, respectively. Right panel: Effect of the evolution of $g_*$ on the critical overdensity $\delta_c$ leading to PBH formation, as a function of the Hubble horizon mass (related to the PBH mass by $M = \gamma m_H$)</figcaption></figure>

constant until the QCD transition at around 200 MeV, when protons and neutrons condense out of the free light quarks and gluons. The number of relativistic degrees of freedom then falls abruptly to $g_* = 17.25$. A little later the pions become non-relativistic and then the muons, giving $g_* = 10.75$. Thereafter, $g_*$ remains constant until $e^+e^-$ annihilation and neutrino decoupling at around 1 MeV, when it drops to $g_* = 3.36$.

Whenever the number of relativistic degrees of freedom suddenly drops, it changes the effective EoS parameter $w$. There are thus four periods in the thermal history of the universe when $w$ decreases. After each of these, $w$ resumes its relativistic value of 1/3 but each sudden drop modifies the probability of gravitational collapse of any large curvature fluctuations present at that time, as shown in Fig. 30.

As illustrative examples, we have computed the PBH mass functions for two models with an (almost) scale-invariant power spectrum and two different values of the spectral index, $n_s = 0.97$ (Model 1) and $n_s = 1$ (Model 2). We assumed $\gamma = 0.8$ in both cases. The imprints of the thermal history on the PBH mass function are clearly visible. It is worth noticing that these features rely on known physics and are therefore unavoidable for any PBH model with a wide mass function. The former case corresponds to the scenario proposed in Carr et al. (2021b), Clesse and Garcia-Bellido (2022) and the latter in Byrnes et al. (2018), De Luca et al. (2021b). They can both account for the totality of the DM and somehow explain some LIGO/Virgo GW events, but produce different abundances in the stellar mass range: $f_{\text{PBH}}(M_{\odot}) \approx 0.8$ in the first case, $f_{\text{PBH}}(M_{\odot}) \approx 10^{-4}$ in the second case where the peak lies in the sublunar range. We stress that the second example avoids the bounds in the LIGO/Virgo range.

## 10.3 Stochastic gravitational-wave background sourced at second order by curvature fluctuations

If PBHs are generated by the collapse of large density perturbations, they are unavoidably associated to the emission of induced GWs at second order by the same scalar perturbations due to the intrinsic nonlinear nature of gravity (Acquaviva et al.
2003; Mollerach et al. 2004). The phenomenological implications have been investigated in various contexts also associated to PBHs Ananda et al. (2007), Baumann et al. (2007), Bugaev and Klimai (2010), Saito and Yokoyama (2010), García-Bellido et al. (2017), Ando et al. (2018), Bartolo et al. (2018), Bartolo et al. (2019c), Bartolo et al. (2019b), Clesse et al. (2018), Ünal (2019), Chatterjee and Mazumdar (2018), Wang et al. (2019), Domènech (2020), Domènech et al. (2020), Pi and Sasaki (2020), Ragavendra et al. (2021), Fumagalli et al. (2021a). If the enhancement of the scalar power spectrum responsible for the generation of PBHs occurs around characteristic scales associated with frequencies between $10^{-7}$ and $10^{-2}$ Hz, this SGWB becomes detectable by GW experiments like LISA. It is worth emphasizing that contrary to the PBH abundance that is exponentially sensitive to the power spectrum, this SGWB depends on the power spectrum amplitude to the second power. This way, LISA will even be able to exclude the existence of an extremely tiny fraction of DM made of PBHs (even a single PBH in our universe) (Clesse et al. 2018), within a wide mass range.

Figure 31 presents the PBH density fraction at formation $\beta_{form}$ (left panel) and the corresponding PBH mass function $f_{PBH}$ today (right panel) for two models with a power-law power spectrum (see the caption of the Figure for details). The SGWB associated with one of these two models is shown in Fig. 32, where it is confronted with several current or forecasted experimental limits. The SGWB covers a wide frequency range. In the ultra-low frequency range, around nHz, PTA experiments like PPTA (Shannon et al. 2015), NANOGrav (Arzoumanian et al. 2018) and EPTA (Lentati et al. 2015) give the most stringent constraints on the GWs abundance. Future experiments like SKA (Dewdney et al. 2009) (see also Moore et al. 2015) will greatly improve the sensitivity. In the LIGO/Virgo frequency range, an additional constraint has been set by the non-observation of a SGWB after O1-O2 (Abbott et al. 2019c) and O3 runs (Abbott et al. 2021g). All these searches can be translated into a constraint on the amplitude of the comoving curvature perturbation at the

<figure><img src="image_35.png" /><figcaption>Fig. 31 PBH density fraction at formation $\beta_{form}$ (left panel) and the corresponding PBH mass function $f_{PBH}$ today (right panel), neglecting the effect of PBH growth by accretion and hierarchical mergers, for two models with a power-law power spectrum and including the effects of thermal history: Model 1 from Carr et al. (2021b), Clesse and García-Bellido (2022) with spectral index $n_s = 0.97$; Model 2 from De Luca et al. (2021b), Byrnes et al. (2018) with $n_s = 1$. and a cut-off mass of $10^{-14} M_{\odot}$. The transition between the large-scale and small-scale power spectrum is fixed at $k = 10^3 \text{ Mpc}^{-1}$. The power spectrum amplitude is normalized such that both models produce an integrated PBH fraction $f_{PBH} = 1$, i.e. PBH constitute the totality of DM. A value of $\gamma = 0.8$ was assumed</figcaption></figure>
<figure><img src="image_36.png" /><figcaption>**Fig. 32** SGWB sourced at second order by the density perturbations at the origin of PBH formation, for Model 2 of Fig. 31. On top of the plot, we show the PBH mass associated to a given GW frequency as in Eq. (149). The LISA sensitivity (Amaro-Seoane et al. 2017) and the hint for a detection by NANOGrav 12.5 yr (Arzoumanian et al. 2020) are represented, as well as the constraints coming EPTA (Lentati et al. 2015), PPTA (Shannon et al. 2015), NANOGrav 11 yrs (Arzoumanian et al. 2018; Aggarwal et al. 2019) and future sensitivity curves for planned experiments like SKA (Zhao et al. 2013), DECIGO/BBO (Yagi and Seto 2011), CE (Abbott et al. 2017b), ET (Hild et al. 2011), Advanced LIGO + Virgo collaboration (Abbott et al. 2017f), Magis-space (MS) and Magis-100 (M100) (Coleman 2019), AEDGE (El-Neaj et al. 2020) and AION (Badurina et al. 2020). Image reproduced with permission from De Luca et al. (2021b)</figcaption></figure>

corresponding scales (Bugaev and Klimai 2010; Byrnes et al. 2019; Inomata and Nakama 2019; Ünal et al. 2021). Those bounds are also affecting the maximum allowed PBHs fraction of DM with the hypothesis that they originate from the collapse of density perturbations. Detailed studies with the LIGO/Virgo data affecting the mass range $M \in [10^{-20}, 10^{-18}]M_{\odot}$ are reported in Kapadia et al. (2020), while very tight bounds in the mass range $M \in [10^{-3}, 1]M_{\odot}$ are obtained in Chen et al. (2020) using the latest NANOGrav data; see also Cai et al. (2019b) where the dependence of the result to non-Gaussianities is also investigated, finding that local non-Gaussianity can for example alleviate the bounds (see Sect. 10.2.3 for details). Finally, the next generation multimessenger experiments, CMB distortion (PIXIE) and PTA-SKA, can test the PBH scenario over solar mass robustly, namely they can conclusively detect or rule out the PBHs over solar-mass and the intriguing proposal that the seeds of the MBHs are formed by PBHs (Ünal et al. 2021) independent of i) statistical properties of perturbations, ii) accretion and merger history and iii) clustering effects.

LISA will be able to provide insights in the intermediate frequencies, and corresponding masses. Since the emission mostly comes when the corresponding scales cross the horizon, one can relate the GWs frequency to the PBHs mass $M$ as (see for example (Saito and Yokoyama 2010; García-Bellido et al. 2017))

$$f \approx 3 \text{ mHz} \left(\frac{\gamma}{0.2}\right)^{1/2} \left(\frac{M}{10^{-12} M_{\odot}}\right)^{-1/2}, \qquad (149)$$

where the factor $\gamma$ is capturing the relation between the horizon mass at formation and
the PBH mass after the collapse. Notice that the peak frequencies fall within the LISA sensitivity band for PBH masses around $M \sim \mathcal{O}(10^{-15} \div 10^{-8}) M_{\odot}$ and for this mass range, the PBHs can constitute the totality of the DM. Hence, García-Bellido et al. (2017) proposed PBHs in this mass range as DM and further found that density perturbations forming PBHs lead to GWs detectable by LISA. This proposal has been studied in more detail in Cai et al. (2019a), Bartolo et al. (2019b), Bartolo et al. (2019c), Ünal (2019).

The computation of the resulting SGWB spectrum was originally performed in Ananda et al. (2007). We provide here the main result, assuming a generic form for the power spectrum of curvature fluctuation. The current GW abundance can then be obtained as

$$ \Omega_{\text{GW}}(\eta_0, k) = \frac{a_f^4 \rho_{\text{GW}}(\eta_f, k)}{\rho_r(\eta_0)} \Omega_{r,0} = \frac{g_*(\eta_f)}{g_*(\eta_0)} \left( \frac{g_{*S}(\eta_0)}{g_{*S}(\eta_f)} \right)^{4/3} \Omega_{r,0} \Omega_{\text{GW}}(\eta_f, k), \quad (150) $$

in terms of the present radiation energy density fraction $\Omega_{r,0}$ if the neutrinos were massless. The crucial quantity is $\Omega_{\text{GW}}(\eta_f, k)$, that is the fractional GW energy density for log interval at the emission epoch $\eta_f$, related to the critical energy density of a spatially flat universe $\rho_c = 3H^2 M_p^2$. Assuming that the scalar perturbations $\zeta$ are Gaussian, it can be calculated as

$$ \begin{aligned} \langle \rho_{\text{GW}}(\eta, \vec{x}) \rangle &\equiv \rho_c(\eta) \int d \ln k \, \Omega_{\text{GW}}(\eta, k) \\ &= \frac{2\pi^4 M_p^2}{81\eta^2 a^2} \int \frac{d^3 k_1 d^3 p_1}{(2\pi)^6} \frac{1}{k_1^4} \frac{\left[ \frac{p_1^2 - (\vec{k}_1 \cdot \vec{p}_1)^2}{k_1^2} \right]^2}{p_1^3 |\vec{k}_1 - \vec{p}_1|^3} \quad (151) \\ \mathcal{P}_{\zeta}(p_1) \mathcal{P}_{\zeta}(|\vec{k}_1 - \vec{p}_1|) &\left( \mathcal{I}_c^2(\vec{k}_1, \vec{p}_1) + \mathcal{I}_s^2(\vec{k}_1, \vec{p}_1) \right), \end{aligned} $$

where the functions $\mathcal{I}_{c,s}$ are found in Espinosa et al. (2018), Kohri and Terada (2018). The integrals need to be done numerically for general power spectra (see Saito and Yokoyama 2010; Bugaev and Klimai 2010 for analytical calculations in the specific case of a monochromatic or Gaussian curvature spectrum).

For the frequencies of interest, using $f \simeq 8$ mHz($g_*/10)^{1/4}(T/10^6$GeV), one can show that the emission of GWs takes place at $\eta_f$ well before the time at which top quarks start annihilating, above which we can assume a RD universe with constant effective degrees of freedom.

The non-linear coupling with the curvature perturbation naturally leads to an intrinsically non-Gaussian GWs signal imprinted in phase correlations. However, the coherence is washed out by the propagation of the waves in the perturbed universe mainly due to time delay effects originated from the presence of large scale variations of gravitational potential (Bartolo et al. 2018, 2019b, c; Margalit et al. 2020). This is simply a consequence of the central limit theorem applied to a number
$N \sim (k_{\star} \eta_0)^2 \gg 1$ of independent lines of sight (Bartolo et al. 2019b), $k_{\star}$ being the characteristic perturbation wave-number roughly proportional to the inverse horizon size at GW emission. Possible small deformations smearing the GWs spectrum can also arise from similar effects (Domcke et al. 2020c).

An interesting primordial signal that is potentially observable is related to scenarios where the scalar power spectrum presents oscillations of sufficiently large amplitude, characteristic of large particle production mechanisms, leading to oscillatory $O(10\%)$ modulations in the frequency profile of the scalar-induced SGWB (Fumagalli et al. 2021a; Braglia et al. 2021), see Sect. 8.3 for details.

## 10.4 Resolved sources and stochastic gravitational-wave backgrounds from primordial black holes binaries and hyperbolic encounters

PBHs can source GWs in several ways, which will be probed as individual and resolved sources, as a SGWB, or continuous wave signals with LISA. One can distinguish:

* GW chirps from the merging of PBH binaries

* SGWB from PBH binaries

* Continuous GWs from PBH binaries far from the merging time.

* Bursts from close BH-BH interactions

The amplitude of these signals depends on the PBH mass function and the resulting merger rates, which themselves depend on the PBH binary formation channel. Two main channels have been identified: at formation in the early universe and by tidal capture in dense PBH clusters. For each we provide an estimation of the corresponding merger rates. We then compute the expected signals for a few example models. For the purpose of calculating the PBH merging/encounter rates and the resulting GW signatures, we use a general mass distribution $f(m_{\text{PBH}})$ that would be specified by the underlying formation model.

### 10.4.1 Binary formation channels and merger rates

**Primordial binaries** may have been created before the epoch of matter-radiation equality, when two PBHs formed sufficiently close to each other for their dynamics to be independent of the expansion the universe. If this binary takes a time of the order of the age of the universe to merge, the resulting GWs could be detected. One can estimate the merger rate of such binaries (Kocsis et al. 2018; Raidal et al. 2019) in the simplest scenario where PBHs are not clustered at formation and only contribute to less than about 10% of the DM as
$$ \frac{dR^{\text{prim}}}{d(\ln m_1)d(\ln m_2)dz} = \frac{1.6 \times 10^6}{\text{Gpc}^3 \text{yr}} f_{\text{PBH}}^{\frac{53}{37}} S(m_1, m_2, f_{\text{PBH}}) \left[ \frac{t m_1 m_2}{t_0 (m_1 + m_2)^2} \right]^{-\frac{34}{37}} \left( \frac{m_1 + m_2}{M_{\odot}} \right)^{-\frac{32}{37}} f(m_1) f(m_2) \quad (152) $$

where $m_1$ and $m_2$ are the BH masses and $S$ accounts for the suppression factor of the merger rate coming from binary disruption in early universe substructures (Raidal et al. 2019). We use the definition of the PBH mass function $f(m)$ normalised to unity as $\int f(m)d\ln m = 1$.

Notice that accretion onto PBHs in binaries can be effective for masses above $\mathcal{O}(10) M_{\odot}$ with an impact on the merger rate (De Luca et al. 2020d, c). When $f_{\text{PBH}} \gtrsim 0.1$, N-body simulations (Raidal et al. 2019) have shown that PBH clusters can rapidly form and change the lifetime of PBH binaries, due to their tidal force. However, the fraction of unperturbed binaries at the present time is at least of the order of $10^{-2}$ for $f_{\text{PBH}} = 1$ (Vaskonen and Veermäe 2020; De Luca et al. 2020a; Jedamzik n.d.). Furthermore, one should also recall that not all the binaries end up inside halos, about $10^{-3}$ to $10^{-2}$ remain isolated. For instance, this gives a primordial merger rate still above the inferred LIGO/Virgo rate, even for large PBH abundances.

It is worth mentioning that, as investigated in (De Luca et al. 2021a), an overall PBH abundance of the order $f_{\text{PBH}} \sim 10^{-4}$ would be enough to explain the rate of events in the pair-instability mass gap as inferred by the LIGO/Virgo observations (Abbott et al. 2020b). This result depends on the rate of PBH accretion which is needed to avoid CMB constraints in that relevant mass window (De Luca et al. 2020b). The redshift dependence of the GW source merger rate at high redshift is a key feature which can be explored to distinguish between the astrophysical BHs and PBHs from LIGO/Virgo (Mukherjee and Silk 2021) and CE and ET (Mukherjee et al. 2022). LISA will be able to observe events in that mass window (Kaiser and McWilliams 2021) and thus potentially help in constraining the population properties of events falling in the pair-instability mass gap (see Abbott et al. 2020c and references therein).

**Capture in PBH clusters:** The second most important PBH binary formation channel is through dynamical capture in dense PBH clusters. Therefore the PBH clustering properties, which are still uncertain and highly model-dependent, mostly determine their overall merger rate. In the absence of additional clustering compared to $\Lambda$CDM expectations, one has a rate $R \approx \mathcal{O}(1) \text{yr}^{-1} \text{Gpc}^{-3}$ for equal-mass binaries (Bird et al. 2016). But in the case of wide mass functions, clustering may be boosted as well as the merger rates, up to the values inferred by LIGO/Virgo if stellar-mass PBHs constitute a significant fraction of the DM. One can encode the clustering uncertainties in a single, rescaling parameter $R^{\text{clust}}$ (Vaskonen and Veermäe 2020; De Luca et al. 2020a). Then, the rate distribution can be approximated as (Bird et al. 2016; Clesse and García-Bellido 2017b; Sasaki et al. 2018)
$$ \frac{dR^{\text{capt}}}{d \ln m_1 d \ln m_2} \approx \text{Gpc}^{-3} \text{yr}^{-1} f_{\text{PBH}} R^{\text{clust}} \frac{(m_1 + m_2)^{10/7}}{(m_1 m_2)^{5/7}} f(m_1) f(m_2). \quad (153) $$

One has to notice though that both in the LIGO/Virgo and LISA mass ranges there are various severe constraints on the PBH abundance (Carr et al. 2021c), which are not sensitive to uncertainties on the PBH clustering, e.g. those coming from the CMB anisotropies (Ali-Haïmoud and Kamionkowski 2017; Poulin et al. 2017; Inman and Ali-Haïmoud 2019; Serpico et al. 2020).

### 10.4.2 Gravitational-wave chirps from primordial black holes binaries

The merger rates expected for the two illustrative PBH models considered here are represented on Fig. 33, for the two binary formation channels detailed in the previous section. From these rates, one can identify three types of GW sources that are relevant for LISA.

Intermediate-mass binaries: The LISA frequency range is ideal to detect IMBBH mergers with mass above $10^3 M_{\odot}$. PBHs can produce such merger events with rates of order $\sim 10^{-2} \text{yr}^{-1} \text{Gpc}^{-3}$ for PBHs in clusters, which could therefore be detected by LISA given the expected astrophysical reach for such events, typically above Gpc. The rate of equal-mass primordial binaries roughly goes like $\sim M^{-1}$ for primordial

<figure><img src="image_37.png" /><figcaption>Fig. 33 Expected merger rates of PBH of mass $m_1$ and $m_2$ for the Model 1 (top panels) and Model 2 (bottom panels) mass distributions displayed in Fig. 31 due to the binary formation channels “primordial binaries” (left panels) and “tidal capture in halos” (right panels) coming from Eqs. (152) and (153) respectively</figcaption></figure>
binaires and is thus relatively low, while it is constant for binaries in clusters. Merger rates are also much larger for Model 1 than for Model 2, due to a much larger PBH density fraction in this range. It must be pointed out that this fraction is restricted by CMB limits. These are not calculated for such wide-mass distributions of highly clustered PBHs, and so they should not be directly applied to the models.

Extreme mass ratios: Another distinctive feature of PBHs with wide mass distributions is to predict large merger rates of binaries with extreme mass ratios, especially those involving a PBH from the QCD peak, merging with an intermediate mass PBH, as shown in Fig. 33. For Model 1 and tidal capture in clusters, these rates can be larger than $10^3 \text{yr}^{-1} \text{Gpc}^{-3}$. Therefore LISA could detect such extreme mass ratio coalescences, because of the ideal sub-Hertz GW frequency of such mergers and the large merger rates that allow for the compensation of the reduced strain sensitivity of LISA compared to the one of ground-based detectors at higher frequencies.

Continuous waves from galactic binaries: Models with an important PBH fraction at planetary-masses and below typically lead to very high merger rates at low mass, which could allow for the detection of such binaries within our galaxy. This is especially motivated for the Model 2 and for primordial binaries. But in the LISA range, such binaries are still far from the time of merger, when the strain evolution is almost constant over the duration of an observing run. For these kind of signals the continuous wave searches (see e.g. Miller et al. 2021) look promising. The sensitivity of LISA to planetary-mass PBHs using such search techniques is still to be determined.

### 10.4.3 Stochastic gravitational-wave background

The overlap of GWs from PBH binaries that are close to merger form together a SGWB that could be detected by LISA. The possible spectral shape of this SGWB for the two possible binary formation channels, the constraints coming from observations and their implications for PBH models, have been discussed in the recent literature (Mandic et al. 2016; Clesse and García-Bellido 2017a; Wang et al. 2018; Raidal et al. 2017; Chen et al. 2019; Wang et al. 2019). Hereafter, we review the basic principle and formula to compute this SGWB and present some predictions for particular PBH models.

The amplitude of the SGWB is given by summing up the energy spectrum of each binary system by taking into account the merger rate distribution as well as its possible evolution with redshift. For binary systems with a circular orbit and for a generic PBH mass function leading to merger rates $R(m_1, m_2)$, e.g. given by Eq. (152) for primordial binaries and by Eq. (153) for capture in halos, one gets (Clesse and García-Bellido 2017a)

$$ \Omega_{\text{GW}}(f) = \frac{8\pi^{5/3} G^{5/3}}{9H_0^2} f^{2/3} \int_0^{z_{\text{max}}} dz \frac{1}{H(z)(1+z)^{4/3}} \times \int d\ln m_1 \int d\ln m_2 R(m_1, m_2, z) M_c^{5/3}(m_1, m_2) \quad (154) $$

where $M_c$ is the chirp mass. The maximum frequency in the observer's frame is
determined by the innermost stable circular orbit which is given by $f_{\text{ISCO}} \approx 4.4\,\text{kHz} M_{\odot}/((m_1 + m_2)(1+z))$. This can be translated to the maximum redshift as $z_{\text{max}} = f/f_{\text{ISCO}} - 1$.

As an example, we used a peaked mass distribution around $2.5 M_{\odot}$ (motivated by the QCD peak) and computed the corresponding SGWB, for binaries formed in halos. This SGWB is clearly detectable by LISA but probably below the sensitivity of LIGO/Virgo. By probing the amplitude and spectral index of this SGWB, LISA would be able to distinguish between different mass functions, different clustering properties, as well as the dominant binary formation channel. Indeed, the lower limit of the frequency range covered by PBH binaries depends on their relative velocity when the binary is formed, which itself depends on the typical cluster mass and radius.

### 10.4.4 Bursts from hyperbolic encounters

If PBHs are at present grouped in clusters, a fraction of BH encounters will not end up producing bound systems, but a single scattering event. This could happen e.g. when the relative velocity is too high for the capture to be possible. The GW signals from hyperbolic encounters can have frequencies, strain amplitudes and characteristic time durations that create GWs signals that can be detected by LISA. These events would have unique signatures (García-Bellido and Nesseris 2018, 2017) and would provide direct information about the orbital parameters and spatial distribution of these BHs, thereby providing complementary information to inspiral binaries, and strong evidence in favor of the scenario of clustered PBHs (Fig. 34).

The waveforms of the GW emission in hyperbolic encounters are very different from those of the inspiraling binaries, since the majority of the energy is released near the point of closest approach, generating a burst of GWs with a characteristic “tear-drop” shape of the emission in the time-frequency domain (García-Bellido and Nesseris 2018, 2017). The burst has a characteristic peak frequency

$$f_{\text{peak}} = 0.32 \, \text{mHz} \times \frac{\beta(e+1) \, \text{AU}}{(e-1) \, b}, \qquad (155)$$

which corresponds to the maximum GW amplitude and depends only on the impact parameter $b$, the total mass of the system $M$ and the eccentricity of the hyperbolic orbit $e = \sqrt{1 + b^2 v_0^4 / G^2 M^2}$, where $v_0$ is the asymptotic relative velocity of the encounter and $\beta \equiv v_0/c$. The burst duration can be of the order of a few milliseconds to several hours (depending on the PBH masses and encounter parameters) and bright enough to be detected at distances up to several Gpc. The maximum strain amplitude of the GW burst is given by

$$h_{\text{c}}^{\text{max}} = 3.24 \times 10^{-23} \frac{R_{\text{S}}(\text{km}) q \beta^2 g_{\text{max}}}{d_{\text{L}}(\text{Gpc}) (1+q)^2} \qquad (156)$$
<figure><img src="image_38.png" /><figcaption>Fig. 34 Example of SGWB expected for PBH binaries formed by tidal capture in clusters, with rates given by Eq. (153) and different virial velocities (2 km/s in green, 20 km/s in blue and 200 km/s in red) and a peaked PBH mass distribution around $30 M_{\odot}$. Image reproduced with permission from Clesse and García-Bellido (2017a)</figcaption></figure>

$$P_{\text{max}}(e) = 5.96 \times 10^{26} L_{\odot} \frac{q^2 \beta^{10}}{(1+q)^4} \frac{(e+1)}{(e-1)^5}, \quad (157)$$

where $L_{\odot}$ is the solar luminosity, $R_S$ is the Schwarzschild radius, $q \equiv m_1/m_2$ is the black hole mass ratio, $m_1 = qm_2 \ge m_2$ and $g_{\text{max}} = 2\sqrt{18(e+1) + 5e^2/(e-1)}$.

GWs detectable by LISA could be generated by close encounters of an intermediate-mass BH and a MBH as expected at the centers of galaxies, as well as from encounters of two MBHs that could occur during galactic collisions at low redshift. In the first case, an intermediate-mass BH of mass $m_2 = 10^3 M_{\odot}$ and a MBH of mass $m_1 = 10^6 M_{\odot}$, with an impact parameter $b = 1$ AU and velocity $v_0 = 0.05 c$, gives an eccentricity parameter $e = 1.031$ with event duration of about 440 s. The maximum stress amplitude would be $h_c^{\text{max}} = 1.02 \times 10^{-19}$ at a distance $d_L = 1$ Gpc, with the peak at frequency $f_{\text{peak}} = 1.05$ mHz, well within the sensitivity band of LISA [17]. In a hyperbolic encounter between two MBHs of equal mass $m_1 = m_2 = 10^6 M_{\odot}$ with impact parameter $b = 10$ AU and relative velocity $v_0 = 0.015 c$, the eccentricity is low $e = 1.01$ and the stress amplitude is huge $h_c^{\text{max}} = 2.22 \times 10^{-17}$ at $f_{\text{peak}} = 1.51 \times 10^{-4}$Hz, again right in the middle of the LISA observational band. Such a close hyperbolic encounter event would be clearly detectable due to a distinctive “burst” waveform (Morrás et al. 2022) that is very different from the usual “chirp” signatures of BBH. Furthermore, even if a single close hyperbolic encounter event may be difficult to detect, their characteristics can be distinguished population-wise, given their expected rates.
## 10.5 Massive black hole seeds from primordial black holes at high redshift

It is challenging to explain how MBHs can exist in only partially reionized environments at redshifts $z \ge 7$ (Banados et al. 2018). The first populations of stars or the direct collapse of gas into BHs are two possible astrophysical mechanisms to generate the seeds of these MBHs. But even if one invokes super-Eddington accretion, it is very challenging for these seeds to reach sufficiently large masses to explain observations. PBHs are an alternative explanation to the existence of MBHs since they can provide seeds of intermediate-mass BHs at higher redshift than for the other astrophysical mechanisms (Duechting 2004; Kawasaki et al. 2012; Clesse and García-Bellido 2015; Bernal et al. 2018). The easiest way to distinguish PBH seeds from other candidates is therefore to observe IMBBHs at $z \ge 20$, before star formation.

The astrophysical range of LISA will allow for the observation of IMBBH mergers at redshifts $z > 20$ with a SNR larger than five, for equal-mass mergers and progenitor masses between $10^3 M_{\odot}$ and $10^6 M_{\odot}$, as shown in Fig. 35. The possible merger rates of PBHs for a broad mass function with the imprints of the thermal history, shown in Fig. 33, can be larger than $O(1)\text{yr}^{-1}$ for primordial IMBBHs that would be formed in PBH clusters at high redshift. The existence of these clusters is relevant since they would also form in the standard Press-Schechter theory. LISA observations will be complementary to those of Earth-based GW detectors, like CE and ET, which will probe mergers with lower masses, and to future PTA limits from SKA, which will probe eventual mergers of MBHs at similar redshifts.

<figure><img src="image_39.png" /><figcaption>**Fig. 35** Redshift range of LISA according to the analysis of Burke-Spolaor et al. (2019) for equal-mass BBH coalescences as a function of the total system mass, and comparison with the range of other detectors and pulsar timing arrays. The color scale represents the expected SNR emerged from the study. Image reproduced with permission from Burke-Spolaor et al. (2019)</figcaption></figure>
# 11 Tools/pipelines for the analysis of transient signal data in cosmology

Section coordinator: L. Lombriser. Contributors: T. Baker, E. Belgacem, G. Calcagni, S. Clesse, G. Congedo, J.M. Ezquiaga, J. Garcia-Bellido, D. Laghi, L. Lombriser, M. Maggiore, M. Mancarella, A. Raccanelli, M. Sakellariadou, N. Tamanini, M. Zumalacarregui.

This section presents a discussion of the tools available and in need of development for the analysis of the LISA transient signal data, including codes, pipelines, algorithms, and methodologies. While the highlighted tools are of specific relevance to the analysis of LISA data, their application is generally not restricted to that. Sect. 11.1 focuses on standard sirens tests of the cosmic expansion (see Sect. 2). Section 11.2 is devoted to the tools for GW lensing (Sect. 3). Tools for tests of modified gravity (Sect. 4) are presented in Sect. 11.3, and finally Sect. 11.4 discusses the tools for PBHs (see Sect. 10).

## 11.1 Tools for testing the cosmic expansion with standard sirens

As presented in Sect. 2, GW observations are natural cosmological probes because their amplitude directly traces the luminosity distance as predicted by GR. When additional redshift information is obtained, GW events become standard sirens. Depending on whether the redshift information comes directly from an EM counterpart or from statistical inference, we will refer to them as bright and dark sirens respectively. In both cases, modern GW data analysis techniques are based on Bayesian inference. In the following we describe them independently.

### 11.1.1 Bright sirens

For bright sirens, the single-event posterior for a set of parameters $\mathcal{H}$ given GW data $D_{\text{GW}}$ and EM-wave data $D_{\text{EM}}$ can be written as (Chen 2020)

$$p(\mathcal{H}|D_{\text{GW}}, D_{\text{EM}}) = p(\mathcal{H}) \frac{\int p(D_{\text{GW}}|\vec{\Theta})p(D_{\text{EM}}|\vec{\Theta})p_{\text{pop}}(\vec{\Theta}|\mathcal{H})d\vec{\Theta}}{\int p_{\text{det}}(\vec{\Theta})p_{\text{pop}}(\vec{\Theta}|\mathcal{H})d\vec{\Theta}}, \quad (158)$$

where $p(\mathcal{H})$ is the prior probability density on $\mathcal{H}$, and $\vec{\Theta}$ represents all the binary parameters. The set of parameters $\mathcal{H}$ can be expressed in full generality as $\mathcal{H} = \{\mathcal{H}_{\text{cosmo}}, \mathcal{H}_{\text{pop}}\}$, where $\mathcal{H}_{\text{cosmo}} = \{H_0, \Omega_m, \Omega_k, \Omega_\Lambda\}$ are the cosmological parameters and $\mathcal{H}_{\text{pop}}$ represent a set of parameters used to model the population such as the mass and spin distribution or the merger rate evolution with redshift. The population prior $p_{\text{pop}}(\vec{\Theta}|\mathcal{H})$ is the probability density of binaries with parameters $\vec{\Theta}$ under our assumption of the merger rate density. The detection probability is given by
$$p_{\det}(\vec{\Theta}) = \int_{\mathcal{D}_{\text{GW}} > \text{GW}_{\text{thr}}} p(\mathcal{D}_{\text{GW}}|\vec{\Theta})p(\mathcal{D}_{\text{EM}}|\vec{\Theta})d\mathcal{D}_{\text{GW}}d\mathcal{D}_{\text{EM}}, \qquad (159) \\ \mathcal{D}_{\text{EM}} > \text{EM}_{\text{thr}}$$

in which the integration is only carried out over the data above the GW and EM detection thresholds, GWthr and EMthr, respectively.

GW and EM-wave observations can both introduce systematic uncertainties to the standard siren measurements. The measurements of the host redshifts can suffer from peculiar motions of the hosts (Howlett and Davis 2020; Mukherjee et al. 2021c; Nicolaou et al. 2020). The measurements of the binary luminosity distances are affected by the detector calibration uncertainty (Karki et al. 2016; Sun et al. 2020), lensing (Holz and Linder 2005; Hirata et al. 2010), and the accuracy of the waveforms (Abbott et al. 2019b). In addition, as shown in Eq. (158), our understanding of the possible observational selection effects (Chen 2020) as well as the astrophysical rate evolution (see e.g. Fig. 12 of Finke et al. (2021) for an application to LIGO-Virgo data) are critical to the accuracy of the standard siren as well.

### 11.1.2 Dark sirens

In the absence of an EM counterpart, one may extract redshift information by putting a prior on potential hosts from a galaxy catalogue.<fnref n="32" /> In this case, it is natural to consider the catalogue as a prior and the most convenient choice is to make a change of variables $d_L \to z(d_L, \mathcal{H})$ in Eq. (158), and use the redshift $z$ instead of the luminosity distance $d_L \in \vec{\Theta}$.<fnref n="33" /> Denoting $\vec{\Theta} = \{d_L, \hat{\Omega}, \vec{\Theta}'\}$, Eq. (158) reads

$$p(\mathcal{H}|\mathcal{D}_{\text{GW}}) = \frac{p(\mathcal{H})}{\alpha(\mathcal{H})} \int dz d\hat{\Omega} d\vec{\Theta}' p(\mathcal{D}_{\text{GW}}|d_L(z, \mathcal{H}), \hat{\Omega}, \vec{\Theta}') p_{\text{pop}}(z, \hat{\Omega}, \vec{\Theta}'|\mathcal{H}), \quad (160)$$

where $\alpha(\mathcal{H}) \equiv \int p_{\det}(\vec{\Theta})p_{\text{pop}}(\vec{\Theta}|\mathcal{H})d\vec{\Theta}$ and we set $p(\mathcal{D}_{\text{EM}}|\vec{\Theta}) = 1$, being in the absence of a counterpart. Note that the GW likelihood now carries an explicit dependence on the parameters $\mathcal{H}$. The population prior on redshift and position is expressed as

$$p_{\text{pop}}(z, \hat{\Omega}|\mathcal{H}) \propto \sum_i w_i \frac{dV_c}{dz d\hat{\Omega}}(z, \mathcal{H}) p(z_i, \hat{\Omega}_i|z, \hat{\Omega}). \qquad (161)$$

In the above equation, $p(z_i, \hat{\Omega}_i|z, \hat{\Omega})$ is the likelihood of the redshift and position of the i-th galaxy in the catalogue, with observed values $z_i, \hat{\Omega}_i$. The angular position can be assumed to be known with very good precision—since the GW likelihood does not vary significantly over the scale of the galaxy location—and the corresponding

<footnote n="32">³² This method assumes that the true source is not of primordial origin.</footnote>

<footnote n="33">³³ Note that the GW likelihood $p(\mathcal{D}_{\text{GW}}|\vec{\Theta})$ is unaffected by this transformation since it is a PDF with respect to the data $\mathcal{D}_{\text{GW}}$.</footnote>
distribution is just a delta function about the observed values; the redshift-dependent part of the likelihood can be instead approximated by a Gaussian centered around the observed redshift $z_i$ and with standard deviation given by the redshift uncertainty.<fnref n="34" /> The term $dV_c/dz d\Omega$ is the differential of the comoving volume element $V_c$, which is a uniform-in-volume prior on the redshift $z$ (Chen et al. 2018; Soares-Santos et al. 2019; Palmese et al. 2020). A more refined choice is to assume a prior that follows the expected evolution with redshift of the rate of GW sources (Gray et al. 2020; Finke et al. 2021). However, since this is not known a priori, this choice also adds more parameters to the analysis. Finally, the factors $w_i$ in Eq. (161) represent the possibility of assigning different weights to the galaxies. Common choices are the B-band luminosity, which traces the star formation rate, or K-band luminosity that traces the total stellar mass.

One has then to deal with the fact that galaxy catalogues are in general not complete. This requires one to include information on the missing galaxies, and the prior in Eq. (161) must be supplemented by a suitable “completion” term (Chen et al. 2018; Fishbach et al. 2019; Gray et al. 2020; Finke et al. 2021), which requires two pieces of information. The first is the knowledge on the “completeness” of the catalogue, that is usually computed by comparing the luminosity distribution in the catalogue to the one given by a Schechter function (Chen et al. 2018; Fishbach et al. 2019; Finke et al. 2021).<fnref n="35" /> The second is to specify a “completion” procedure, i.e. to determine how the missing galaxies are spatially distributed. The simplest options can be to distribute them uniformly (Chen et al. 2018; Fishbach et al. 2019), or to assume that the catalogue traces well the actual structures, hence assigning higher weight to the galaxies that are in low completeness regions (Finke et al. 2021), or even assuming that completion is not needed (Soares-Santos et al. 2019; Palmese et al. 2020; Laghi et al. 2021). Eventually, the best choice depends on the galaxy survey, its sky coverage and on the completeness of the catalogue. A more refined procedure is to use the uniform, uninformative distribution in regions of low completeness, while switching to the second option where the galaxy survey traces the actual structures fairly well (Finke et al. 2021). This aspect will become more important as the precision on the localisation increases but will remain a limiting factor until catalogues with very large completeness are available. Since the completeness drops drastically with increasing redshift, this can be limiting for LISA sources at high redshift when using catalogues that have a large sky coverage but limited redshift range. Another option is to use deeper catalogues with limited sky coverage, but this requires to take into account the variation of the completeness with the angular position in the sky (Finke et al. 2021). In any case, the availability of complete galaxy catalogues, accurate determination of the redshifts, and small GW localisation regions will be crucial in order for the statistical method to give competitive constraints on cosmological parameters.

<footnote n="34">³⁴ However, the true PDF should be used when available, and this choice can have a non-negligible impact on the result (Palmese et al. 2020).</footnote>

<footnote n="35">³⁵ A slightly different, though comparable, approach is to include the possibility that the true host is not in the catalogue in the formalism (Gray et al. 2020).</footnote>
Finally, a good understanding and an accurate computation of the selection bias is required (Chen 2020; Finke et al. 2021). This is related to other relevant systematic uncertainties that deserve further investigation, related to our limited knowledge of the population properties of BBHs. In particular, the population prior should encode information on the mass distribution of the system’s component and on the rate of BBH formation as a function of redshift. At the moment, very mild constraints exist on both in the case of SOBBH (Fishbach et al. 2018; Abbott et al. 2021e), but their impact—especially of the evolution of the BBH merger rate—can be relevant (Finke et al. 2021).

A public Python code implementing the completion procedures described above, as well as calculation of selection effects, is available at https://github.com/CosmoStatGW/DarkSirensStat. The code has been applied to the latest data release from the LIGO-Virgo collaboration for constraining $H_0$ and modified GW propagation and supports the combination with the GLADE, DES and GWENS galaxy catalogues, but it could easily be adapted for LISA and other galaxy catalogues.

## 11.2 Tools for gravitational-wave lensing

GW lensing has attracted considerable attention recently, leading to different groups developing tools to compute amplification factors and lensed waveforms. In order to maximise the science yield of gravitational lensing with LISA, it is important to develop codes to compute amplification factors for general lenses and including wave effects.

LensingGW is a publicly available tool for GW lensing in geometric optics. It is a Python package designed to compute amplification factors for general lenses, adapted to the needs of GW observations (Pagano et al. 2020). Currently, LensingGW operates in the geometric optics regime, using a selective sub-tiling to find microimages, their associated magnifications and time delays, cf. Eq. (14). The code is based on the lenstronomy package (Birrer and Amara 2018), allowing it to draw from its library of algorithms and lens models.

It is necessary to develop wave optics codes able to solve Eq. (15). The diffraction integral is highly oscillatory, which makes numerical solutions challenging. Proposed methods to address this problem use the time delay function (Ulmer and Goodman 1995) and Fourier transforms to the frequency domain. Variants of these ideas have been implemented in Diego et al. (2019), Cheung et al. (2021) and others. These codes are not publicly available. An alternative approach employs Picard-Lefschetz theory (Feldbrugge et al. 2023). This is valid only for lensing potentials that are meromorphic functions, which may be used to approximate non-analytic lensing potentials. Finally, corrections to geometric optics rely on higher derivatives of the Fermat potential around images (Takahashi 2004) and could be easily implemented in a geometric optics formalism.

Predictions on GW propagation beyond GR are notoriously complicated. (Ezquiaga and Zumalacárregui 2020) develops a Python package to study GW birefringence. The code is able to (1) compute static spherical backgrounds in Horndeski theory; (2) compute the local propagation eigenstates and their speed as a function of position and direction; (3) solve the geodesics in the Born approximation
to compute time delays and deflection angles for each eigenstate; and (4) obtain the waveform, including birefringence. Future tools to study GW lensing should improve by (a) including higher-order corrections in the Wentzel-Kramers-Brillouin (WKB) expansion to compute the amplitude and diffraction effects (or optimally, arrive at a full wave optics expression); (b) include more general theories of gravity; and (c) work with more general backgrounds (for instance, performing GW ray tracing on modified gravity simulations). These codes could be interfaced with a statistical sampler to analyse real or mock data.

## 11.3 Tools for testing modified gravity

The sensitivity of LISA to modified GW propagation introduced by a frequency independent modification of the damping term was investigated in Belgacem et al. (2019c). A MCMC analysis was conducted on LISA mock catalogues of MBBHs with EM counterparts, where luminosity distances were evaluated according to a fiducial ΛCDM cosmology. As discussed in (Belgacem et al. 2018b, 2019c), in the context of late-time modifications of GR invoked for DE studies, the most relevant parameters affecting GW propagation are the high-redshift ratio of GW to EM luminosity distances $\Xi_0$, defined in Eq. (21), the DE EoS $w_0$, the Hubble parameter $H_0$, and the current matter energy density fraction $\Omega_m$. The exponent $n$ of the parametrisation in Eq. (21) plays a less important role and was fixed to a reference value. In order to reduce the degeneracies among those parameters and assess the potential of LISA for cosmological parameter inference, in addition to LISA mock catalogues, the analysis in Belgacem et al. (2019c) also included further currently available datasets such as CMB measurements from Planck 2015, Type Ia Supernovae (JLA catalog), and a collection of BAO measurements. The likelihood of the LISA mock data given the parameters {$H_0, \Omega_m, w_0, \Xi_0$} was assumed to follow

$$ \ln(L(H_0, \Omega_m, w_0, \Xi_0)) = -\frac{1}{2} \sum_{i=1}^{N_s} \frac{[d_L^{\text{GW}}(z_i; H_0, \Omega_m, w_0, \Xi_0) - d_i]^2}{\sigma_i^2}, \quad (162) $$

where an additive constant in the logarithm of the likelihood from normalization is omitted. Here $N_s$ is the number of mock sources in the catalogue, $d_L^{\text{GW}}(z_i; H_0, \Omega_m, w_0, \Xi_0)$ is the theoretical value of the GW luminosity distance for the $i$-th source (at redshift $z_i$), $\sigma_i$ is the error on luminosity distance taking also into account the error on redshift determination, and $d_i$ is the “measured” value of the luminosity distance of each event contained in the catalogue (obtained by scattering the fiducial ΛCDM prediction with a Gaussian distribution with variance $\sigma_i^2$). The MCMC code of Belgacem et al. (2019c) explores the cosmological parameter space, accepting or rejecting sampled points following a Metropolis-Hastings algorithm. The total likelihood used in the algorithm is obtained by multiplying the independent contributions from each of the four datasets employed, namely from CMB, Type Ia supernovae, BAO, and LISA. The cosmological evolution, both at the level of the background and perturbations, is computed using the CLASS Boltzmann code. The modified version of CLASS, implementing the GW luminosity distance parametrisation in Eq. (21), is available at https://github.com/enisbelgacem/class_public. The
MCMC code is built on MontePython with the inclusion of the LISA mock source catalogues and the likelihood from Eq. (162), which uses the modified CLASS version. The code is available at https://github.com/enisbelgacem/montepython_public.

In addition to the test of the model-independent parameters {$Ξ_0, n$}, similar analyses have been carried out using the modified luminosity distances of specific gravity models, e.g. Horndeski gravity (Lombriser and Taylor 2016; Baker and Harrison 2021) or $f(Q)$ gravity (Frusciante 2021). These make use of the modified Boltzmann codes hi_CLASS (Zumalacárregui et al. 2017; Bellini et al. 2020) and EFTCAMB (Hu et al. 2014; Raveri et al. 2014), respectively. A wrapper to interface hi_CLASS with the MCMC code Cosmosis (Zuntz et al. 2015) is available at https://github.com/itrharrison/hi_class.

Current tools are predominantly specialised to probing the effects of frequency-independent modifications of the GW propagation, although frequency-dependent effects have been studied in the context of GW oscillations in bigravity (Belgacem et al. 2019c), where the associated numerical tools have however not been made publicly available yet. The desire for UV-completion of modified gravity theories, whilst maintaining IR phenomenology, also motivates frequency-dependent modifications of the GW propagation equation (de Rham and Melville 2018). Such modifications could be detectable in the LISA band whilst being suppressed to irrelevant levels in the band of ground-based detectors (~10–1000 Hz). Their basic phenomenology corresponds to a frequency-dependent GW luminosity distance or a frequency-dependent GW propagation speed.

To study the detectability of these effects, one must first construct a motivated and manageable parameterisation for the frequency and redshift dependence of $c_T(f, z)$ and $\delta(f, z)$ that appear in Eq. (19). One then needs to solve for the amplitude and phase evolution of the waveform under this parameterisation, up to some appropriate post-Newtonian order. A Fisher forecast code, ideally calibrated against a smaller, more rigorous (but computationally expensive) MCMC forecast, can analyse the ability of LISA to constrain the frequency-dependent modifications to $c_T$ and $d_{GW}$. The parameterisation must be carefully chosen to minimise degeneracies between the MG and standard source parameters. The creation of these tools is in progress (Baker et al. 2022).

Finally, we note that sources of anisotropic stress in the universe act as source terms on the right-hand side of Eq. (19). Neutrinos are one such example, although their effects on GW propagation in the late-time universe are expected to be extremely small. However, their effects on GWs travelling through the early universe may be more significant. Analytical tools for solving GW evolution in the presence of neutrino anisotropic stress are presented in Wren et al. (2017) and could likely be adapted for a generic imperfect fluid with anisotropic stress.

## 11.4 Tools for PBHs

The recent developments in the field of PBHs have revealed a rather rich phenomenology, e.g. related to PBH formation with a broad variety of primordial power spectra and new classes of models relying on non-trivial modifications of the
statistical distribution of curvature fluctuations; related to the history of PBH clustering, accretion and mergings, etc. The computations of GW observables therefore need to become more accurate and have to take into account model dependencies and various astrophysical uncertainties. It becomes increasingly difficult to integrate all the recent development in new analyses. Therefore, there is need in the community for an advanced numerical tool that progressively integrates all these recent developments at various levels (theoretical models, PBH formation, clustering, accretion, merging rates, SGWBs) in a unified and modular manner. This is the main objective of an ongoing project of the LISA Cosmology Working Group that is developing the *PrimBholes* toolbox for the computation of model-dependent gravitational observables that are relevant for LISA and other GW experiments. The *PrimBholes* toolbox is currently under development and will be made publicly available soon along with both a code companion paper and a review paper on PBHs.

The different computations performed with *PrimBholes* and the possible options and effects that will be implemented for selection are the following:

1. For the primordial curvature fluctuations this will include (a) several phenomenological models of power spectra (power-law, log-normal, Gaussian, broken power-law, etc.); (b) Gaussian and non-Gaussian perturbations with a generic distribution function and some specific realisations; and (c) generic power spectra and/or distribution functions imported from a file provided by the user.

2. The SGWB from second-order perturbations will be calculated for general as well a number of specific shapes of the primordial power spectrum.

3. The computation of the PBH formation and density distribution $\beta(m_{\text{PBH}})$ at formation will be performed following (a) the standard formalism of Sect. 10.2.2 and (b) an advanced method and algorithm presented in Musco et al. (2021) accounting for non-linear effects and the shape of the primordial power spectrum. This will include (c) effects of the evolution of the EoS through the thermal history of the universe and (d) a *reversed* method to rescale the amplitude of primordial curvature fluctuations in order to get a fixed value of $f_{\text{PBH}}$.

4. To obtain the late-time PBH mass function $f(m_{\text{PBH}})$, *PrimBholes* will rely on the standard formalism and include additional effects from accretion and hierarchical mergers.

5. The distribution of merging rates $R(m_1, m_2)$ will be computed for (a) primordial binaries, accounting for the effects from the formation of early clusters due to the Poisson noise, from matter inhomogeneities and from nearby BHs, and (b) for tidal capture in clusters, accounting for different clustering histories and halo mass functions. This will include (c) the dependence on the redshift of the merging rates.

6. Finally, *PrimBholes* will include the spin distribution, based on (De Luca et al. 2019; 7) the rate of hyperbolic encounters, based on (García-Bellido and Nesseris 2018); and (8) the SGWB from primordial binaries, for the two aforementioned formation channels.
PrimBholes will also provide plotting modules that will facilitate the production and exportation of key figures. Examples are given in Figs. 31 and 33.

In a first step, the PrimBholes code will focus on the computation of GW observables without including a likelihood module for LISA. However, it will be made modular so that it can be used in combination with other codes such as the LISA analysis code for the SGWBs or the Botzmann code CLASS (Blas et al. 2011). Later, likelihood modules for the MontePython code (Brinckmann and Lesgourgues 2019) could be included for the computation of the various astrophysical and cosmological limits on the PBH abundance.

# 12 Tools/pipelines for the analysis of stochastic gravitational-wave background data

Section coordinators: N. Karnesis, M. Peloso, M. Pieroni. Contributors: N. Bartolo, G. Boileau, N. Christensen, C. Contaldi, V. Desjacques, R. Flauger, N. Karnesis, V. Mandic, S. Mataresse, M. Peloso, M. Pieroni, A. Renzini, A. Ricciardone, J. Romano, M. Sakellariadou, L. Sorbo, J. Torrado.

The SGWB to be measured by LISA consists mainly of the confusion noise created by the overlapping signals from unresolved astrophysical events, and possibly of a cosmological component due to one or more physical mechanisms taking place at different stages of the late or early universe (e.g. the mechanisms discussed in Secs. 8, 7 and 6). LISA will open a completely new window into the parameter space of each of these phenomena and has the potential to produce new discoveries about the physics of the early universe through the characterisation of this primordial SGWB component. A careful characterisation of the total background is also necessary to aid the detection and identification of individual events. In this section, we discuss the features of each contribution to the background signal, the instrumental noise, and the different tools we need to develop to separate and characterise all of them.

We start with Sect. 12.1 which summarises the main characteristics of the noise in LISA. We comment on existing noise models, the assumptions on which they are built, and on possible improvements. In Sect. 12.2 we discuss foregrounds (i.e. astrophysical sources of a SGWB) by describing the possible sources and the characteristics and main properties of the different signals. In Sect. 12.3 we then focus on frequency shape reconstruction methods i.e. methods which aim at recovering the spectral profile of an unknown SGWB. We proceed by discussing anisotropy reconstruction i.e. recovering the angular structure hidden in the signal. For this topic both theoretical and data oriented approaches are covered in Sect. 12.4 and in Sect. 12.5 respectively. Finally, in Sect. 12.6 we frame the different topics touched in this section in the big picture of the so called *global fit* problem. For all of these points the discussion covers both existing techniques and possible further developments.
## 12.1 Noise modelling

The instrumental noise in LISA is expected to be non-stationary, and at the same time we expect noise transients that cause the statistical properties of the noise to depart from Gaussianity. At a first approximation we can assume that the short noise transients (i.e. glitches) will be modelled, identified, and removed from the data streams (Robson and Cornish 2019). However, slow variations of the noise, like for example the slow decrease of the acceleration noise of the test-masses due to outgassing (Armano et al. 2018), are due to known effects that can be properly modelled and considered in the analysis. We can tackle those problems by adopting a piecewise analysis, where we can consider that each data segment can be assumed as stationary and Gaussian. In order to eliminate the dominant noise sources (i.e. fluctuations of the laser central frequency and displacements of the optical benches), LISA will employ time domain interferometry (TDI) techniques (Tinto and Armstrong 1999; Estabrook et al. 2000; Tinto et al. 2002; Armstrong et al. 2003; Tinto et al. 2004; Tinto and Dhurandhar 2021; Vallisneri et al. 2021). In the simplified scenario considered in most of the literature (Bartolo et al. 2019a; Caprini et al. 2019; Contaldi et al. 2020; Smith et al. 2019; Pieroni and Barausse 2020; Flauger et al. 2021; Orlando et al. 2021) the residual noise for each arm link has two main components: the “acceleration” and the “interferometric” noise components. The acceleration component is associated with the random force noise acting on the test masses inside each of the three satellites, due for example to local environmental disturbances, and it dominates the low frequency part of the LISA band. The interferometric noise is directly connected to the interferometry metrology system (IMS). It describes the random readout noise of the optical system (mostly due to shot noise), and dominates the spectrum at high frequencies. For the 1.5 TDI variables (Bayle et al. 2019; Babak et al. 2021), one can follow the recipe of Flauger et al. (2021) and construct a likelihood based on the { X, Y, Z } TDI channels. In principle, analogous procedures can be defined for any kind of TDI variable combination. In order to remove noise correlations between different channels, it is customary to introduce an alternative TDI basis, typically dubbed the *AET* basis, which diagonalises the noise matrix. The total noise spectra for both XYZ and AET basis are shown in the right panel of Fig. 36. On the other hand, the left panel of Fig. 36 shows the acceleration and IMS noise power spectra.

Naturally, searching for SGWB signals requires a sufficient knowledge of the noise spectral shape, and it is possible that a simplified model such as the two-parameter model described above might not suffice for the signal search. Any instrument calibration mismatch will impact the joint fit of the noise and signal parameters, which may jeopardise our abilities to measure the underlying SGWB spectrum. Therefore, a possible solution to this problem would be to incorporate the calibration of the instrument in the analysis, by allowing a flexible model of the noise spectral shape. The development of methods (and of corresponding tools) which do not explicitly depend on analytical parameterisation for the noise spectra, will be an interesting research line to be explored over the years before the launch of LISA.

Another important aspect for searching for SGWB signals is assessing the quality of the data before they reach the designated pipelines. The use of signal subtraction,
<figure><img src="image_40.png" /><figcaption>**Fig. 36** Left panel: IMS and acceleration noise power spectra expressed in the simplified scenario where all test-mass and laser noises are equal for all space-crafts. Right panel: LISA noise spectra in the XYZ and AET basis as given in e.g. in Flauger et al. (2021). Images reproduced with permission from Flauger et al. (2021), copyright by IOP/SISSA</figcaption></figure>

if implemented for LISA data, could result in imperfect subtraction and hence residuals, which in turn would affect the statistical properties of the noise (Thrane 2013; Romano and Cornish 2017; Smith and Thrane 2018; Ginat et al. 2020). The impact of these residuals on the search for SGWBs with LISA is yet to be investigated. A better framing of this problem in the broader picture of the global fit problem for LISA is presented in Sect. 12.6.

## 12.2 Expected astrophysical foregrounds

Due to the richness of the sky in the milliHertz band LISA will be sensitive to a wide variety of sources such as MBBHs with masses ~$10^4–10^7 M_\odot$, SOBBHs, EMRs, and GBs. Beyond the resolvable sources, measurements by LISA will also be affected by a huge number of unresolvable events which will sum up incoherently, leading to the generation of a SGWB (Farmer and Phinney 2003; Regimbau and Hughes 2009; Regimbau 2011; Lamberts et al. 2019). At least two SGWB components are guaranteed to be present in the LISA band: a contribution due to signals originating from compact GBs is expected to dominate the LISA band at low frequencies (up to ~$10^{-3}$ Hz), and a contribution from extragalactic BBH mergers is expected to be present at slightly larger frequencies (~$10^{-3} – 10^{-2}$ Hz). A plot of the impact of these signals as described in Flauger et al. (2021) on data measured in the self correlations of the *A* and *T* channels (denoted by *AA* and *TT*, respectively; we note that the self correlation *EE* is identical to *AA*) is shown in Fig. 37. While the latter is expected to be isotropic and stationary, the GB contribution is expected to be anisotropic (since the binaries mostly lie on the galactic plane). Moreover, as discussed in Adams and Cornish (2014), due to the yearly rotation of the satellite constellation, this signal is expected to present a yearly modulation. Both these features are thus expected to be, at least partially, present in the total SGWB measured by LISA.

The characteristics of the GBs, i.e. being almost monochromatic, and the majority of them being located in our galactic neighbourhood, will allow us to build reliable
<figure><img src="image_41.png" /><figcaption>**Fig. 37** Theoretical models for the LISA sensitivity in the self-correlation of the *A* and *T* TDI data channels (*AA* and *TT*, respectively) and for the SGWBs due to GBs and SOBBHs as specified in Flauger et al. (2021). Mock data for the sum of all these theoretical templates are shown for reference</figcaption></figure>

models for their residual foreground noise contribution (Crowder and Cornish 2007; Sachdev et al. 2020). The same applies to the case of the isotropic signal due to SOBBH events, where the spectral model of the residual confusion noise will be constructed based on the priors from ground observations (Abbott et al. 2021d)<fnref n="36" /> and the actual measurements with LISA.

There is a possibility to measure a foreground signal component due to EMRIs (Barack and Cutler 2005), but a more detailed study on this type of source is needed in order to make robust predictions on their expected level. A generic method to make a first-level characterisation of the stochastic signals originating from compact binaries populations is studied in Karnesis et al. (2021).

The classical central-limit theorem is violated when the single source variance is infinite, which is the case of e.g. astrophysical backgrounds produced by compact mergers. As a consequence, the convergence to a Gaussian distribution can be much slower depending on the event rates (which determines the average number of sources overlapping in the detector frequency band(s)) and the properties of the source distribution (Ginat et al. 2020). If the sources are a spatial Poisson process, the resulting (non-Gaussian) distribution of the observed strain can be predicted (both in the time and frequency domains) using the techniques outlined in Ginat et al. (2020). This approach can be extended to include the subtraction of bright sources, which (partly) gaussianises the signal (Timpano et al. 2006), and to furnish predictions for the distribution of the background of unresolved binary mergers. Quantifying deviations from Gaussianity is necessary to optimise searches for this confusion

<footnote n="36">³⁶ This can be done by generating catalogues of events with the estimated population properties as measured by Earth-based detectors (Dvorkin et al. 2016b; Babak et al. 2023). Further investigation is needed to assess the potential uncertainties that arise from extrapolating the population properties as measured by Earth-based detectors to predictions in the LISA band.</footnote>
background (Thrane 2013; Smith and Thrane 2018) and backgrounds of cosmological origin (Bartolo et al. 2019a), and for their shape reconstruction.

Models for the residual foreground contributions can be provided when attempting a frequency shape reconstruction (see Sect. 12.3), so that the reconstructed spectrum can be expected to consist of contributions from unknown, possibly cosmological sources. The energy spectrum of the foreground due to unresolved GBs can be modelled by a broken power law. Indeed, at sufficiently large frequencies the number of sources decreases, which produces a break in the $\alpha = 2/3$ power law behaviour (where $\alpha$ is the spectral index in $h^2\Omega_{GX} \propto f^\alpha$). On the other hand, the foreground due to SOBBHs and BNSs is expected to be a power law with slope $\alpha = 2/3$. Studies of SOBBH and BNS populations, such as (Périgois et al. 2021), predict $\Omega_{GW}(25 \text{ Hz}) \approx 4.97 \times 10^{-9} - 2.58 \times 10^{-8}$. Ground detectors (Abbott et al. 2021g) set an upper limit for this foreground at $\Omega_{2/3}(f_{\text{ref}} = 25 \text{ Hz}) \lesssim 5.8 \times 10^{-9}$. Clearly, these residual foreground models must account for the right amount of statistical uncertainty, or we risk overestimating the significance of a detection or of a spectral shape reconstruction. This includes, at descending levels in the Bayesian hierarchy, the possibility of different population models being used (mass or spin distribution, merger rate at different redshifts, etc), the uncertainty on the parameters of a given population when inferred from resolved events by LIGO/Virgo and LISA, and the possibility of residuals being left by the imperfect subtraction of these resolved events, either because of a sub-optimal SNR that leads to imperfect characterisation of the waveform to be subtracted, or to the very waveform not representing the actual events with perfect accuracy. Beyond these problems, as discussed in the previous paragraph, the statistical properties of the SGWBs of astrophysical origin (either before and after the removal of loud sources) may be dominated by strong non-Gaussianities which have to be appropriately modelled. A careful study is thus required to assess the impact of each of these sources of statistical uncertainty. A computational pipeline that, for each component of the residual confusion noise, can take a number of resolved events with their associated parameter uncertainties, and assuming a Bayesian-hierarchical population model (or family of them), can generate a probability distribution for each component of the foreground confusion noise needs to be developed. The final result would, of course, rely on theoretical and/or observational uncertainties on the statistical properties of the populations in the parameter regions where the binaries cannot be individually resolved. However, the integrated contribution from all these unresolved sources would substantially affect the foreground (e.g. the SOBBH binaries at redshift $z \approx 2$ impact the SOBBH LISA foreground (Babak et al. 2023)). Improvements in our understanding of the LISA astrophysical populations (Amaro-Seoane et al. 2023) will thus be crucial to strengthen this kind of pipeline. Preliminary work on similar pipelines has been developed for current- and next-generation ground-based detectors (Sachdev et al. 2020; Périgois et al. 2021; Abbott et al. 2021g; Callister et al. 2020), but a comprehensive study is needed in the framework of LISA.
## 12.3 Spectral shape reconstruction: existing tools and future developments

The richness of sources in the LISA frequency band requires the development of techniques and tools to disentangle the different contributions. These may include both a residual background of unresolved astrophysical sources (see Sect. 12.2), a stochastic cosmological component sourced by topological defects, phase transitions, or inflationary mechanisms, among others. Different cosmological sources are expected to produce SGWBs of characteristic frequency shapes (Bartolo et al. 2016b; Caprini et al. 2016, 2020; Auclair et al. 2020). The reconstruction of the spectral shape of the SGWBs, regardless of its origin, is expected to play a special role in the separation of its different components. The correct estimation of the SGWB and its associated uncertainty is also essential for the characterisation and subtraction of resolved events, for which this background plays the role of an additional noise component, together with the instrumental noise.

Most of the SGWB detectability studies and ground-based searches have so far only focused on power-law templates (Thrane and Romano 2013; Adams and Cornish 2014; Lentati et al. 2015; Arzoumanian et al. 2018; Abbott et al. 2019c) (and a few on more complicated, but fixed templates (Kuroyanagi et al. 2018; Barish et al. 2021)). In the context of preparations for the LISA mission, some recent works (Karnesis et al. 2020; Caprini et al. 2019; Pieroni and Barausse 2020; Flauger et al. 2021) have attempted an agnostic template-free frequency shape reconstruction. One such reconstruction can take into account more complicated scenarios where the overall signal can be the superposition of several unknown signals, and can be extremely useful both for the disentanglement of the GW sources and for the constraining power that they have on astrophysical and cosmological parameters of the underlying theories (see previous sections). An example of reconstruction obtained with the method of Caprini et al. (2019), Flauger et al. (2021) in the presence of the foreground due to GBs is shown in Fig. 38.

<figure><img src="image_42.png" /><figcaption>Fig. 38 Reconstruction of a broken power law signal in presence of GB foreground. Image reproduced with permission from Flauger et al. (2021)</figcaption></figure>
In Karnesis et al. (2020), Caprini et al. (2019), Flauger et al. (2021), under the assumption of stationary, Gaussian, and isotropic signals, TDI data is processed into a set of frequency bins whose amplitudes are jointly estimated in the presence of noise. In the state-of-the art iteration (Flauger et al. 2021), the number of bins is dynamically chosen using Bayesian criteria, and in each of the bins, together with the signal amplitude, a power-law tilt is also considered, leading to a more accurate reconstruction with fewer bins. On the other hand, the procedure employed in Pieroni and Barausse (2020) first expresses the signal as a linear combination of an arbitrarily dense basis of Gaussians and then selects the components associated with large eigenvalues of the Fisher information matrix. While technically different, this strategy follows the same philosophical approach set by Karnesis et al. (2020), Caprini et al. (2019), Flauger et al. (2021) and offers an interesting alternative to the methods developed therein. For example, in Karnesis et al. (2020), after proper assumptions about the Gaussianity of the instrument noise, one can construct a posterior distribution for the level of the excess of spectral power, caused by an unmodelled SGWB signal at each frequency of the analysis.

All of the methods discussed in this section are focused on spectral shape reconstruction at the power spectrum level only. While in the most general set up this could be performed by only assuming the noise to be Gaussian, the techniques described here all assume (at the likelihood level) approximate Gaussianity of the SGWB as well. As long as the SNR per frequency mode is small, this is a good approximation but should be revisited for a bright signal. As discussed in Sect. 12.2 the assumption of Gaussianity can be violated for SGWBs of astrophysical origin. On the other hand, for SGWBs of cosmological origin any intrinsic non-Gaussianity in the metric perturbations is expected to be washed out by propagation effects (Bartolo et al. 2019b, c; Margalit et al. 2020). Notice however that, as pointed out in Bartolo et al. (2020b), non-Gaussianities in the energy density of the SGWB could be non-zero. For a better discussion of propagation effects see Sect. 5.3.2. If the backgrounds are expected to be significantly non-Gaussian, spectral shape reconstruction methods based on the observed strain distributions will have to be modified accordingly.

It is also worth mentioning that all the methods discussed in this section are applied to time-averaged data. While it could be difficult to induce any intrinsic time modulation in SGWBs of either cosmological or astrophysical origin, the motion of the detector with respect to the SGWB source frame is expected to naturally induce a yearly modulation in the amplitude of SGWB. While for isotropic SGWBs the effect is expected to be very small (of order of the velocity of the detector in units of the speed of light), for intrinsically anisotropic SGWBs (as for unresolved GBs) this effect could become sizeable (Adams and Cornish 2014). In order to model these effects it would be necessary to modify the existing pipelines to keep track (both at the generation and at the data analysis level) of the measurement time associated with each data segment. Such a modification in the existing tools would simultaneously improve their suitability to be applied to real data and possibly their capability of disentangling the different components contributing to the SGWB.

The methods discussed in this section do not consider possible angular structure in the signal. While techniques for reconstruction of anisotropies are more extensively
discussed in the next section, it is worth mentioning that for most SGWB sources different angular structure is expected to be associated with different spectra (namely, with a different frequency-dependence). As a consequence, by simultaneously tackling the frequency shape and the angular reconstruction it could be possible to break degeneracies between different components contributing to the observed signal. This could potentially lead to a significant improvement of the existing component separation techniques. For this reason embedding this possible method in the present pipelines may be an interesting possibility for future developments.

We conclude this section by commenting on the impact of foregrounds on the minimal SGWB intensity that would yield a detectable signal and allow us to attempt reconstruction. While detection thresholds are clearly model dependent, i.e. different amplitudes and SNRs are required in order to make quantitative statements for different SGWB shapes, some studies (see Sect. 12.6 for details) have already demonstrated the possibility of performing component separation with LISA. This implies that SGWB reconstruction can typically be performed even in the presence of foregrounds. However, it is worth mentioning that a crucial ingredient for component separation is the absence of degeneracies between the foregrounds and other SGWB components. Otherwise, unless this degeneracy is broken by the prior knowledge on the foreground parameters, the detection of the SGWB originating from cosmological sources could be negatively affected.

## 12.4 Anisotropy reconstruction

The intensity of astrophysical and cosmological SGWBs are in general anisotropic. These anisotropies could allow one to differentiate between these two backgrounds. They can be generated both at the moment of the GW production and by the GW propagation in the perturbed universe. In the case of the cosmological SGWB the anisotropies contain information about the primordial generation mechanism and moreover they can be a new probe of the primordial non-Gaussianity of the large-scale cosmological perturbations (Bartolo et al. 2019a, 2020a, b). In the case of the astrophysical SGWB they contain information about the angular distribution of the sources and astrophysical properties, and they can be used as a tracer of large-scale structure; in both cases they also allow one to test the particle physics content of the universe (Valbusa Dall'Armi et al. 2021). The LIGO/Virgo collaboration has already produced upper limits on the SGWB anisotropy in the 20–500 Hz band (Abbott et al. 2019a). The poor angular resolution of second-generation ground-based interferometers will probably not allow for the detection of the anisotropies of the SGWB, however future third-generation ground-based experiments like ET and CE will be sensitive to such a signal especially if it is characterised by a large monopole amplitude. For future space missions, like LISA, DECIGO, and BBO, prospects for reconstruction and measurement of the anisotropies in the SGWB, both astrophysical and cosmological, have been explored in the literature (see e.g. Giampieri and Polnarev 1997; Allen and Ottewill 1997; Cornish 2001; Seto 2004; Kudoh and Taruya 2005; Kudoh et al. 2006). An updated analysis for LISA, using the most up-to-date specifications, is ongoing (Bartolo et al. 2022), and in the next sections we report some main results. While in this section we discuss anisotropy reconstruction
from a more theoretical point of view, in Sect. 12.5 we present data analysis techniques to tackle this problem.

### 12.4.1 LISA response function

The reconstruction of an isotropic SGWB has been studied for instance in Flauger et al. (2021), combining measurements in the A,E,T channels. This computation can be readily extended to an anisotropic SGWB. Details of this extension can be found in Bartolo et al. (2022). Here we simply report the main result.

For a Gaussian SGWB, the object of our study is the two point correlation function of the GW signal. From the standard decomposition

$$h_{ab}(\mathbf{x}, t) = \int_{-\infty}^{+\infty} df \int d\Omega_{\hat{k}} e^{2\pi i f(t-\hat{k}\cdot\mathbf{x})} \sum_{\lambda} \tilde{h}_{\lambda}(f, \hat{k}) e_{ab}^{\lambda}(\hat{k}), \quad (163)$$

where, for definiteness, $\lambda$ denotes the GW polarisation in the chiral basis. We are interested in the intensity of an unpolarised signal, for which

$$\langle \tilde{h}_{\lambda}(f, \hat{k}) \tilde{h}_{\lambda'}(f', \hat{k}') \rangle = \frac{1}{4\pi} \delta(f+f') \delta^{(2)}(\hat{k}-\hat{k}') \delta_{\lambda,-\lambda'} I(f, \hat{k}). \quad (164)$$

We further decompose

$$I(f, \hat{k}) = \sum_{\ell m} \tilde{I}_{\ell m}(f) \tilde{Y}_{\ell m}(\hat{k}), \quad (165)$$

where $\tilde{Y}_{\ell m}(\hat{k}) = \sqrt{4\pi} Y_{\ell m}(\hat{k})$, are rescaled spherical harmonics, normalised so that $\tilde{Y}_{00} = 1$. In this way the isotropic case corresponds to $I(f, \hat{k}) = \tilde{I}_{00}(f)$.

The multipole coefficients $\tilde{I}_{\ell m}$ of the intensity are related to those of the fractional density by (Bartolo et al. 2022)

$$\tilde{I}_{\ell m}(f) = \frac{1}{\sqrt{4\pi}} \frac{3H_0^2}{4\pi^2} \frac{\Omega_{GW}(f)}{f^3} \delta_{GW,\ell m}, \quad (166)$$

where the latter are defined from

$$\frac{\omega_{GW}(f, \hat{k}) - \Omega_{GW}(f)}{\Omega_{GW}(f)} = \sum_{\ell m} \delta_{GW,\ell m}(f) Y_{\ell m}(\hat{k}), \quad (167)$$

with $\omega_{GW}(f, \hat{k})$ being the fractional energy density before angular integration,

$$\Omega_{GW}(f) = \int d^2\hat{k} \, \omega_{GW}(f, \hat{k}) / 4\pi. \quad (168)$$

In Eq. (166), $H_0$ is the present value of the Hubble rate.
For definiteness, we consider a statistically isotropic signal, so that the unbiased estimators $\hat{\delta}_{GW,\ell m}$ for coefficients of the decomposition $\delta_{GW,\ell m}$, have expectation value$^{37}$

$$ \langle \delta_{GW,\ell m} \delta_{GW,\ell' m'}^* \rangle = C_{\ell}^{GW} \delta_{\ell\ell'} \delta_{mm'} \quad (169) $$

The expectation value of the SNR can be then expressed as the sum $\langle \text{SNR} \rangle = \sqrt{\sum_{\ell} \langle \text{SNR}_{\ell}^2 \rangle}$ over multipoles. Each term is formally of the type

$$ \langle \text{SNR}_{\ell}^2 \rangle = T_{\text{obs}} \int_{0}^{\infty} df \left[ \frac{\sqrt{C_{\ell}^{\text{GW}} \Omega_{\text{GW}}(f) h^2}}{\Omega_{\text{GW},n}^{\ell}(f) h^2} \right]^2, \quad (170) $$

where $T_{\text{obs}}$ is the observation time, and $\Omega_{\text{GW},n}^{\ell}(f) h^2$ denotes the LISA sensitivity to the $\ell$-th multipole. The sensitivity to the first few multipoles are shown in Fig. 39. The mathematical form of the sensitivity, and its derivation, are given in Bartolo et al. (2022).

### 12.4.2 Cosmic dipole

One would expect that LISA sensitivity degrades as the multipole order increases. While this is generally true, as can be also seen from Fig. 39, due to the symmetries of the system, LISA turns out to be more sensitive to the quadrupole than to the dipole. Nevertheless, once the isotropic component of the SGWB is detected, a naturally expected value of the amplitude and the direction of the dipole will be determined by the velocity of the Earth with respect to the cosmic rest frame, providing a first target for the detection of SGWB anisotropies.

There are only two correlators sensitive to the dipole, namely the one between the A and the T channel, and the one between the E and the T channel. Remarkably, moreover, these AT and ET correlators vanish for the SGWB monopole. The SNR for the dipole induced by the Earth's motion with speed $\beta$ with respect to the cosmic rest frame is given, in the noise dominated regime, by the formula

$$ \text{SNR}_{A/ET}^2 = 6 \beta^2 \left( \frac{3 H_0^2 \Omega_{\text{GW}}}{4 \pi^2} \right)^2 \int_0^{T_{\text{obs}}} dt \int \frac{df}{f^6} \sum_{J=A,E} \frac{\left( \mathcal{R}_{\ell=1}^{JT} \hat{n}(t) \cdot \hat{\beta} \right)^2}{N_J N_T}, \quad (171) $$

where $N_J$ is the noise in mode J (see the right panel of Fig. 36) and $\mathcal{R}_{\ell=1}^{AT} = \mathcal{R}_{\ell=1}^{ET}$ measures the amplitude of the $\ell=1$ mode of the angular response (see Bartolo et al. 2022), evaluating to $\sim 10^{-2}$ for $f = 0.1$ Hz and to $\sim 10^{-3}$ for $f = 0.01$ Hz. Also, $\hat{\beta}$

<footnote>$^{37}$ Notice that the brackets in Eq. (164) denote an ensemble average. Since we have a single realisation of the observable universe, an estimator $\hat{C}_{\ell}^{GW}$ for $C_{\ell}^{GW}$ is typically built by averaging the measured value of $\delta_{GW,\ell m} \delta_{GW,\ell' m'}$ over the different m indices.</footnote>
<figure><img src="image_43.png" /><figcaption>Fig. 39 LISA Sensitivity (170) to various multipoles of the SGWB. Image reproduced with permission from Bartolo et al. (2022)</figcaption></figure>

is the direction of the dipole, whereas $\hat{n}(t)$ is the unit normal to the LISA configuration, whose direction changes over the course of the year.

Assuming a flat energy spectrum in the LISA band, the SNR takes the approximate value

$$ \text{SNR} \approx 0.5 \times \left( \frac{\beta}{10^{-3}} \right) \left( \frac{\Omega_{\text{GW}}}{6 \times 10^{-8}} \right) \sqrt{\frac{T_{\text{obs}}}{1 \text{ year}}} \quad (172) $$

where $\beta$ is normalised to its value with respect to the CMB rest frame, and $\Omega_{\text{GW}}$ is normalised to its upper bound from LIGO/Virgo (Abbott et al. 2019c). While Eq. (172) suggests a small SNR already for a SGWB that saturates the current constraint, it is important to note that the bound (Abbott et al. 2019c) holds at LIGO/Virgo frequencies, and there is no reason for the SGWB spectrum to be flat all the way to LISA frequencies.

Remarkably, due to the higher sensitivity to the quadrupole, LISA will have comparable sensitivity to the dipole and to the quadruple induced by the motion of the Solar System with respect to the cosmic rest frame, despite the fact that the quadrupole has an amplitude that is a factor of $\beta \approx 10^{-3}$ smaller than the dipole.
It is also worth pointing out that LISA will be much more sensitive to the dipole if the SGWB contains a chiral component. In this case, in fact, the AE correlator is not vanishing, leading to a SNR (Domcke et al. 2020b)

$$ \text{SNR} \approx 10^3 \delta\chi \left(\frac{\beta}{10^{-3}}\right) \left(\frac{\Omega_{\text{GW}}}{6 \times 10^{-8}}\right) \sqrt{\frac{T_{\text{obs}}}{1 \text{ year}}} \quad (173) $$

where $0 \le \delta\chi \le 1$ measures the degree of chirality of the SGWB, see Sect. 8.7.1.

## 12.5 Gravitational-wave map making with LISA

An available mapping algorithm tailored to the LISA detector is based on an optimal quadratic estimator (Contaldi et al. 2020). The mapper takes advantage of the time-dependent sky response of the LISA constellation to scan the sky over a long observation times. It reconstructs full-sky stochastic signals and their anisotropies. A prime example of anisotropic SGWB signal in the LISA band is the statistical signal from unresolved white dwarf GBs. This signal is expected to trace out the distribution of white dwarfs in the Milky Way. Beyond this, the astrophysical SGWB from extragalactic stellar mass compact binaries is also expected to have some degree of angular anisotropy, as do many of the cosmological backgrounds reviewed in this paper. More details on anisotropic SGWB sources potentially detectable by LISA may be found in Bartolo et al. (2022).

The sky response used assumes equal arms and heliocentric circular orbits, which are considered to be good approximations until full, time-dependent flight solutions become available after launch. These are the same assumptions made throughout Sect. 12.4.1. The frequency transfer function included in the response induces a frequency dependence in the sky modulation, as may be seen in Fig. 40. This drives the effective angular resolution as a function of frequency sensitivity and determines the resolution of final maps. The map-making can be based on a broad-band integration of frequencies or the integration can be split into narrower frequency bands (Contaldi et al. 2020).

The estimator is based on the assumption that the signal and noise components in the timestream data are Gaussian and maximises the standard Gaussian likelihood. This approach is similar to optimal map-making steps in CMB analysis (see e.g. Bond et al. 1998). In this case, the likelihood for the data is parametrised by the sky-signal covariance as function of direction. The maximum-likelihood solution for the

<figure><img src="image_44.png" /><figcaption>Fig. 40 Normalised auto-correlated response of TDI channel X at time $t = 0$ in the Solar System barycentre reference frame, at frequencies $f = 10^{-4}$ Hz, $f = 10^{-1}$ Hz, $f = 5 \times 10^{-1}$ Hz from left to right respectively. Estimates based on (Contaldi et al. 2020)</figcaption></figure>
covariance (signal intensity) is obtained using a quadratic, iterative estimator which effectively inverts the time-integrated projection of the sky-signal onto the data. A noise model is added to the total covariance in the likelihood. The algorithm is tailored to and applied in the pixel domain but it can also be applied in the spherical harmonic domain. The two choices differ in the use of regularisation methods that need to be applied to the ill-conditioned problem.

In practice, the data are segmented into short-duration frames, throughout which the sky response of the detector is considered to be constant. The time-frame length is a key element in the analysis: it sets the lower limit of frequency space probed, while also setting the maximum pixelisation scale to ensure a smooth transition on the sky, frame by frame. The covariance of the signal is initially estimated by subtracting the noise model from the data covariance, and is then noise-weighted over frequencies and projected onto the sky via the response operator to obtain a map. This operation is performed over all time-segments and subsequently averaged. The pixel-pixel Fisher matrix relative to the measurement is similarly calculated and inverted, using singular-value-decomposition techniques to regularise the inversion. The inverse Fisher matrix is then applied to the sky-and-noise weighted data to extract the optimal map solution. This is performed iteratively, until convergence requires ments are satisfied.

Figure 41 shows the $\ell$-mode transfer function $T_\ell$ for sky-map reconstruction presented in Contaldi et al. (2020) in the case of very high SNR injections. This shows the effective angular resolution of LISA at different frequency pivot points. This is calculated by combining sets of 50 simulations for each frequency window examined. Each simulation consists of the injection and reconstruction of an $\ell^2 C_\ell$-flat Gaussian realisation for the sky GW intensity. The procedure works under the assumption that the frequency and directional dependence can be separated:

<figure><img src="image_45.png" /><figcaption>**Fig. 41** Transfer functions $T_\ell$ for the average reconstructed $C_\ell$s obtained with different frequency cutoffs $f_{max}$, drawn from Contaldi et al. (2020). Each simulation set consists of 50 maps, each a different realisation of the same $C_\ell$ input (dashed grey baseline). There appears to be a clear one-to-one relation between the resolution $\ell_{max}$ of the instrument and the frequency cutoff</figcaption></figure>
$$ \Omega_{GW}(f, \mathbf{n}) = \Omega(f)\mathcal{P}(\mathbf{n}). \qquad (174) $$

For simplicity, the spectral dependence of the signal is taken to be a simple power law, with spectral index $\alpha = 3$. The transfer function is defined as $T_{\ell} = C_{\ell}^{\text{out}}/C_{\ell}^{\text{in}}$. In this extremely high signal scenario, the solution converges after a single iteration. As shown in Fig. 41, the higher angular modes are well preserved when allowing the reconstructor to integrate up to higher frequencies, where there is finer structure in the response pattern, whereas they are aliased into lower modes when integrating over lower frequencies only. In the best case scenario of a strong signal at high frequencies, the map-maker can recover anisotropies up to scales of $\ell_{\text{max}} \sim 15$. The reconstruction of a simplified model for the GB background is also tested; it is found that even though the signal has a relatively high SNR, it peaks above the noise around $10^{-3}$ Hz and hence may be recovered at best at $\ell_{\text{max}} \sim 5$.

The estimator currently assumes a fixed instrumental noise model, however it is possible to include independent noise estimates for each time segment. A possible extension of this method would fit for both signal and noise components simultaneously. The current noise model employed in the mapping pipeline is the one presented in Sect. 12.1. It is also possible to extend the estimator to include different spectral shapes for the signal, beyond a simple power law, or alternatively solve for maps in narrow frequency bands to attempt a model-independent directional search. Furthermore, it is important to note that the Fisher matrix is highly singular and requires heavy conditioning before inversion; this is a delicate process which may bias the outcome of the mapping, and further investigation is required to assess its impact. The implications of the inversion for certain high and low SNR stationary SGWB signals are discussed in Contaldi et al. (2020).

Going a step further, a recently developed Bayesian formalism (Banagiri et al. 2021) allows for simultaneous estimation of the frequency and directional content of a SGWB with LISA.<fnref n="38" /> One again starts with the assumption in Eq. (174), where a frequency dependence of choice could be inserted, for example a simple power law (as for example in Adams and Cornish 2010, 2014; Caprini et al. 2019; Flauger et al. 2021), while the angular dependence can be decomposed with respect to a basis on the sphere, such as spherical harmonics.

One can then divide LISA data into segments of duration $T_{\text{seg}}$ and compute the Fourier transforms of the LISA TDI channels in each segment (in either the X-Y-Z or A-E-T configuration). It is then straightforward to define the likelihood function

$$ \mathcal{L}(\tilde{d}|N_p, N_d, \alpha, \Omega_{\alpha}, \{b_{l,m}\}) = \prod_{t,f} \frac{1}{2\pi T_{\text{seg}} |C(t,f)|} \times \exp\left(-\frac{2\tilde{d}_{t,f}^* C(t,f)^{-1} \tilde{d}_{t,f}}{T_{\text{seg}}}\right), \qquad (175) $$

where $\tilde{d}_{t,f} = [\tilde{d}_X(t,f), \tilde{d}_Y(t,f), \tilde{d}_Z(t,f)]$ is the array of data in the Fourier domain for the three LISA channels (in the time segment $t$ and at frequency $f$), and $C(t, f)$ is the $3 \times 3$ covariance matrix for the three channels. The covariance matrix can be modelled to include both instrumental noise (here described by position and

<footnote n="38">³⁸ The code for this analysis pipeline can be found at https://github.com/sharanbngr/blip.</footnote>
acceleration noise parameters $N_p$ and $N_a$), as well as the astrophysical and cosmological contributions to the GWB. In this example, we use the amplitude $\Omega_\alpha$ and spectral index $\alpha$ to describe the frequency dependence of the astrophysical/cosmological contributions, and the $b_{l,m}$ coefficients of the decomposition onto spherical harmonics, to describe the directional dependence. Using this likelihood in Bayesian parameter estimation, with suitable choices for prior distributions on the free parameters, allows for the recovery of the free parameters. In particular, the recovery of simulated noise and a cosmological background up to $l=2$ is achievable (Banagiri et al. 2021). The same technique can be applied to study the galactic foreground due to white dwarf binaries, as shown in Fig. 42. Of course, additional model complexity could be added to this method: for example a more complex instrument noise model, or frequency dependence of $b_{l,m}$ coefficients, or polarisation dependence of the cosmological model.

## 12.6 Global fit

One of the main data analysis challenges that LISA will face is to disentangle the various sources that overlap in time and in frequency. To tackle this problem, the LISA Consortium will employ a *global fit* scheme (Teukolsky and Vallisneri 2019; Littenberg et al. 2020) for detecting, separating and eventually identifying/classifying many overlapping signals of astrophysical and hopefully early-universe origin. Different strategies could be employed to tackle this compelling task. A fully consistent Bayesian framework, ideally the most accurate approach to this problem, would perform a simultaneous fit of all resolvable sources and SGWBs. From a computational point of view this would be extremely difficult due to very high dimensionality of the parameter space that would need to be explored.

A possible approach (for a first step in this direction see Karnesis et al. 2021) would alternate between a step in which, given some current estimation of co-added background noise (instrumental noise, astrophysical confusion noise and a possible stochastic cosmological component), resolved events are detected, characterised and subtracted from the data streams, and a step in which the resulting “residual” data is used to update the knowledge on this background noise. In this step the residual is treated as a component which simultaneously models instrumental noise as well as the expected confusion noise inferred from the current catalogues of resolved sources (see Sect. 12.2), and a free-form SGWB spectrum (see Sect. 12.3). The advantage of this approach is that the dimensionality of the parameter inference problem at each of these steps is significantly lower (since technically we are performing conditional— or Gibbs—sampling at each step) and could be easily dealt with using standard techniques. Hereafter we discuss the state of the art on integrating SGWB detection, characterisation and component separation into a LISA *global fit* pipeline, and the further developments which are crucial for LISA to match its scientific goals regarding this problem.

The detection of an astrophysical SGWB with an accurate modelling, in the presence of the instrumental noise has been established in e.g. Christensen (1997), Adams and Cornish (2010), Meyers et al. (2020). In what concerns the simultaneous
<figure><img src="image_46.png" /><figcaption>**Fig. 42** Simulation (left) and recovery (right) of the white dwarf GB foreground energy density in GWs at 1 mHz, using one year of data in the Solar System barycentric frame. The bright spot in the map corresponds to the galactic central bulge. Image reproduced with permission from Banagiri et al. (2021)</figcaption></figure>

estimation of a cosmological component and one or more astrophysical confusion backgrounds, a number of approaches have been explored. We anticipate that their conclusions must be considered as preliminary since they all rely on some simplifying assumptions on the noise modelling, response functions and contamination from resolved-event misreconstructions, among others. Properly addressing these aspects is one of the main priorities in the forthcoming years.

Flauger et al. (2021) demonstrates the capability of LISA to separate between general astrophysical templates and a free-form cosmological component of sufficient intensity. In this work the three LISA TDI channels are fully exploited, using the *AET* basis (Hogan and Bender 2001; Adams and Cornish 2010), in particular the ability of the closed-path *T* channel to place strong, independent constraints on the instrumental noise, using the noise model of Sect. 12.1. The proposed pipeline determines an optimal binning for the reconstruction of the cosmological signal using maximum likelihood estimation. This binning is then used in a fully Bayesian pipeline implemented in the Cosmological sampling framework Cobaya (Torrado and Lewis 2021), to jointly estimate the parameters of the noise model, the amplitude of an astrophysical foreground template, and the free-form binned cosmological component. Using the Monte Carlo Nested Sampler PolyChord (Handley et al. 2015), it turns out that a sufficiently high-SNR broken-power-law signal can be efficiently reconstructed in the presence of either a power-law model for extragalactic binaries consistent with LIGO/Virgo observations (Abbott et al. 2019c), or a foreground of GBs modelled as in Robson et al. (2019). On the contrary, SGWBs with SNR $\lesssim 10$ such as flat signal with $\Omega_{GW} \lesssim 10^{-13}$ are likely to escape the LISA searches.

Another step in this direction is taken in Boileau et al. (2021a) and Boileau et al. (2021b) which progressively increase the complexity of the foreground. The analyses adopt a Bayesian strategy based on an Adaptive MCMC algorithm (Christensen and Meyer 1998; Cornish and Littenberg 2007) analysing the *A*, *E* and *T* channels given by the LISA model of Smith et al. (2019). The Adaptive MCMC results are independently confirmed by an analysis from the Fisher information matrix. Also with these analyses, it turns out that the LISA *T* channel helps to efficiently estimate the LISA noise parameters and to thus measure the SGWB in the *A* and *E* channels. Ultimately, it comes out that LISA can detect (power-law) cosmological backgrounds
in the presence of astrophysical foregrounds. In particular, given the expected LISA noise and the astrophysical foreground consistent with LIGO/Virgo observations, it is possible to observe a flat power-law SGWB with amplitude larger than $\Omega_{GW,\text{Cosmo}} \approx (1-10) \times 10^{-13}$ after 4 years of LISA observations (Boileau et al. 2021a). Boileau et al. (2021b) makes this lower bound more robust by adding the foreground due to GBs simulated with the binary catalogues of Lamberts et al. (2019). The study takes into account that the GB foreground measured by LISA has a yearly modulation (see Fig. 43). It implements a Bayesian analysis for each week in the year, and the corresponding variation of the LISA pattern antenna. Based on this analysis, the reconstruction of the amplitude of a flat SGWB has uncertainties below 50% when $\Omega_{GW,\text{Cosmo}} \gtrsim 8 \times 10^{-13}$ (see Fig. 44).

Similar efforts are carried out in the context of ground-based experiments. Martinovic et al. (2021) explores a method for the simultaneous estimation of astrophysical and cosmological SGWBs in the frequency band of current and future terrestrial interferometers. As an example, the detection of cosmological signals with GW fractional energy density $4.5 \times 10^{-13}$ at 25 Hz, for cosmic strings, and $2.2 \times 10^{-13}$, for a broken power-law model of an early universe phase transition, turns out to be possible in the presence of astrophysical confusion noise from compact binary coalescences, assuming a detector network containing CE and ET. In these approaches, reasonable levels of individual source subtraction are always assumed (Cutler and Harms 2006; Harms et al. 2008; Sharma and Harms 2020; Sachdev et al. 2020), given that future detector networks should be sensitive enough

<figure><img src="image_47.png" /><figcaption>**Fig. 43** Measurement of the orbital modulation of the white dwarf binary foreground. In grey: 1500 estimates of $\Omega_{\text{Mod},i} = \frac{4\pi^2}{3\hbar^2} \left(\frac{c}{2\pi L}\right)^2 A_i^2$ ($A_i$ amplitude of the characteristic strain). In red: 50 MCMC results with 8 parameters (2 parameters for BBHs, 4 parameters for white dwarf binaries, 2 parameters for the LISA noise). In green, fit to the 50 MCMC run results to estimate the modulation from the LISA antenna pattern amplitude at 3 mHz. Modulation model: $\Omega_{\text{Mod},i} = \Omega_{\text{DWD,LF}}^\mu (F_{+,i}^2 + F_{\times,i}^2)$. Image reproduced with permission from Boileau et al. (2021b)</figcaption></figure>
<figure><img src="image_48.png" /><figcaption>**Fig. 44** Cosmological amplitude uncertainty estimates from the Fisher Information study denoted by solid lines and from the MCMC by crosses. The upper horizontal dashed line represents the error level 50%. Image reproduced with permission from Boileau et al. (2021b)</figcaption></figure>

to accurately resolve a large fraction of the compact binary mergers in the visible
universe.

Biscoveanu et al. (2020) explores an alternative method to source subtraction. It proposes a Bayesian framework for simultaneously measuring individual compact binary mergers and a confusion foreground of these sources, together with a cosmological background. This method, in which the parameters of individual events are approximately marginalised-over, takes into account the non-Gaussian nature of the astrophysical signal (see Sect. 12.2), increasing the sensitivity to the Gaussian cosmological background. The capability of the method is demonstrated with a combination of an astrophysical foreground of merging BBHs and a cosmological background with a power-law spectrum.

The methods of both Martinovic et al. (2021) and Biscoveanu et al. (2020) could be adapted for the space-based detector LISA. In doing so, one should keep in mind that the corresponding spectral separation study for LISA would be more involved due to the nature of the TDI and the necessity to simultaneously estimate the LISA noise. A successful contribution separation will give us valuable information on the one hand about merger rates and population models and the astrophysics of exotic objects, and on the other about phenomenology in early universe models. This is indeed a research line that should be incentivated over the next years.

**Acknowledgements** This work is partly supported by: A.G. Leventis Foundation; Academy of Finland Grants 328958 and 345070; Alexander S. Onassis Foundation, Scholarship ID: FZO 059-1/2018-2019; Amaldi Research Center funded by the MIUR program "Dipartimento di Eccellenza" (CUP: B81118001170001); ASI Grants No. 2016-24-H.0 and No. 2016-24-H.1-2018; Atracción de Talento Grant 2019-T1/TIC-15784; Atracción de Talento contract no. 2019-T1/TIC-13177 granted by the Comunidad de Madrid; Ayuda 'Beatriz Galindo Senior' by the Spanish 'Ministerio de Universidades', Grant BG20/00228; Basque Government Grant (IT-979-16); Belgian Francqui Foundation; Centre national
d'Etudes spatiales; Ben Gurion University Kreitman Fellowship, and the Israel Academy of Sciences and Humanities (IASH) & Council for Higher Education (CHE) Excellence Fellowship Program for International Postdoctoral Researchers; Centro de Excelencia Severo Ochoa Program SEV-2016-0597; CERCA program of the Generalitat de Catalunya; Cluster of Excellence "Precision Physics, Fundamental Interactions, and Structure of Matter" (PRISMA+ EXC 2118/1); Comunidad de Madrid, Contrato de Atracción de Talento 2017-T1/TIC-5520; Czech Science Foundation GAČR, Grant No. 21-16583M; Delta ITP consortium; Department of Energy under Grant No. DE-SC0008541, DE-SC0009919 and DE-SC0019195; Deutsche Forschungsgemeinschaft (DFG), Project ID 438947057; Deutsche Forschungs-gemeinschaft under Germany's Excellence Strategy - EXC 2121 Quantum Universe - 390833306; European Structural and Investment Funds and the Czech Ministry of Education, Youth and Sports (Project CoGraDS - CZ.02.1.01/0.0/0.0/15 003/0000437); European Union's H2020 ERC Consolidator Grant "GRavity from Astrophysical to Microscopic Scales" (Grant No. GRAMS-815673); European Union's H2020 ERC, Starting Grant Agreement No. DarkGRA-757480; European Union's Horizon 2020 programme under the Marie Skłodowska-Curie Grant Agreement 860881 (ITN HIDDeN); European Union's Horizon 2020 Research and Innovation Programme Grant No. 796961, "AxiBAU" (K.S.); European Union's Horizon 2020 Research Council grant 724659 MassiveCosmo ERC-2016-COG; FCT through national funds (PTDC/FIS-PAR/31938/2017) and through project "BEYLA - BEYond LAmbda" with Ref. Number PTDC/FIS-AST/0054/2021; FEDER-Fundo Europeu de Desenvolvimento Regional through COMPETE2020 - Programa Operacional Competitividade e Internacionalização (POCI-01-0145-FEDER-031938) and research Grants UIDB/04434/2020 and UIDP/04434/2020; Fondation CFM pour la Recherche in France; Foundation for Education and European Culture in Greece; French ANR project MMUniverse (ANR-19-CE31-0020); FRIA Grant No.1.E.070.19F of the Belgian Fund for Research, F.R. S.-FNRS Fundação para a Ciência e a Tecnologia (FCT) through Contract No. DL 57/2016/CP1364/CT0001; Fundação para a Ciência e a Tecnologia (FCT) through Grants UIDB/04434/2020, UIDP/04434/2020, PTDC/FIS-OUT/29048/2017, CERN/FIS-PAR/0037/2019 and "CosmoTests - Cosmological tests of gravity theories beyond General Relativity" CEECIND/00017/2018; Generalitat Valenciana Grant PROMETEO/2021/083; Grant No. 758792, project GEODESI; Government of Canada through the Department of Innovation, Science and Economic Development and Province of Ontario through the Ministry of Colleges and Universities; Grants-in-Aid for JSPS Overseas Research Fellow (No. 201960698); I+D Grant PID2020-118159GB-C41 of the Spanish Ministry of Science and Innovation; INFN *iniziativa specifica* TEONGRAV; Israel Science Foundation (Grant No. 2562/20); Japan Society for the Promotion of Science (JSPS) KAKENHI Grant Nos. 20H01899 and 20H05853; IFT Centro de Excelencia Severo Ochoa Grant SEV-2; Kavli Foundation and its founder Fred Kavli; Minerva Foundation; Ministerio de Ciencia e Innovacion Grant PID2020-113644GB-I00; NASA Grant 80NSSC19K0318; NASA Hubble Fellowship grants No. HST-HF2-51452.001-A awarded by the Space Telescope Science Institute with NASA contract NAS5-26555; Netherlands Organisation for Science and Research (NWO) Grant Number 680-91-119; new faculty seed start-up grant of the Indian Institute of Science, Bangalore, the Core Research Grant CRG/2018/002200 of the Science and Engineering; NSF Grants PHY-1820675, PHY-2006645 and PHY-2011997; Polish National Science Center Grant 2018/31/D/ST2/02048; Polish National Agency for Academic Exchange within the Polish Returns Programme under Agreement PPN/PPO/2020/1/00013/U/00001; Pró-Reitoria de Pesquisa of Universidade Federal de Minas Gerais (UFMG) under Grant No. 28359; Ramón y Cajal Fellowship contract RYC-2017-23493; Research Project PGC2018-094773-B-C32 [MINECO-FEDER]; Research Project PGC2018-094773-B-C32 [MINECO-FEDER]; ROMFORSK Grant Project. No. 302640; Royal Society Grant URF/R1/180009 and ERC STG 949572: SHADE; Shota Rustaveli National Science Foundation (SRNSF) of Georgia (Grant FR/18-1462); Simons Foundation/SFARI 560536; SNSF Ambizione grant; SNSF professorship Grant (No. 170547); Spanish MINECO's "Centro de Excelencia Severo Ochoa" Programme Grants SEV-2016-0597 and PID2019-110058GB-C22; Spanish Ministry MCIU/AEI/FEDER Grant (PGC2018-094626-B-C21); Spanish Ministry of Science and Innovation (PID2020-115845GB-I00/AEI/10.13039/501100011033); Spanish Proyectos de I+D via Grant PGC2018-096646-A-I00; STFC Consolidated Grant ST/T000732/1; STFC Consolidated Grants ST/P000762/1 and ST/T000791/1; STFC Grant ST/S000550/1; STFC Grant ST/T000813/1; STFC Grants ST/P000762/1 and ST/T000791/1; STFC under the research Grant ST/P000258/1; Swiss National Science Foundation (SNSF), project *The Non-Gaussian Universe and Cosmological Symmetries*, Project Number: 200020-178787; Swiss National Science Foundation Professorship Grants No. 170547 and No. 191957; SwissMap National Center for Competence in Research; "The Dark Universe: A Synergic Multi-messenger Approach" Number 2017X7X85K under the MIUR program PRIN 2017; UK Space Agency; UKSA Flagship Project, Euclid.
**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

**Author Contribution** D. Bacon, G. Calcagni, S. Clesse, E. Dimastrogiovanni, J. M. Ezquiaga, J. García-Bellido, N. Karnesis, J. Kozaczuk, S. Kuroyanagi, M. Lewicki, L. Lombriser, M. Peloso, M. Pieroni, A. Raccanelli, L. Sousa, N. Tamanini, D. Vernieri and M. Zumalacárregui: Section coordinators. More details at the beginning of each section. R. Caldwell, G. Nardini, M. Peloso and N. Tamanini: Document coordinators.

## References

Abbott BP et al (2016a) GW150914: implications for the stochastic gravitational wave background from binary black holes. Phys Rev Lett 116(13):131102. https://doi.org/10.1103/PhysRevLett.116.131102. arXiv:1602.03847 [gr-qc]

Abbott BP et al (2016b) Observation of gravitational waves from a binary black hole merger. Phys Rev Lett 116(6):061102. https://doi.org/10.1103/PhysRevLett.116.061102. arXiv:1602.03837 [gr-qc]

Abbott BP et al (2017a) A gravitational-wave standard siren measurement of the Hubble constant. Nature 551(7678):85–88. https://doi.org/10.1038/nature24471. arXiv:1710.05835 [astro-ph.CO]

Abbott BP et al (2017b) Exploring the sensitivity of next generation gravitational wave detectors. Class Quantum Grav 34(4):044001. https://doi.org/10.1088/1361-6382/aa51f4. arXiv:1607.08697 [astro-ph.IM]

Abbott BP et al (2017c) Gravitational waves and gamma-rays from a binary neutron star merger: GW170817 and GRB 170817A. Astrophys J Lett 848(2):L13. https://doi.org/10.3847/2041-8213/aa920c. arXiv:1710.05834 [astro-ph.HE]

Abbott BP et al (2017d) GW170817: observation of gravitational waves from a binary neutron star inspiral. Phys Rev Lett 119(16):161101. https://doi.org/10.1103/PhysRevLett.119.161101. arXiv:1710.05832 [gr-qc]

Abbott BP et al (2017e) Multi-messenger observations of a binary neutron star merger. Astrophys J Lett 848(2):L12. https://doi.org/10.3847/2041-8213/aa91c9. arXiv:1710.05833 [astro-ph.HE]

Abbott BP et al (2017f) Upper limits on the stochastic gravitational-wave background from advanced LIGO’s first observing run. Phys Rev Lett 118(12):121101. https://doi.org/10.1103/PhysRevLett.118.121101, [Erratum: Phys. Rev. Lett. 119, 029901 (2017)]. arXiv:1612.02029 [gr-qc]

Abbott BP et al (2018a) Constraints on cosmic strings using data from the first Advanced LIGO observing run. Phys Rev D 97(10):102002. https://doi.org/10.1103/PhysRevD.97.102002. arXiv:1712.01168 [gr-qc]

Abbott BP et al (2018b) GW170817: implications for the stochastic gravitational-wave background from compact binary coalescences. Phys Rev Lett 120(9):091101. https://doi.org/10.1103/PhysRevLett.120.091101. arXiv:1710.05837 [gr-qc]

Abbott BP et al (2019a) Directional limits on persistent gravitational waves using data from Advanced LIGO’s first two observing runs. Phys Rev D 100(6):062001. https://doi.org/10.1103/PhysRevD.100.062001. arXiv:1903.08844 [gr-qc]

Abbott BP et al (2019b) Properties of the binary neutron star merger GW170817. Phys Rev X 9(1):011001. https://doi.org/10.1103/PhysRevX.9.011001. arXiv:1805.11579 [gr-qc]

Abbott BP et al (2019c) Search for the isotropic stochastic background using data from Advanced LIGO’s second observing run. Phys Rev D 100(6):061101. https://doi.org/10.1103/PhysRevD.100.061101. arXiv:1903.02886 [gr-qc]
Abbott BP et al (2020a) GW190425: observation of a compact binary coalescence with total mass ~3.4M⊙. Astrophys J Lett 892(1):L3. https://doi.org/10.3847/2041-8213/ab75f5. arXiv:2001.01761 [astro-ph.HE]

Abbott BP et al (2021a) A gravitational-wave measurement of the Hubble constant following the second observing run of advanced LIGO and Virgo. Astrophys J 909(2):218. https://doi.org/10.3847/1538-4357/abdcb7. arXiv:1908.06060 [astro-ph.CO]

Abbott R et al (2020b) GW190521: a binary black hole merger with a total mass of 150M⊙. Phys Rev Lett 125(10):101102. https://doi.org/10.1103/PhysRevLett.125.101102. arXiv:2009.01075 [gr-qc]

Abbott R et al (2020c) Properties and astrophysical implications of the 150 M⊙ binary black hole merger GW190521. Astrophys J Lett 900(1):L13. https://doi.org/10.3847/2041-8213/aba493. arXiv:2009.01190 [astro-ph.HE]

Abbott R et al (2021b) Constraints on cosmic strings using data from the third advanced LIGO-Virgo observing run. Phys Rev Lett 126(24):241102. https://doi.org/10.1103/PhysRevLett.126.241102. arXiv:2101.12248 [gr-qc]

Abbott R et al (2021c) Constraints on the cosmic expansion history from GWTC-3. arXiv e-prints arXiv:2111.03604 [astro-ph.CO]

Abbott R et al (2021d) GWTC-2: compact binary coalescences observed by LIGO and Virgo during the first half of the third observing run. Phys Rev X 11:021053. https://doi.org/10.1103/PhysRevX.11.021053. arXiv:2010.14527 [gr-qc]

Abbott R et al (2021e) Population properties of compact objects from the second LIGO-Virgo gravitational-wave transient catalog. Astrophys J Lett 913(1):L7. https://doi.org/10.3847/2041-8213/abe949. arXiv:2010.14533 [astro-ph.HE]

Abbott R et al (2021f) Tests of general relativity with binary black holes from the second LIGO-Virgo gravitational-wave transient catalog. Phys Rev D 103(12):122002. https://doi.org/10.1103/PhysRevD.103.122002. arXiv:2010.14529 [gr-qc]

Abbott R et al (2021g) Upper limits on the isotropic gravitational-wave background from Advanced LIGO and Advanced Virgo's third observing run. Phys Rev D 104(2):022004. https://doi.org/10.1103/PhysRevD.104.022004. arXiv:2101.12130 [gr-qc]

Abbott R et al (2023) Population of merging compact binaries inferred using gravitational waves through GWTC-3. Phys Rev X 13(1):011048. https://doi.org/10.1103/PhysRevX.13.011048. arXiv:2111.03634 [astro-ph.HE]

Abbott TMC et al (2018c) Dark energy survey year 1 results: a precise h0 estimate from DES Y1, BAO, and D/H data. Mon Not R Astron Soc 480(3):3879-3888. https://doi.org/10.1093/mnras/sty1939. arXiv:1711.00403 [astro-ph.CO]

Abell PA et al (2009) LSST science book, version 2.0. arXiv e-prints arXiv:0912.0201 [astro-ph.IM]

Accetta FS, Krauss LM (1989) The stochastic gravitational wave spectrum resulting from cosmic string evolution. Nucl Phys B 319:747-764. https://doi.org/10.1016/0550-3213(89)90628-7

Achucarro A, Gong JO, Hardeman S, Palma GA, Patil SP (2011) Mass hierarchies and non-decoupling in multi-scalar field dynamics. Phys Rev D 84:043502. https://doi.org/10.1103/PhysRevD.84.043502. arXiv:1005.3848 [hep-th]

Ackermann M et al (2018) The search for spatial extension in high-latitude sources detected by the *Fermi* large area telescope. Astrophys J Suppl 237(2):32. https://doi.org/10.3847/1538-4365/aacdf7. arXiv:1804.08035 [astro-ph.HE]

Acquaviva V, Bartolo N, Matarrese S, Riotto A (2003) Second order cosmological perturbations from inflation. Nucl Phys B 667:119-148. https://doi.org/10.1016/S0550-3213(03)00550-9. arXiv:astro-ph/0209156

Adachi Y, Maru N, Maru N (2020) Strong first order electroweak phase transition in Gauge-Higgs unification at finite temperature. Phys Rev D 101(3):036013. https://doi.org/10.1103/PhysRevD.101.036013. arXiv:1911.05963 [hep-ph]

Adams FC, Bond JR, Freese K, Frieman JA, Olinto AV (1993) Natural inflation: particle physics models, power law spectra for large scale structure, and constraints from COBE. Phys Rev D 47:426-455. https://doi.org/10.1103/PhysRevD.47.426

Adams MR, Cornish NJ (2010) Discriminating between a stochastic gravitational wave background and instrument noise. Phys Rev D 82:022002. https://doi.org/10.1103/PhysRevD.82.022002. arXiv:1002.1291 [gr-qc]

Adams MR, Cornish NJ (2014) Detecting a stochastic gravitational wave background in the presence of a galactic foreground and instrument noise. Phys Rev D 89(2):022001. https://doi.org/10.1103/PhysRevD.89.022001. arXiv:1307.4116 [gr-qc]
Addazi A, Marcianò A, Pasechnik R (2019) Probing trans-electroweak first order phase transitions from gravitational waves. Physics 1(1):92-102. https://doi.org/10.3390/physics1010010. arXiv:1811.09074 [hep-ph]

Addazi A et al (2022) Quantum gravity phenomenology at the dawn of the multi-messenger era—a review. Prog Part Nucl Phys 125:103948. https://doi.org/10.1016/j.ppnp.2022.103948. arXiv:2111.05659 [hep-ph]

Ade PAR et al (2014) Planck 2013 results. XXV. Searches for cosmic strings and other topological defects. Astron Astrophys 571:A25. https://doi.org/10.1051/0004-6361/201321621. arXiv:1303.5085 [astro-ph.CO]

Ade PAR et al (2018) BICEP2 / Keck Array x: constraints on primordial gravitational waves using Planck, WMAP, and New BICEP2/Keck observations through the 2015 season. Phys Rev Lett 121:221301. https://doi.org/10.1103/PhysRevLett.121.221301. arXiv:1810.05216 [astro-ph.CO]

Ade PAR et al (2021) Improved constraints on primordial gravitational waves using Planck, WMAP, and BICEP/Keck observations through the 2018 observing season. Phys Rev Lett 127(15):151301. https://doi.org/10.1103/PhysRevLett.127.151301. arXiv:2110.00483 [astro-ph.CO]

Adshead P, Wyman M (2012) Chromo-natural inflation: natural inflation on a steep potential with classical non-Abelian gauge fields. Phys Rev Lett 108:261302. https://doi.org/10.1103/PhysRevLett.108.261302. arXiv:1202.2366 [hep-th]

Adshead P, Giblin JT, Scully TR, Sfakianakis EI (2015) Gauge-preheating and the end of axion inflation. JCAP 12:034. https://doi.org/10.1088/1475-7516/2015/12/034. arXiv:1502.06506 [astro-ph.CO]

Adshead P, Giblin JT, Weiner ZJ (2018) Gravitational waves from gauge preheating. Phys Rev D 98 (4):043525. https://doi.org/10.1103/PhysRevD.98.043525. arXiv:1805.04550 [astro-ph.CO]

Adshead P, Giblin JT, Pieroni M, Weiner ZJ (2020a) Constraining axion inflation with gravitational waves across 29 decades in frequency. Phys Rev Lett 124(17):171301. https://doi.org/10.1103/PhysRevLett.124.171301. arXiv:1909.12843 [astro-ph.CO]

Adshead P, Giblin JT, Pieroni M, Weiner ZJ (2020b) Constraining axion inflation with gravitational waves from preheating. Phys Rev D 101(8):083534. https://doi.org/10.1103/PhysRevD.101.083534. arXiv:1909.12842 [astro-ph.CO]

Adshead P, Afshordi N, Dimastrogiovanni E, Fasiello M, Lim EA, Tasinato G (2021) Multimessenger cosmology: correlating cosmic microwave background and stochastic gravitational wave background measurements. Phys Rev D 103(2):023532. https://doi.org/10.1103/PhysRevD.103.023532. arXiv:2004.06619 [astro-ph.CO]

Agashe K, Du P, Ekhterachian M, Kumar S, Sundrum R (2020) Cosmological phase transition of spontaneous confinement. JHEP 05:086. https://doi.org/10.1007/JHEP05(2020)086. arXiv:1910.06238 [hep-ph]

Aggarwal K et al (2019) The NANograv 11-year data set: limits on gravitational waves from individual supermassive black hole binaries. Astrophys J 880:2. https://doi.org/10.3847/1538-4357/ab2236. arXiv:1812.11585 [astro-ph.GA]

Aghamousa A et al (2016) The DESI experiment part I: science, targeting, and survey design. arXiv e-prints arXiv:1611.00036 [astro-ph.IM]

Aghanim N et al (2020) Planck 2018 results. VI. Cosmological parameters. Astron Astrophys 641:A6. https://doi.org/10.1051/0004-6361/201833910, [Erratum: Astron. Astrophys. 652, C4 (2021)]. arXiv:1807.06209 [astro-ph.CO]

Agrawal A, Fujita T, Komatsu E (2018) Large tensor non-Gaussianity from axion-gauge field dynamics. Phys Rev D 97(10):103526. https://doi.org/10.1103/PhysRevD.97.103526. arXiv:1707.03023 [astro-ph.CO]

Ahonen J, Enqvist K (1998) Magnetic field generation in first order phase transition bubble collisions. Phys Rev D 57:664-673. https://doi.org/10.1103/PhysRevD.57.664. arXiv:hep-ph/9704334

Ahriche A, Hashino K, Kanemura S, Nasri S (2019) Gravitational waves from phase transitions in models with charged singlets. Phys Lett B 789:119-126. https://doi.org/10.1016/j.physletb.2018.12.013. arXiv:1809.09883 [hep-ph]

Aiola S et al (2020) The Atacama cosmology telescope: DR4 maps and cosmological parameters. JCAP 12:047. https://doi.org/10.1088/1475-7516/2020/12/047. arXiv:2007.07288 [astro-ph.CO]

Akhshik M (2015) Clustering fossils in solid inflation. JCAP 05:043. https://doi.org/10.1088/1475-7516/2015/05/043. arXiv:1409.3004 [astro-ph.CO]

Akhshik M, Emami R, Firouzjahi H, Wang Y (2014) Statistical anisotropies in gravitational waves in solid inflation. JCAP 09:012. https://doi.org/10.1088/1475-7516/2014/09/012. arXiv:1405.4179 [astro-ph.CO]
Akrami Y et al (2020) Planck 2018 results. X. Constraints on inflation. Astron Astrophys 641:A10. https://doi.org/10.1051/0004-6361/201833887. arXiv:1807.06211 [astro-ph.CO]

Alam MF et al (2021) The NANOGrav 12.5 yr data set: wideband timing of 47 millisecond pulsars. Astrophys J Suppl 252(1):5. https://doi.org/10.3847/1538-4365/abc6a1. arXiv:2005.06495 [astro-ph. HE]

Alam U, Sahni V, Saini TD, Starobinsky AA (2003) Exploring the expanding universe and dark energy using the Statefinder diagnostic. Mon Not R Astron Soc 344:1057. https://doi.org/10.1046/j.1365-8711.2003.06871.x. arXiv:astro-ph/0303009

Alanne T, Hugle T, Platscher M, Schmitz K (2020) A fresh look at the gravitational-wave signal from cosmological phase transitions. JHEP 03:004. https://doi.org/10.1007/JHEP03(2020)004. arXiv:1909.11356 [hep-ph]

Alba V, Maldacena J (2016) Primordial gravity wave background anisotropies. JHEP 03:115. https://doi.org/10.1007/JHEP03(2016)115. arXiv:1512.01531 [hep-th]

Albrecht A, Steinhardt PJ (1982) Cosmology for grand unified theories with radiatively induced symmetry breaking. Phys Rev Lett 48:1220-1223. https://doi.org/10.1103/PhysRevLett.48.1220

Alcock C et al (1997) The MACHO project LMC microlensing results from the first two years and the nature of the galactic dark halo. Astrophys J 486:697-726. https://doi.org/10.1086/304535. arXiv:astro-ph/9606165

Alexander S, Yunes N (2009) Chern-Simons modified general relativity. Phys Rep 480:1-55. https://doi.org/10.1016/j.physrep.2009.07.002. arXiv:0907.2562 [hep-th]

Alexander S, Finn LS, Yunes N (2008) A gravitational-wave probe of effective quantum gravity. Phys Rev D 78:066005. https://doi.org/10.1103/PhysRevD.78.066005. arXiv:0712.2542 [gr-qc]

Ali-Haïmoud Y (2018) Correlation function of high-threshold regions and application to the initial small-scale clustering of primordial black holes. Phys Rev Lett 121(8):081304. https://doi.org/10.1103/PhysRevLett.121.081304. arXiv:1805.05912 [astro-ph.CO]

Ali-Haïmoud Y, Kamionkowski M (2017) Cosmic microwave background limits on accreting primordial black holes. Phys Rev D 95(4):043534. https://doi.org/10.1103/PhysRevD.95.043534. arXiv:1612.05644 [astro-ph.CO]

Allahverdi R, Brandenberger R, Cyr-Racine FY, Mazumdar A (2010) Reheating in inflationary cosmology: theory and applications. Annu Rev Nucl Part Sci 60:27-51. https://doi.org/10.1146/annurev.nucl.012809.104511. arXiv:1001.2600 [hep-th]

Allahverdi R et al (2020) The first three seconds: a review of possible expansion histories of the early universe. arXiv e-prints https://doi.org/10.21105/astro.2006.16182. arXiv:2006.16182 [astro-ph.CO]

Allahyari A, Nunes RC, Mota DF (2022) No slip gravity in light of LISA standard sirens. Mon Not R Astron Soc 514(1):1274-1281. https://doi.org/10.1093/mnras/stac1445. arXiv:2110.07634 [astro-ph.CO]

Allen B (1996) The stochastic gravity wave background: sources and detection. In: Les Houches School of Physics: Astrophysical Sources of Gravitational Radiation, pp 373-417. arXiv:gr-qc/9604033

Allen B, Ottewill AC (1997) Detection of anisotropies in the gravitational wave stochastic background. Phys Rev D 56:545-563. https://doi.org/10.1103/PhysRevD.56.545. arXiv:gr-qc/9607068

Allys E, Peter P, Rodriguez Y (2016) Generalized SU(2) Proca theory. Phys Rev D 94(8):084041. https://doi.org/10.1103/PhysRevD.94.084041. arXiv:1609.05870 [hep-th]

Alonso D, Bellini E, Ferreira PG, Zumalacárregui M (2017) Observational future of cosmological scalar-tensor theories. Phys Rev D 95(6):063502. https://doi.org/10.1103/PhysRevD.95.063502. arXiv:1610.09290 [astro-ph.CO]

Alonso D, Contaldi CR, Cusin G, Ferreira PG, Renzini AI (2020a) Noise angular power spectrum of gravitational wave background experiments. Phys Rev D 101(12):124048. https://doi.org/10.1103/PhysRevD.101.124048. arXiv:2005.03001 [astro-ph.CO]

Alonso D, Cusin G, Ferreira PG, Pitrou C (2020b) Detecting the anisotropic astrophysical gravitational wave background in the presence of shot noise through cross-correlations. Phys Rev D 102 (2):023002. https://doi.org/10.1103/PhysRevD.102.023002. arXiv:2002.02888 [astro-ph.CO]

Alves A, Ghosh T, Guo HK, Sinha K (2018) Resonant Di-Higgs production at gravitational wave benchmarks: a collider study using machine learning. JHEP 12 (2018) 070. https://doi.org/10.1007/JHEP12 (2018) 070. arXiv:1808.08974 [hep-ph]

Alves A, Gonçalves D, Ghosh T, Guo HK, Sinha K (2020) Di-Higgs production in the 4b channel and gravitational wave complementarity. JHEP 03:053. https://doi.org/10.1007/JHEP03(2020)053. arXiv:1909.05268 [hep-ph]
Alves A, Gonçalves D, Ghosh T, Guo HK, Sinha K (2021) Di-Higgs blind spots in gravitational wave signals. Phys Lett B 818:136377. https://doi.org/10.1016/j.physletb.2021.136377. arXiv:2007.15654 [hep-ph]

Amaro-Seoane P, Santamaria L (2010) Detection of IMBHs with ground-based gravitational wave observatories: a biography of a binary of black holes, from birth to death. Astrophys J 722:1197–1206. https://doi.org/10.1088/0004-637X/722/2/1197. arXiv:0910.0254 [astro-ph.CO]

Amaro-Seoane P et al (2017) Laser Interferometer Space Antenna. arXiv e-prints arXiv:1702.00786 [astro-ph.IM]

Amaro-Seoane P et al (2023) Astrophysics with the Laser Interferometer Space Antenna. Living Rev Relativ 26:2. https://doi.org/10.1007/s41114-022-00041-y. arXiv:2203.06016 [gr-qc]

Amendola L (2000) Coupled quintessence. Phys Rev D 62:043511. https://doi.org/10.1103/PhysRevD.62.043511. arXiv:astro-ph/9908023

Amendola L, Ballesteros G, Pettorino V (2014) Effects of modified gravity on B-mode polarization. Phys Rev D 90:043009. https://doi.org/10.1103/PhysRevD.90.043009. arXiv:1405.7004 [astro-ph.CO]

Amendola L, Sawicki I, Kunz M, Saltas ID (2018a) Direct detection of gravitational waves can measure the time variation of the Planck mass. JCAP 08:030. https://doi.org/10.1088/1475-7516/2018/08/030. arXiv:1712.08623 [astro-ph.CO]

Amendola L et al (2018b) Cosmology and fundamental physics with the Euclid satellite. Living Rev Relativ 21:2. https://doi.org/10.1007/s41114-017-0010-3. arXiv:1606.00180 [astro-ph.CO]

Amin MA, Hertzberg MP, Kaiser DI, Karouby J (2014) Nonperturbative dynamics of reheating after inflation: a review. Int J Mod Phys D 24:1530003. https://doi.org/10.1142/S0218271815300037. arXiv:1410.3808 [hep-ph]

Amin MA, Braden J, Copeland EJ, Giblin JT, Solorio C, Weiner ZJ, Zhou SY (2018) Gravitational waves from asymmetric oscillon dynamics? Phys Rev D 98:024040. https://doi.org/10.1103/PhysRevD.98.024040. arXiv:1803.08047 [astro-ph.CO]

Ananda KN, Clarkson C, Wands D (2007) The cosmological gravitational wave background from primordial density perturbations. Phys Rev D 75:123518. https://doi.org/10.1103/PhysRevD.75.123518. arXiv:gr-qc/0612013

Anber MM, Sorbo L (2010) Naturally inflating on steep potentials through electromagnetic dissipation. Phys Rev D 81:043534. https://doi.org/10.1103/PhysRevD.81.043534. arXiv:0908.4089 [hep-th]

Anber MM, Sorbo L (2012) Non-Gaussianities and chiral gravitational waves in natural steep inflation. Phys Rev D 85:123537. https://doi.org/10.1103/PhysRevD.85.123537. arXiv:1203.5849 [astro-ph.CO]

Andersson L, Joudioux J, Oancea MA, Raj A (2021) Propagation of polarized gravitational waves. Phys Rev D 103(4):044053. https://doi.org/10.1103/PhysRevD.103.044053. arXiv:2012.08363 [gr-qc]

Ando K, Inomata K, Kawasaki M, Mukaida K, Yanagida TT (2018) Primordial black holes for the LIGO events in the axionlike curvaton model. Phys Rev D 97(12):123512. https://doi.org/10.1103/PhysRevD.97.123512. arXiv:1711.08956 [astro-ph.CO]

Ando S, Kusenko A (2010) Evidence for gamma-ray halos around active galactic nuclei and the first measurement of intergalactic magnetic fields. Astrophys J Lett 722:L39. https://doi.org/10.1088/2041-8205/722/1/L39. arXiv:1005.1924 [astro-ph.HE]

Andreev V et al (2018) Improved limit on the electric dipole moment of the electron. Nature 562 (7727):355–360. https://doi.org/10.1038/s41586-018-0599-8

Angelescu A, Huang P (2019) Multistep strongly first order phase transitions from new fermions at the TeV scale. Phys Rev D 99(5):055023. https://doi.org/10.1103/PhysRevD.99.055023. arXiv:1812.08293 [hep-ph]

Antusch S, Cefala F, Orani S (2017) Gravitational waves from oscillons after inflation. Phys Rev Lett 118 (1):011303. https://doi.org/10.1103/PhysRevLett.118.011303, [Erratum: Phys. Rev. Lett. 120, 219901 (2018)]. arXiv:1607.01314 [astro-ph.CO]

Antusch S, Cefala F, Orani S (2018) What can we learn from the stochastic gravitational wave background produced by oscillons? JCAP 03:032. https://doi.org/10.1088/1475-7516/2018/03/032. arXiv:1712.03231 [astro-ph.CO]

Apreda R, Maggiore M, Nicolis A, Riotto A (2002) Gravitational waves from electroweak phase transitions. Nucl Phys B 631:342–368. https://doi.org/10.1016/S0550-3213(02)00264-X. arXiv:gr-qc/0107033

Arai S, Nishizawa A (2018) Generalized framework for testing gravity with gravitational-wave propagation. II. Constraints on Horndeski theory. Phys Rev D 97(10):104038. https://doi.org/10.1103/PhysRevD.97.104038. arXiv:1711.03776 [gr-qc]