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

HAL
open science

# Cosmology with the Laser Interferometer Space Antenna

Pierre Auclair, David Bacon, Tessa Baker, Tiago Barreiro, Nicola Bartolo, Enis Belgacem, Nicola Bellomo, Ido Ben-Dayan, Daniele Bertacca, Marc Besancon, et al.

► To cite this version:

Pierre Auclair, David Bacon, Tessa Baker, Tiago Barreiro, Nicola Bartolo, et al. Cosmology with the Laser Interferometer Space Antenna. Living Rev.Rel., 2023, 26 (1), pp.5. 10.1007/s41114-023-00045-2. http://arxiv.org/abs/2302.03669

HAL Id: hal-03669190

https://hal.science/hal-03669190

Submitted on 15 Feb 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Cosmology with the Laser Interferometer Space Antenna

Pierre Auclair · David Bacon · Tessa Baker · Tiago Barreiro · Nicola Bartolo · Enis Belgacem et al. [full author details at the end of the article]

Received: 19 April 2022 / Accepted: 13 June 2023 / Published online: 28 August 2023

© The Author(s) 2023

## Abstract

The Laser Interferometer Space Antenna (LISA) has two scientific objectives of cosmological focus: to probe the expansion rate of the universe, and to understand stochastic gravitational-wave backgrounds and their implications for early universe and particle physics, from the MeV to the Planck scale. However, the range of potential cosmological applications of gravitational-wave observations extends well beyond these two objectives. This publication presents a summary of the state of the art in LISA cosmology, theory and methods, and identifies new opportunities to use gravitational-wave observations by LISA to probe the universe.

**Keywords** Laser Interferometer Space Antenna (LISA) · Cosmology

## Acronyms

BAO Baryon acoustic oscillations

BBH Binary black hole

BBN Big bang nucleosynthesis

BBO Big Bang Observer

BH Black hole

BNS Binary neutron star

BSM Beyond the standard model of particle physics

CE Cosmic Explorer

CP Charge parity

CMB Cosmic microwave background

CPL Chevalier-Polarski-Linder

DECIGO DECihertz Interferometer Gravitational wave Observatory

DE Dark energy

DES Dark Energy Survey

DESI Dark Energy Spectroscopic Instrument

EDM Electric dipole moment

<footnote>Extended author information available on the last page of the article</footnote>
<table>
  <tr>
    <td>E-ELT</td>
    <td>European-Extremely Large Telescope</td>
  </tr>
  <tr>
    <td>EFT</td>
    <td>Effective field theory</td>
  </tr>
  <tr>
    <td>EM</td>
    <td>Electromagnetic</td>
  </tr>
  <tr>
    <td>EMRI</td>
    <td>Extreme mass ratio inspiral</td>
  </tr>
  <tr>
    <td>EoS</td>
    <td>Equation of state</td>
  </tr>
  <tr>
    <td>ET</td>
    <td>Einstein Telescope</td>
  </tr>
  <tr>
    <td>EWPT</td>
    <td>Electroweak phase transition</td>
  </tr>
  <tr>
    <td>FLRW</td>
    <td>Friedmann-Lemaître-Robertson-Walker</td>
  </tr>
  <tr>
    <td>FOPT</td>
    <td>First-order phase transition</td>
  </tr>
  <tr>
    <td>GB</td>
    <td>Galactic binary</td>
  </tr>
  <tr>
    <td>GW</td>
    <td>Gravitational wave</td>
  </tr>
  <tr>
    <td>GR</td>
    <td>General relativity</td>
  </tr>
  <tr>
    <td>IMBBH</td>
    <td>Intermediate-mass binary black hole</td>
  </tr>
  <tr>
    <td>IMS</td>
    <td>Interferometry metrology system</td>
  </tr>
  <tr>
    <td>IR</td>
    <td>Infrared</td>
  </tr>
  <tr>
    <td>KAGRA</td>
    <td>Kamioka Gravitational-Wave Detector</td>
  </tr>
  <tr>
    <td>KiDS</td>
    <td>Kilo-Degree Survey</td>
  </tr>
  <tr>
    <td>ΛCDM</td>
    <td>Cosmological constant plus cold dark matter</td>
  </tr>
  <tr>
    <td>LIGO</td>
    <td>Laser Interferometer Gravitational-Wave Observatory</td>
  </tr>
  <tr>
    <td>LISA</td>
    <td>Laser Interferometer Space Antenna</td>
  </tr>
  <tr>
    <td>LSS</td>
    <td>Large scale structure</td>
  </tr>
  <tr>
    <td>MBBH</td>
    <td>Massive binary black hole</td>
  </tr>
  <tr>
    <td>MBH</td>
    <td>Massive black hole</td>
  </tr>
  <tr>
    <td>MCMC</td>
    <td>Markov Chain Monte Carlo</td>
  </tr>
  <tr>
    <td>MHD</td>
    <td>Magnetohydrodynamic</td>
  </tr>
  <tr>
    <td>NG</td>
    <td>Nambu Goto</td>
  </tr>
  <tr>
    <td>PBH</td>
    <td>Primordial black hole</td>
  </tr>
  <tr>
    <td>PISN</td>
    <td>Pair-instability supernova</td>
  </tr>
  <tr>
    <td>PLS</td>
    <td>Power law sensitivity</td>
  </tr>
  <tr>
    <td>ppE</td>
    <td>Parameterized post-Einsteinian</td>
  </tr>
  <tr>
    <td>PTA</td>
    <td>Pulsar timing array</td>
  </tr>
  <tr>
    <td>RD</td>
    <td>Radiation domination</td>
  </tr>
  <tr>
    <td>QCD</td>
    <td>Quantum chromodynamics</td>
  </tr>
  <tr>
    <td>SGWB</td>
    <td>Stochastic gravitational-wave background</td>
  </tr>
  <tr>
    <td>SKA</td>
    <td>Square Kilometer Array</td>
  </tr>
  <tr>
    <td>SM</td>
    <td>Standard model of particle physics</td>
  </tr>
  <tr>
    <td>SNR</td>
    <td>Signal-to-noise ratio</td>
  </tr>
  <tr>
    <td>SOBH</td>
    <td>Stellar-origin black hole</td>
  </tr>
  <tr>
    <td>SOBBH</td>
    <td>Stellar-origin binary black hole</td>
  </tr>
  <tr>
    <td>TDI</td>
    <td>Time domain interferometry</td>
  </tr>
  <tr>
    <td>UV</td>
    <td>Ultraviolet</td>
  </tr>
</table>
<table>
  <tbody>
    <tr>
      <td>5.3.1</td>
      <td>Frequency profile of the stochastic gravitational-wave background</td>
      <td>56</td>
    </tr>
    <tr>
      <td>5.3.2</td>
      <td>Anisotropies and propagation effects</td>
      <td>57</td>
    </tr>
    <tr>
      <td>5.3.3</td>
      <td>Polarisations</td>
      <td>59</td>
    </tr>
    <tr>
      <td>5.3.4</td>
      <td>Non-Gaussianity</td>
      <td>60</td>
    </tr>
    <tr>
      <td>6</td>
      <td>First-order phase transitions</td>
      <td>61</td>
    </tr>
    <tr>
      <td>6.1</td>
      <td>Introduction</td>
      <td>61</td>
    </tr>
    <tr>
      <td>6.2</td>
      <td>Determination of the relevant parameters</td>
      <td>62</td>
    </tr>
    <tr>
      <td>6.3</td>
      <td>GW sources</td>
      <td>64</td>
    </tr>
    <tr>
      <td>6.4</td>
      <td>Discovery prospects beyond the standard model</td>
      <td>67</td>
    </tr>
    <tr>
      <td>6.4.1</td>
      <td>First-order phase transitions in the LISA window</td>
      <td>69</td>
    </tr>
    <tr>
      <td>6.4.2</td>
      <td>Connection with baryogenesis</td>
      <td>71</td>
    </tr>
    <tr>
      <td>6.5</td>
      <td>LISA complementarity to other experimental tests</td>
      <td>73</td>
    </tr>
    <tr>
      <td>7</td>
      <td>Cosmic strings</td>
      <td>75</td>
    </tr>
    <tr>
      <td>7.1</td>
      <td>Introduction</td>
      <td>76</td>
    </tr>
    <tr>
      <td>7.2</td>
      <td>Network modelling</td>
      <td>78</td>
    </tr>
    <tr>
      <td>7.2.1</td>
      <td>Model 1</td>
      <td>78</td>
    </tr>
    <tr>
      <td>7.2.2</td>
      <td>Model 2</td>
      <td>80</td>
    </tr>
    <tr>
      <td>7.2.3</td>
      <td>Model 3</td>
      <td>81</td>
    </tr>
    <tr>
      <td>7.2.4</td>
      <td>Field theory</td>
      <td>82</td>
    </tr>
    <tr>
      <td>7.3</td>
      <td>Computation of the gravitational-wave spectrum from loops</td>
      <td>83</td>
    </tr>
    <tr>
      <td>7.3.1</td>
      <td>Method I</td>
      <td>84</td>
    </tr>
    <tr>
      <td>7.3.2</td>
      <td>Method II</td>
      <td>85</td>
    </tr>
    <tr>
      <td>7.3.3</td>
      <td>Templates for the stochastic gravitational-wave background generated by cosmic string loops</td>
      <td>87</td>
    </tr>
    <tr>
      <td>7.4</td>
      <td>LISA detection prospects</td>
      <td>89</td>
    </tr>
    <tr>
      <td>7.4.1</td>
      <td>Anisotropies in the stochastic gravitational-wave background</td>
      <td>91</td>
    </tr>
    <tr>
      <td>7.5</td>
      <td>Extended and alternative considerations</td>
      <td>93</td>
    </tr>
    <tr>
      <td>7.5.1</td>
      <td>Long strings</td>
      <td>94</td>
    </tr>
    <tr>
      <td>7.5.2</td>
      <td>Agnostic approach to loop size</td>
      <td>96</td>
    </tr>
    <tr>
      <td>7.5.3</td>
      <td>Global strings</td>
      <td>97</td>
    </tr>
    <tr>
      <td>7.5.4</td>
      <td>Cosmic superstrings</td>
      <td>99</td>
    </tr>
    <tr>
      <td>7.5.5</td>
      <td>Metastable strings</td>
      <td>100</td>
    </tr>
    <tr>
      <td>8</td>
      <td>Inflation</td>
      <td>102</td>
    </tr>
    <tr>
      <td>8.1</td>
      <td>Introduction</td>
      <td>102</td>
    </tr>
    <tr>
      <td>8.2</td>
      <td>Gravitational waves sourced by additional fields during inflation</td>
      <td>104</td>
    </tr>
    <tr>
      <td>8.2.1</td>
      <td>Axion-gauge field inflation</td>
      <td>105</td>
    </tr>
    <tr>
      <td>8.2.2</td>
      <td>Non-minimally coupled (spinning) fields</td>
      <td>107</td>
    </tr>
    <tr>
      <td>8.2.3</td>
      <td>Models that produce peaks in the power spectrum of scalar perturbations</td>
      <td>109</td>
    </tr>
    <tr>
      <td>8.3</td>
      <td>Small-scale primordial features</td>
      <td>111</td>
    </tr>
    <tr>
      <td>8.4</td>
      <td>Effective field theory of broken space-reparametrisations</td>
      <td>114</td>
    </tr>
    <tr>
      <td>8.5</td>
      <td>UV complete models of R<sub>α</sub>R<sup>2</sup> model</td>
      <td>116</td>
    </tr>
    <tr>
      <td>8.6</td>
      <td>Preheating</td>
      <td>118</td>
    </tr>
    <tr>
      <td>8.7</td>
      <td>Summary of distinctive gravitational-wave observables from inflation</td>
      <td>120</td>
    </tr>
    <tr>
      <td>8.7.1</td>
      <td>Chirality</td>
      <td>122</td>
    </tr>
    <tr>
      <td>8.7.2</td>
      <td>Frequency profile</td>
      <td>122</td>
    </tr>
    <tr>
      <td>8.7.3</td>
      <td>Anisotropies</td>
      <td>123</td>
    </tr>
    <tr>
      <td>8.7.4</td>
      <td>Cross-correlations</td>
      <td>125</td>
    </tr>
    <tr>
      <td>9</td>
      <td>Tests of non-standard pre-Big-Bang nucleosynthesis cosmology via the SGWB</td>
      <td>126</td>
    </tr>
    <tr>
      <td>9.1</td>
      <td>Introduction</td>
      <td>126</td>
    </tr>
    <tr>
      <td>9.2</td>
      <td>Non-standard expansion histories</td>
      <td>128</td>
    </tr>
    <tr>
      <td>9.2.1</td>
      <td>Inflation</td>
      <td>128</td>
    </tr>
    <tr>
      <td>9.2.2</td>
      <td>Phase transitions</td>
      <td>130</td>
    </tr>
    <tr>
      <td>9.2.3</td>
      <td>Cosmic strings</td>
      <td>132</td>
    </tr>
    <tr>
      <td>9.3</td>
      <td>SGWB in quantum gravity</td>
      <td>134</td>
    </tr>
    <tr>
      <td>10</td>
      <td>Primordial black holes</td>
      <td>138</td>
    </tr>
  </tbody>
</table>

<figure><img src="image_4.png" /></figure>
<table><tr><td>10.1 Introduction</td><td>140</td></tr><tr><td>10.2 Formation scenarios</td><td>142</td></tr><tr><td>10.2.1 Origin of curvature fluctuations</td><td>143</td></tr><tr><td>10.2.2 Standard theory of primordial black holes formation from Gaussian curvature fluctuations</td><td>145</td></tr><tr><td>10.2.3 Beyond Gaussianity</td><td>147</td></tr><tr><td>10.2.4 The threshold for primordial black holes and the impact of thermal history</td><td>147</td></tr><tr><td>10.3 Stochastic gravitational-wave background sourced at second order by curvature fluctuations</td><td>149</td></tr><tr><td>10.4 Resolved sources and stochastic gravitational-wave backgrounds from primordial black holes binaries and hyperbolic encounters</td><td>151</td></tr><tr><td>10.4.1 Binary formation channels and merger rates</td><td>153</td></tr><tr><td>10.4.2 Gravitational-wave chirps from primordial black holes binaries</td><td>155</td></tr><tr><td>10.4.3 Stochastic gravitational-wave background</td><td>156</td></tr><tr><td>10.4.4 Bursts from hyperbolic encounters</td><td>157</td></tr><tr><td>10.5 Massive black hole seeds from primordial black holes at high redshift</td><td>159</td></tr><tr><td>11 Tools/pipelines for the analysis of transient signal data in cosmology</td><td>160</td></tr><tr><td>11.1 Tools for testing the cosmic expansion with standard sirens</td><td>160</td></tr><tr><td>11.1.1 Bright sirens</td><td>160</td></tr><tr><td>11.1.2 Dark sirens</td><td>161</td></tr><tr><td>11.2 Tools for gravitational-wave lensing</td><td>163</td></tr><tr><td>11.3 Tools for testing modified gravity</td><td>164</td></tr><tr><td>11.4 Tools for PBHs</td><td>165</td></tr><tr><td>12 Tools/pipelines for the analysis of stochastic gravitational-wave background data</td><td>167</td></tr><tr><td>12.1 Noise modelling</td><td>168</td></tr><tr><td>12.2 Expected astrophysical foregrounds</td><td>169</td></tr><tr><td>12.3 Spectral shape reconstruction: existing tools and future developments</td><td>171</td></tr><tr><td>12.4 Anisotropy reconstruction</td><td>173</td></tr><tr><td>12.4.1 LISA response function</td><td>175</td></tr><tr><td>12.4.2 Cosmic dipole</td><td>176</td></tr><tr><td>12.5 Gravitational-wave map making with LISA</td><td>178</td></tr><tr><td>12.6 Global fit</td><td>180</td></tr><tr><td>References</td><td>186</td></tr></table>

# 1 Introduction

Contributors: R. Caldwell, G. Nardini.

The Laser Interferometer Space Antenna (LISA) (Amaro-Seoane et al. 2017) is a planned space-borne gravitational-wave (GW) detector that will open a new frontier on astrophysics and cosmology in the mHz frequency band. This European Space Agency-led mission includes participation by ESA member countries and significant contributions from NASA and space agencies of other countries. Phase A work is on track for mission adoption in mid 2020s, and is compatible with a launch in the mid 2030s.

LISA will consist of a trio of satellites, located at the vertices of an equilateral triangle, in an Earth-trailing heliocentric orbit. The 2.5-million km distances between the satellites will be monitored using precision laser interferometry to detect passing GWs. Here we consider a nominal mission of six years with a duty cycle of around 75%, although we understand that the instruments will be engineered to a specification that will enable a possible extension.
LISA will be sensitive to GWs from a wide array of sources (Amaro-Seoane et al. 2017). A primary target will be the inspiral and merger of massive binary black holes (MBBHs), ranging in masses $10^4–10^7 M_⊙$, at redshifts out to $z \sim 10$. A significant foreground signal will be the many galactic white dwarf binaries, each effectively a monotone source.

By the time LISA flies, the state of GW observation will have evolved. The extended Advanced Laser Interferometer Gravitational Wave Observatory (LIGO), Advanced Virgo, and Kamioka Gravitational Wave Detector (KAGRA) family of ground-based GW detectors will have begun implementing third-generation technology demonstration upgrades. The network of pulsar-timing radio telescopes will have grown to include the Square Kilometer Array (SKA). Yet, LISA will be different from its predecessors. The size of the detector will enable access to a completely fresh part of the GW spectrum, leading to observations of new astrophysical sources as well as a new window on primordial stochastic gravitational-wave backgrounds (SGWBs). Many sources will produce overlapping signals, owing to the improved sensitivity. Extracting individual sources and events, and discriminating from an unresolved hum, will be part of the challenge.

According to the mission proposal (Amaro-Seoane et al. 2017), LISA has two main scientific objectives of purely cosmological bearing. The first is to probe the expansion rate of the universe, with specific requirements to measure the dimensionless Hubble parameter by means of GW observations alone, and further to constrain cosmological parameters through joint GW and electromagnetic (EM) observations. The second such objective is to understand SGWBs and their implications for early universe and particle physics. This will entail the characterisation of the astrophysical SGWB, and subsequently a measurement or bound on the amplitude and spectral shape of a cosmological SGWB. There are further scientific imperatives to use LISA to explore the fundamental nature of gravity and to search for unforeseen sources with relevance for cosmology. There is a wealth of cosmological information that may be extracted from LISA observations.

We start with Sects. 2 and 3 on standard sirens and weak gravitational lensing; these are “sure bets” for LISA, based on our current understanding of source populations. These sections are directly related to LISA science objective SO6 “probe the rate of the expansion universe” (LISA Science Study Team 2018). They also identify new opportunities to derive cosmological information from GW astrophysical sources, in connection with LISA science objectives SO1, SO2, SO3 and SO4, which are devoted to understanding the galactic and extragalactic astrophysical source populations (LISA Science Study Team 2018). We follow with sections on more exploratory topics, which are potentially profound and revolutionary. Sect. 4 discusses the constraints on modified gravity theories that may be achieved through measurement of GW sources at cosmological distances. Results on this research subject are aligned with LISA science objectives SO5 “explore the fundamental nature of gravity and black holes” as well as the aforementioned SO6. Section 5 introduces the theoretical foundations, observables, and conventions relevant for subsequent sections. Sections 6, 7 and 8 describe predictions of SGWBs sourced by first-order phase transitions, cosmic strings, and inflationary processes.
Section 9 explores how these diverse SGWB signals convey unique information on the expansion rate of the universe at redshift ~ 1000 or higher. These latter four sections touch on topics that are crucial for LISA science objective SO7 “understand SGWB and their implications”. Inflation not only leads to a SGWB but also to density perturbations which may give rise to the formation of primordial black holes (PBHs), which is the subject of Sect. 10. Finally, Sects. 11 and 12 present existing or planned tools and methods to analyse the GW signals discussed in the previous sections. Such tools and methods potentially constitute key deliverables for SO1-SO7 as well as LISA science objective SO8 “search for GW bursts and unforeseen sources”.

We emphasize that the analyses presented in these final sections need to take into account the superposition between cosmological and astrophysical sources, so that the identification of a potential cosmological contribution necessarily involves modeling of the astrophysical signal. While we briefly discuss the current status of these analyses, we acknowledge that the further development of the astrophysical modeling is by itself a crucial line of reasearch that, while impacting the cosmological searches, goes well beyond what is described in the present document. More in general, certain topics of cosmological interest are intentionally omitted from this document: dark matter particles, some tests of general relativity (GR), waveform uncertainties, astrophysical backgrounds, and the astrophysics of discrete sources such as MBBHs. These and related topics are covered by the Living Reviews maintained by the LISA Astrophysics (Amaro-Seoane et al. 2023), Data Challenge (LISA Data Challenge Working Group n.d.), Fundamental Physics (Arun et al. 2022) and Waveform Working Groups (LISA Waveform Working Group n.d.). Such reviews complement the picture presented herein, by the Cosmology Working Group. The goal of all these documents is to both identify LISA science objectives and corresponding work packages, and to alert the scientific community about novel research opportunities, or potential gaps. The tests of general relativity of Sect. 4 and PBH science of Sect. 10 are exemplary cases of why these Living Reviews are needed. The original LISA proposal (Amaro-Seoane et al. 2017) makes no mention of these science investigations. But in recent years, as the subjects have evolved, a set of new science objectives have been proposed to cover them. The possibility that similar situations arise again justifies the effort and interest for living reviews that report on the thrilling, blooming, and fast-evolving LISA science.

## 2 Tests of cosmic expansion and acceleration with standard sirens

Section coordinators: J.M. Ezquiaga, A. Raccanelli, N. Tamanini. Contributors: D. Bacon, T. Baker, T. Barreiro, E. Belgacem, N. Bellomo, D. Bertacca, C. Caprini, C. Carbone, R. Caldwell, H-Y. Chen, G. Congedo, M. Crisostomi, G. Cusin, C. Dalang, W. Del Pozzo, J.M. Ezquiaga, N. Frusciante, J. Garcia-Bellido, D. Holz, D. Laghi, L. Lombriser, M. Maggiore, M. Mancarella, A. Mangiagli, S. Mukherjee, A. Raccanelli, A. Ricciardone, O. Sergijenko, L. Speri, N. Tamanini, G. Tasinato, M. Volonteri, M. Zumalacarregui.
## 2.1 Introduction

Broadly speaking, to learn about the universe and its cosmic expansion we need to measure distances and times. GW astronomy offers a unique perspective in this matter, since the signal emitted by a compact binary coalescence is well predicted by GR. Namely, the amplitude of the GW is inversely proportional to its luminosity distance and it only depends on the masses and orbital inclination of the binary system source. Since cosmological propagation at the background level (namely, excluding the effect of perturbations over the Friedmann–Lemaître–Robertson–Walker (FLRW) geometry) only changes the overall strain amplitude, one can use the frequency evolution of the GW to unveil the masses of the compact binary and the relative amplitude of the two GW polarisations to estimate the orbital inclination, obtaining thus a direct and absolute measurement of the luminosity distance. However, GW signals alone do not provide a way to relate the time (of merger, for example) in the observer frame to the one in the source frame. To access this information, one needs an independent determination of the redshift of the source. In such a case the GW signal from compact binary coalescence can be considered a standard siren (Schutz 1986), namely a cosmological event for which a distance measurement and complementary redshift information are both available. For example, the binary neutron star (BNS) merger GW170817, observed by the Advanced LIGO and Advanced Virgo detectors jointly with several EM facilities which spotted associated EM emissions, has already been used as a proof-of-principles, low-redshift measurement of $H_0$ (Abbott et al. 2017a). On the other hand, massive black holes (MBHs) seen by LISA with an EM counterpart could be used to map the cosmic expansion up to high redshift (Holz and Hughes 2005; Tamanini et al. 2016; Tamanini 2017; Belgacem et al. 2019c).

In this section we will present the different standard sirens that LISA will detect and the information about the cosmological model that they will provide. The section is organized as follows. We begin by describing the concept of standard siren in Sect. 2.2, detailing the expected LISA bright and dark sirens. We then consider the constraints that could be placed in the standard cosmological constant plus cold dark matter (ΛCDM) cosmological model in Sect. 2.3. Subsequently, we explore LISA capabilities to probe different dark energy (DE) models in Sect. 2.4. Next, we show the synergies of LISA with other EM and GW observatories in Sect. 2.5. Finally, we describe the benefit of cross-correlating LISA data with large-scale structure surveys in Sect. 2.6.

## 2.2 Standard sirens

GW signals from compact binary coalescences are natural cosmic rulers because of the inverse dependence of the strain with the GW luminosity distance, $h \propto 1/D_L^{GW}$. In GR and over a Friedman–Lemaître–Robertson–Walker (FLRW) background, the GW luminosity distance is equal to the one of EM radiation and given by
$$D_L^{\text{GW}} = \frac{c}{H_0} \frac{(1+z)}{\sqrt{|\Omega_k|}} \text{sinn} \left[ H_0 \int_0^z \frac{\sqrt{|\Omega_k|}}{H(z')} dz' \right], \quad (1)$$

where sinn(x) = sin(x), x, sinh(x) for a positive, zero and negative spatial curvature respectively. Assuming a ΛCDM cosmology, the Hubble parameter is a function of the matter content Ω_m, the curvature Ω_k and the amount of DE Ω_Λ (radiation at present time is negligible)

$$H(z) = H_0 \sqrt{\Omega_m(1+z)^3 + \Omega_k(1+z)^2 + \Omega_\Lambda}. \quad (2)$$

LISA will attempt to measure H₀, Ω_m, Ω_k and Ω_Λ using sirens.

GW observations themselves, however, do not provide direct information about the redshift. Therefore, in order to be able to probe the cosmological evolution we need additional input. In the case in which the redshift of the GW source is directly obtained from an EM counterpart, we will refer to the source as a *bright siren*. A beautiful example of this kind of multi-messenger event was the LIGO/Virgo event GW170817 (Abbott et al. 2017d), which provided the first standard siren measurement of H₀ (Abbott et al. 2017a) (see Abbott et al. 2021c for an update on the measurement of H₀ from standard sirens). As we present in Sect. 2.2.1, LISA will be sensitive to very different bright sirens, but the concept remains the same. On the other hand, when the redshift information is obtained from an analysis that does not include EM counterparts, we will refer to the GW sources as *dark sirens*. In Sect. 2.2.2 we will present different dark siren classes that LISA will detect and which can be used to obtain cosmological information by cross-matching the sources with galaxy catalogues and looking for correlated features in the mass distribution.

Modern analyses of standard sirens are based on Bayesian inference. We refer the reader to Sect. 11.1 for a glimpse at the actual statistical tools LISA will use and a detailed discussion of their associated systematic uncertainties. In what follows we focus mostly on the different GW sources and their potential as standard sirens. LISA will detect three types of potential standard siren populations: MBBHs at $1 \lesssim z \lesssim 8$, extreme mass ratio inspirals (EMRs) at $0.1 \lesssim z \lesssim 1$ and SOBBHs at $z \lesssim 0.1$. An example of the expected Hubble diagrams from these three different standard sirens populations can be found in Fig. 1. The details regarding each of these populations will be presented in the following.

### 2.2.1 Bright sirens: MBBHs with electromagnetic counterpart

LISA will detect the coalescence of MBBHs up to redshift $z \simeq 15-20$. However, the mass and redshift distributions of the events are still uncertain. Currently, our knowledge of MBHs is limited to cases where either an active galactic nucleus is present or to quiescent MBHs in nearby galaxies (Kormendy and Ho 2013; Graham 2016; Greene et al. 2020). The population of MBBHs accessible by LISA might be considerably different from our current expectations. Several groups have attempted to address this question with hydrodynamics simulations (Salcido et al. 2016; Katz et al. 2020; Volonteri et al. 2020) or semianalytic formation models (Volonteri et al.
<figure><img src="image_5.png" /><figcaption>**Fig. 1** Examples of LISA standard siren data sets at different redshift ranges. Data from low-redshift dark sirens, defined by SOBBHs, are reported in the bottom-left plot. Data at intermediate redshifts correspond to EMRIs and are reported in the top plot. High-redshift data are provided by standard sirens from MBBHs and are reported in the bottom-right plot. Note how SOBBHs and EMRIs, being dark sirens, can only provide broad likelihood regions in the Hubble diagram, while MBBHs, being bright sirens, provide precise data points thanks to the unique redshift association coming from the EM counterpart identification. Images reproduced with permission from [top] Laghi et al. (2021), [bottom left] Del Pozzo et al. (2018), and [bottom right] Speri et al. (2021), copyright by the author(s)</figcaption></figure>

2003; Barausse et al. 2020; Ricarte and Natarajan 2018; Trinca et al. 2022; Valiante et al. 2018). While the former are able to handle more naturally hydrodynamical, thermodynamical, and dynamical processes, the latter are computationally efficient and can be used to explore a larger parameter space.

Two main sources of uncertainties affect the expected redshift and mass distribution of merging MBBHs: black hole seeding and delay time prescription (Klein et al. 2016). If MBHs grow from the remnants of metal-poor population-III stars, the population of MBBHs accessible by LISA is expected to peak at the total mass $M_{\text{tot}} \simeq 10^3 M_{\odot}$. However if MBHs arise from the monolithic collapse of gas in protogalaxies, the mass distribution is expected to range from $10^4 M_{\odot}$ up to few $10^7 M_{\odot}$. We note that additional formation mechanisms have been proposed and that they would further modulate the distribution of merging MBHs. Moreover, delay times, between the merger of two galaxies and the merger of their central MBHs,
shape the redshift distribution, with short delays leading to more mergers at higher redshift. Further uncertainties arise from the gas inflow to the halo centre, its efficiency and the geometry of accretion. Even if LISA will be able to distinguish different formation scenarios (Sesana et al. 2011), the aforementioned uncertainties reflect in a broad interval for the number of events detected and their distributions (e. g. mass, mass ratio, spins, redshift).

Combining these uncertainties, LISA should be able to detect between a few and several tens of MBBH events per year (Klein et al. 2016). Multiple-body interactions among a MBH binary and one or more intruder MBHs, arising naturally from the hierarchical galaxy formation process, might still produce $\simeq 10-20$ events per year (Bonetti et al. 2019).

LISA will also provide exquisite accuracy on MBBH parameters. For the search of a possible EM counterpart, the sky position accuracy is of paramount importance. Taking into account the full inspiral-merger-ringdown GW signal, MBBHs from few $10^5 M_{\odot}$ to few $10^6 M_{\odot}$ can be localised within $0.4 \text{ deg}^2$ up to $z \simeq 3$ (Mangiagli et al. 2020), but with high accuracy obtained only at merger. For these sources the posterior on the sky position is expected to be Gaussian; however, for more massive and distant sources, the recovered sky position is expected to present multimodalities (Marsat et al. 2021). For heavy systems with total mass $M_{\text{tot}} > 10^7 M_{\odot}$ the ringdown portion of the signal might carry most of the information for the source localisation (Baibhav et al. 2020). For cosmology applications, also the estimate on the luminosity distance plays a fundamental role: due to the typical large signal-to-noise ratio (SNR) value for these sources, LISA should be able to constrain the luminosity distance to better than 10% for most of the events at $z<3$ (Tamanini et al. 2016).

If MBBHs evolve in gas-rich environment, EM radiation might be produced by the accretion of gas onto the MBHs close to or after merger. The orbital motion of the binary is expected to open a cavity in the circumbinary disk. Hydrodynamical simulations show that minidisks generally form around each BH from the stream of gas from the circumbinary disk (Farris et al. 2015; Tang et al. 2018).

This leads to pre-merger EM emission across all wavelengths (D’Ascoli et al. 2018), which can be identified if the pre-merger sky localisation is good (Dal Canton et al. 2019). In the optical band, the Vera C. Rubin Observatory (formerly known as Legacy Survey of Space and Time) (Abell et al. 2009) will reach a magnitude limit of 24.5 in 30 s of pointing over a field of view of $\simeq 10 \text{ deg}^2$. This survey speed enables the Vera C. Rubin Observatory to cover a sky area of $\simeq 100 \text{ deg}^2$ allowing for possible pre-merger EM detection, though the number of detected counterparts is expected to be low (Tamanini et al. 2016). At or after merger, several transients have been proposed, from spectral changes and brightening (Schnittman and Krolik 2008; Rossi et al. 2010) to jets (Palenzuela et al. 2010; Khan et al. 2018; Yuan et al. 2021). In X-ray, Athena (Nandra et al. 2013), with a field of view of $0.4 \text{ deg}^2$ and a limiting flux of $\approx 3 \times 10^{-16} \text{ erg cm}^{-2} \text{ s}^{-1}$ in 100 ks, will be optimal to search for possible post-merger signatures. Similarly, in radio, SKA (Dewdney et al. 2009) will observe the launch of putative post-merger radio jets with an initial field of view of $1 \text{ deg}^2$.

A counterpart detection strategy has been proposed (Tamanini et al. 2016) consisting in first localising the source in the radio with the SKA, and subsequently
to proceed to a redshift determination of the host galaxy in the optical with the Extremely Large Telescope.<fnref n="1" /> This strategy is more promising than direct optical identification with the Vera C. Rubin Observatory. Depending on the seed and dynamical evolution models, LISA will detect between ≈ 10 and ≈ 25 MBBH events with EM counterpart during a mission assumed to be of five year (Tamanini et al. 2016). These estimates are however affected by several astrophysical uncertainties that have as yet not been properly characterised. It is, however, robust to expect that MBBH bright sirens will all be detected at relatively high redshift: one study finds standard sirens distributions peaking at around $z \sim 2-3$ and extending up to $z \sim 8$ (Tamanini et al. 2016; Tamanini 2017). Therefore, MBBH bright sirens will be of great relevance to test the ΛCDM model, and possible deviations from it, in a redshift range so far scarcely probed by EM observations. In Sects. 2.3 and 2.4 we review the cosmological constraints that LISA will be able to impose both at low and high redshift.

### 2.2.2 Dark sirens: SOBBH, EMRs, IMBBHs

In addition to observing the EM counterpart of individual events, there are other methods that one can employ to obtain information about the redshift of the GW source. The most common and widely-used among these methods relies on statistical matching the inferred position of the GW source with a catalogue of galaxies with known redshift (see Abbott et al. 2021c for the application of this method to the Advanced LIGO-Advanced Virgo-KAGRA data). GW events for which an EM counterpart cannot be identified, but which can still be used to extract cosmological information statistically, are usually referred to as dark standard sirens. In this section we outline how dark sirens are treated by correlating galaxy catalogues with the localisation of the source, *statistical dark sirens*, and using properties of their mass distribution, *spectral sirens*. Generally speaking dark standard sirens have the advantage to be applicable to all kind of GW sources for which a distance measurement can be retrieved (not exclusively those emitting EM counterpart signals), although a large number of them are necessary to achieve precise measurements through solid statistics.

In the absence of an EM counterpart, redshift information can be extracted by putting a prior on potential hosts from a galaxy catalogue, assuming that galaxies are good tracers of binary black hole (BBH) mergers. In order to do that (see more details on statistical methods in Sect. 11.1.2), one associates the GW event with every galaxy within the 3D localisation error of the event, assigning to each of them a certain probability of being the true host galaxy of the GW event. In this way, by stacking together the information gathered from several dark sirens, one can statistically infer the true values of the cosmological parameters. For this method to be effective, one requires a large number of events to combine statistics or a very small localisation volume. Errors in the luminosity distance scale with the SNR, $\Delta D_L/D_L \sim 1/\rho$, while the localisation depends also on the duration of the source, since LISA orbital modulation can be used to help disentangle the location in the sky. One complication

<footnote n="1">¹ https://elt.eso.org/.</footnote>
is given by the fact that galaxy catalogues are in general not complete, especially at
high redshift. This requires one to include information on the missing galaxies within
the catalogue, and the catalogue prior must be supplemented by a suitable “comple-
tion” term (Chen et al. 2018; Fishbach et al. 2019; Gray et al. 2020; Finke et al. 2021)
(again see Sect. 11.1.2 for more details). This aspect will remain a limiting factor
until catalogues with very large completeness are available. The availability of
complete galaxy catalogues, small GW localisation regions and accurate redshift
determination—including characterisation of the uncertainty due to peculiar
velocities for low-redshift binaries (Howlett and Davis 2020; Mukherjee et al.
2021c; Nicolaou et al. 2020)—will be crucial in order for the statistical method to
give competitive constraints on cosmological parameters. At the time when LISA
will be taking data, galaxy catalogues will be available from a plethora of current and
future experiments, providing observations of different types of galaxies, with
varying number density and sky coverage, over different redshift ranges. In
particular, the completed observations from Euclid and Dark Energy Spectroscopic
Instrument (DESI), with the addition of redshift-deep catalogues from the Subaru
Prime Focus Spectrograph project and the Roman Space Telescope, should be
available before LISA is launched. Additionally, the Vera C. Rubin Observatory
should have observed hundreds of millions of galaxies, and deep and wide catalogues
should be available from SPHEREx. Finally, on time-scales comparable with LISA,
the full SKA2 and the ATLAS satellite should provide extremely deep and full-sky
catalogues of the sky. Moreover, there are plans to build a next-generation billion-
galaxies survey as a successor of DESI. These future surveys are expected to provide
well-complete galaxy catalogs up to $z = 1$ and possibly beyond.

Stellar-origin black holes (SOBHs) are guaranteed dark sirens for LISA. From the
observations with current LIGO/Virgo detectors (Abbott et al. 2021e), we know that
there is a population of BBHs with masses between ~5 $M_\odot$ and ~50 $M_\odot$ that LISA
will see in their early inspiral. Some of those events will be subsequently detected by
ground-based detectors becoming in this way “multi-band” events (see Sect. 2.5.2).
SOBBHs with good enough localisation can be used as dark sirens (Kyutoku and
Seto 2017; Del Pozzo et al. 2018). Because of their low masses, cosmologically
useful SOBHs could only be seen by LISA up to $z \sim 0.1$, probing essentially the local
expansion rate $H_0$. In practice, only those events with better than 20% accuracy in the
luminosity distance and with a typical sky-localisation error better than ~1deg² will
be useful as dark sirens. According to current forecasts, LISA could detect ~10–
100 dark sirens (Del Pozzo et al. 2018), though uncertainties on the LISA noise level
at the high-frequency end of its band (Amaro-Seoane et al. 2017), on the detection
threshold (Moore et al. 2019) and merger rate of these systems (Abbott et al. 2021e),
might well invalidate the most optimistic expectations. In Sects. 2.3 and 2.4 we will
review the cosmological constraints that LISA can impose with SOBBHs as dark
sirens.

LISA will also be able to use EMRIs as statistical dark sirens. They will in fact be
detected at cosmological distances, possibly in high numbers, though the rates are so
far extremely uncertain. EMRIs detected by LISA are expected to be broadly peaked
around $0.5 < z < 2$, potentially reaching $z \sim 4$ (Babak et al. 2017). A recent
investigation showed that events up to $z \sim 0.7$ can safely be used to estimate $H_0$ (Laghi et al. 2021), provided they can be detected with high SNR. These events can reach relative uncertainties on $\Delta D_L/D_L$ below 0.05 and typical sky-localisation errors less than 1 deg$^2$, representing the best well-localised events suited for a statistical approach on the inference of the GW redshift, comparable to, if not better than, what future third-generation ground-based detectors can achieve (Iacovelli et al. 2022). As we will see in Sects. 2.3 and 2.4, according to the analysis of Laghi et al. (2021), LISA will be able to use from a few to several tens of EMRs to extract useful cosmological information.

Redshift information on a GW source can also be obtained in a statistical way performing a population analysis when there are distinctive features in the source mass distribution. This is simply because GW observatories are only sensitive to the redshifted or detector frame masses, which directly relate to the source masses via

$$m_z = (1 + z)m. \qquad (3)$$

Therefore, if the mass distribution presents a known feature, e.g. a peak or a drop, at a reference scale which is invariant under cosmic evolution, by observing this feature in the GW events at different luminosity distance bins one can infer their redshift, becoming spectral sirens. If such a feature exists, this would be a very convenient probe of the cosmic expansion because it only requires GW data. In case of a time evolution of the reference mass scale, using the full mass distribution helps calibrating against astrophysical uncertainty (Ezquiaga and Holz 2022).

A good example of such features occurs in the mass spectrum of SOBBHs. This is because as stars become more massive, a runaway process induced by electron-positron pair production known as pair-instability supernova (PISN) is triggered (Barkat et al. 1967; Fowler and Hoyle 1964; Heger and Woosley 2002; Fryer et al. 2001; Heger et al. 2003; Belczynski et al. 2016). These PISN result in complete disruption of the stars, preventing the formation of remnant BHs and thus inducing a gap in the mass spectrum starting at around $50 M_\odot$. Nonetheless, for sufficiently massive stars the PISN process is insufficient to prevent direct collapse, and a population of BBHs with masses larger than $\sim 120 M_\odot$ could arise. Therefore, the theory of PISN predicts a gap in the BBH mass spectrum with two edges that act as reference scales.$^2$ While the lower edge of the gap lies within the main sensitivity of present LIGO/Virgo detectors and could lead to precise measurements of $H(z)$ (Farr et al. 2019; Ezquiaga and Holz 2022), LISA will be more sensitive to the upper edge if a population of “far side” binaries in fact exists (Ezquiaga and Holz 2021). These alternative methodologies are currently under development and will need to be further investigated in the future, especially in the framework of GW cosmology with LISA.

Finally, we mention another effect that in principle allows one to access the redshift information directly from the GW signal alone. The variation of the background expansion of the universe during the time of observation of the binary induces an effectively $-4$ post-Newtonian term in the waveform phase, whose
amplitude directly depends on the redshift and on the value of the Hubble parameter both at the source and at the observer (Seto et al. 2001; Nishizawa et al. 2012). Unfortunately, redshift perturbations due to the inhomogeneous distribution of matter between the source and the observer also depend on time, and therefore also contribute to the extra terms in the phase (Bonvin et al. 2017). Among these, the time-varying peculiar velocity of the GW source centre of mass may dominate the signal, effectively preventing the extraction of the redshift information from the amplitude of the dephasing—but possibly allowing a measurement of the binary’s peculiar acceleration (Tamanini et al. 2020; Inayoshi et al. 2017).

### 2.2.3 Systematic uncertainties on standard sirens

Bright and dark sirens will suffer from some common systematic uncertainties. The measurements of the binary luminosity distances are affected by the detector calibration uncertainty (Karki et al. 2016; Chen et al. 2021a) and the accuracy of the waveform models (see Abbott et al. 2019b for discussions in the context of LIGO/Virgo observations). Accurate waveforms will be particularly needed for the high SNR sources that LISA will detect. Moreover, high redshift sources will be affected by weak lensing uncertainties (Holz and Linder 2005; Hirata et al. 2010; Cusin and Tamanini 2021) (see Sect. 3 for more details). In addition, the parameter estimation of the luminosity distance will be subject to degeneracies with the orbital plane inclination and other parameters, though for long duration signals or when higher harmonics are measured (Baibhav et al. 2020), this degeneracy can be broken. In general, these parameter degeneracies will lead to larger statistical errors but not systematic biases. Finally, our understanding of the possible observational selection effect (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 standard siren analysis as well. Not many investigations have so far assessed the systematic uncertainties affecting standard siren measurements with LISA. A thorough exploration of all these effects, needed to consolidate our confidence on LISA cosmological observations, will be necessary in the future.

## 2.3 Constraints on $\Lambda$CDM

In this subsection we present how well LISA will be able to constrain the cosmological parameters of the standard $\Lambda$CDM model, by using different classes of standard sirens as presented in Sect. 2.2. We first focus on the Hubble constant $H_0$ and consider constraints on additional parameters afterwards. We conclude the subsection with a discussion on consistency tests of $\Lambda$CDM at high-redshift with LISA MBBH standard sirens.

<footnote n="2">² Recent analyses have shown however that both the lower end of the gap (Farmer et al. 2019) and its width (Farmer et al. 2020) are robust against ambient factors and nuclear reaction rates.</footnote>
### 2.3.1 $H_0$ tension and standard sirens

The standard model of cosmology is extremely successful and allows one to describe the universe from the time of BBN to the present time of cosmic acceleration. Remarkably, it contains only six parameters, one of which, the present-day Hubble constant $H_0$ describes the expansion rate of the universe. At small redshifts, it relates the luminosity distance and the redshift of a source such that $D_L(z \ll 1) = cz/H_0$. Consistency of the model requires the inferred value of $H_0$ to be independent of the probe and any deviations should be seen as a sign of unaccounted for systematics or more excitingly, new physics.

In recent years, two sets of values for $H_0$ have emerged from so-called early or late measurements of $H_0$ depending on the origin of the calibration. As the error bars shrink, it becomes increasingly clear that the two values are in tension, reaching $4-5\sigma$ disagreements (Verde et al. 2019). The most precise measurement of $H_0$ from the early universe comes from the cosmic microwave background (CMB) (at recombination redshift $z \sim 1100$) with an inferred value of $H_0 = 67.4 \pm 0.5 \text{ km s}^{-1} \text{ Mpc}^{-1}$ at 68% C.L. assuming a flat $\Lambda$CDM cosmology (Aghanim et al. 2020). In this case the $H_0$ value is inferred from the angle upon which the scale associated to the horizon at the last scattering surface is projected, which is obtained from the measurement of the density fluctuations. Compatible values of $H_0$ are obtained also from the Atacama Cosmology Telescope and WMAP5 for which $H_0 = 67.6 \pm 1.1 \text{ km s}^{-1} \text{ Mpc}^{-1}$ at 68% C.L. (Aiola et al. 2020) and from the joint analysis of Dark Energy Survey (DES) clustering and weak lensing data with baryon acoustic oscillations (BAO) and big bang nucleosynthesis (BBN), $H_0 = 67.2_{-1.0}^{+1.2} \text{ km s}^{-1} \text{ Mpc}^{-1}$ at 60% confidence (Abbott et al. 2018c).

In contrast, several teams have measured a significantly higher value of $H_0$ in the local universe with redshifts $z \le 1$ in a model independent fashion. For example, the SH0ES team used Cepheid calibrated supernovae type Ia to measure $H_0 = 74.03 \pm 1.42 \text{ km s}^{-1} \text{ Mpc}^{-1}$ (Riess et al. 2019) (see also recent updates, Riess et al. 2022). The H0LiCOW collaboration used strong lensing time delays of background quasars to infer $H_0 = 73.3_{-1.8}^{+1.7} \text{ km s}^{-1} \text{ Mpc}^{-1}$ (Wong et al. 2020). The Megamaser Cosmology project used very long baseline interferometric observations of water masers in Keplerian orbits around MBHs to measure $H_0 = 73.9 \pm 3.0 \text{ km s}^{-1} \text{ Mpc}^{-1}$ (Pesce et al. 2020). The Carnegie-Chicago Hubble Program collaboration used tip of the red giant branch measurements in the large Magellanic cloud to calibrate 18 supernovae type Ia, instead of Cepheids, and found a slightly lower $H_0 = 69.8 \pm 2.6 \text{ km s}^{-1} \text{ Mpc}^{-1}$ (Freedman et al. 2019).

GWs offer an independent test of the tension using bright or dark sirens as described in Sect. 2.2.1 and Sect. 2.2.2. The LIGO-Virgo collaboration used the first bright siren, namely GW170817, to infer $H_0 = 70_{-8}^{+12} \text{ km s}^{-1} \text{ Mpc}^{-1}$ at $1\sigma$ (Abbott et al. 2017a), which lies somewhat in between the early and late universe values but with a worse precision if compared to current EM results. Furthermore, current dark siren measurements reached an inferred value of $H_0 = 75_{-22}^{+25} \text{ km s}^{-1} \text{ Mpc}^{-1}$ at $1\sigma$ (Fishbach et al. 2019; Finke et al. 2021; Abbott et al. 2021a). The rather large error bars are expected to shrink as $1/\sqrt{N}$, where $N$ indicates the number of well-localized
events. Percent-level precision on $H_0$ is expected in the 2020s with BNS mergers detected by Advanced LIGO-Virgo and their EM counterpart observations (Nissanke et al. 2013; Chen et al. 2018). Bright sirens in the era of the third generation of ground-based GW detectors could further constrain other cosmological parameters, such as $\Omega_m$ and $w_0$ (Sathyaprakash et al. 2010; Zhao et al. 2011; Cai and Yang 2017; Jin et al. 2020; Chen et al. 2021a) and usher us into the era of precise GW cosmology.

LISA will offer an alternative and complementary probe of ΛCDM which might as well provide useful information on the Hubble tension. In the following, we explore the potential of LISA to probe $H_0$ and beyond.

### 2.3.2 LISA forecast for $H_0$

LISA will be able to contribute measurements of $H_0$ coming from different classes of standard siren sources (see Sect. 2.2). For the time being, the literature has described only measurements coming from individual classes of sources. A complete analysis combining the constraining power of different LISA GW sources is still missing.

By considering SOBBHs, and by cross-matching with simulated galaxy catalogues, Kyutoku and Seto (2016), Del Pozzo et al. (2018) found that constraints on $H_0$ can reach the few % level. In particular, in the study presented in Del Pozzo et al. (2018), several different instrumental configurations for LISA were investigated, as well as several coalescence rate models within the range allowed by the LIGO-Virgo observations from O1. This analysis also had $\Omega_m$ as a free parameter, although due to the low-redshift of the sources results are not informative. The SOBBHs entering the analysis were selected to have SNR > 8, a cosmological redshift < 0.1 and an uncertainty on the luminosity distance smaller than 20%. No other selection criteria were applied. With the aforementioned selections, the number of SOBBHs considered ranged from a pessimistic case of 7, yielding an accuracy on $H_0$ of 7% to a most optimistic case of 259, yielding an accuracy of 1%, see left panel in Fig. 2.

A similar analysis can be done also with EMRIs. A first investigation (MacLeod and Hogan 2008) pointed out that ~20 EMRIs detected at $z \sim 0.5$ could be enough to constrain $H_0$ at the 1% level. The analysis provided in MacLeod and Hogan (2008) employed however a simplified approach to estimate cosmological forecasts with LISA, and moreover assumed the more optimistic mission design considered at the time. A recent, more detailed analysis has been performed with LISA EMRIs (Laghi et al. 2021). Using only the most informative, high-SNR ($> 100$) EMRIs up to redshift $z \le 0.7$, cross-matching with the galaxy catalogue obtained from the simulated sky of Henriques et al. (2013), it is shown that constraints at the few % can be forecast for $H_0$. An analysis of three different EMRI population models taken from Babak et al. (2017), representing a pessimistic, a fiducial, and an optimistic scenario, points out that in a 4-year LISA mission lifetime constraints are expected to be at 3.6%, 2.5%, and 1.6% (68% CL), respectively, while in case of a 10-year mission $H_0$ can be constrained at the 2.6%, 1.5%, and 1.1% accuracy (68% CL). The different accuracy in the various scenarios reflects the different number of useful EMRIs available in each model, which in case of 10 years of observation and after the SNR
<figure><img src="image_6.png" /><figcaption>**Fig. 2** Left panel: 90% (black) and 68% (red) credible intervals for $H_0/\text{km s}^{-1}\text{Mpc}^{-1}$ for each of the LISA configurations considered in Del Pozzo et al. (2018). The credible regions are averaged over the galaxy hosts realisations. Right panel: 90% (black) and 68% (red) percentiles, together with the median (red dot) of $h$ for the pessimistic, fiducial, and optimistic EMRI models (M5, M1, M6, respectively) and for two different LISA observational scenario (4 and 10 years). Here $h = H_0/(100\text{km s}^{-1}\text{Mpc}^{-1})$ and the blue dashed horizontal line denotes the true cosmology (set at $h = 0.73$ in Del Pozzo et al. 2018; Laghi et al. 2021). For each data point, we also report the average number $N$ of EMRIs considered in the analysis. Images reproduced with permission from [left] Del Pozzo et al. (2018) and [right] Laghi et al. (2021), copyright by the author(s)</figcaption></figure>

selection, ranges from ~5 (in the worst scenario), passing to ~30 (in the fiducial scenario), up to ~70 (optimistic scenario), see right panel in Fig. 2.

Finally the last standard siren class that LISA can employ to constrain $H_0$ are MBBHs. Although these events are expected to be detected at high-redshift ($z > 1$), by assuming $\Lambda$CDM one can set bounds on $H_0$, i.e. on the cosmic evolution at low-redshift. By simulating different populations of MBBH mergers, performing (simple) parameter estimations over the expected GW signal, and by simulating the emission and detection of possible EM counterparts, recent works showed that constraints at the few % level can be imposed on $H_0$ (Tamanini et al. 2016; Tamanini 2017; Belgacem et al. 2019c). These studies predict around four useful standard sirens with observed EM counterpart, per year. Their redshift distribution peaks between redshift 2 and 4, with tails up to $z \approx 8$. The dominant contribution on the distance uncertainty of these events is not the LISA measurement error, but rather the systematic effect due to week lensing which dominates at high-$z$ providing an estimated average uncertainty of up to 5-10% (Hirata et al. 2010; Cusin and Tamanini 2021). Nevertheless even if LISA will detect only a few MBBH standard sirens, the fact that these will be at high redshift, with relatively precise distance determination, and of course with the single redshift value identified with the EM counterpart, will allow for interesting constraints on $H_0$ at the few percent level (Tamanini et al. 2016; Tamanini 2017; Belgacem et al. 2019c). This is comparable with what is expected for low-redshift, more numerous LISA dark sirens such as SOBBHs and EMRIs. As we will show below, MBBHs will however be more interesting for cosmological analyses beyond the simple measurement of $H_0$.

The joint-inference on $H_0$ resulting from the combination of the analyses described above for SOBBHs, EMRIs and MBBHs, is expected to provide interesting constraints, possibly reaching the 1% level, or better. Such a combined investigation however has not yet been performed and will be the focus of future
studies. Note also that a further class of standard sirens for LISA could be provided by intermediate-mass binary black holes (IMBBHs), ranging from $O(100)$ to $O(1000)$ solar masses. Although recent LIGO/Virgo observations may point towards the existence of this class of BBHs (Abbott et al. 2020b), their merger rate and population properties are completely unknown at the moment (Toubiana et al. 2021), making current LISA forecasts too uncertain to be seriously considered. Given their high-mass and redshift range, IMBBHs could nevertheless well represent the best class of LISA dark sirens if they can be detected in high numbers. Moreover the association with possible EM counterparts, which for these systems are being realistically considered (Graham et al. 2020), could well turn the LISA IMBBH detections that merge within the LIGO/Virgo band in a relatively short time, into useful multi-band bright sirens with a great potential to yield precise cosmological measurements (Muttoni et al. 2022).

By surveying the considerations made above, we can expect LISA to deliver constraints on $H_0$ at the few percent level or better, similarly to what ground-based detectors may produce on comparable timescales (Chen et al. 2021a). The robustness of these constraints against uncertainties in the overall analysis, including for example calibration and waveform modelling issues, will be something to carefully assess in the future. In any case, a precise and accurate measurement of $H_0$ with LISA will provide useful insights on the Hubble tension, should the tension still persist. On the other hand, a further independent and complementary measurement of $H_0$ will help strengthen our confidence in the value of the Hubble constant, especially if hints of physics beyond ΛCDM appear.

### 2.3.3 ΛCDM beyond $H_0$

The relatively high-redshift reach for some of the standard sirens sources detected by LISA, will allow the inference on cosmological parameters beyond $H_0$. In particular MBBHs will be extremely useful to test the cosmic expansion at high-redshift ($z \lesssim 8$) while EMRs could be useful for cosmological applications up to $z \sim 1$, or more generally up to the redshift at which we will be able to employ fairly complete galaxy catalogues.

As shown in Tamanini et al. (2016), Tamanini (2017), MBBH mergers with an identified EM counterpart could be used by LISA to infer the values of $\Omega_m$ and $\Omega_k$, albeit with large uncertainties. Assuming ΛCDM, constraints on $\Omega_m$ are forecast to reach the ~10% level (68% C.L.) only in the most optimistic scenarios (Tamanini et al. 2016; Belgacem et al. 2019c; Spergel et al. 2015), while allowing for spatial curvature degrades these estimates for $\Omega_m$ to ~25% and yields a measurement of $\Omega_k$ at similar precision (Tamanini et al. 2016). Needless to say these results will certainly not be competitive with EM observations, but will at least provide an independent and complementary measurement.

Similarly, EMRs will be able to provide information on $\Omega_m$ only in the most optimistic scenarios. In fact recent estimates (Laghi et al. 2021) indicate that using the loudest $O(70)$ events, as predicted in an optimistic population scenario, it is possible to get constraints with ~20% accuracy (68%CL) on $\Omega_m$ when jointly inferred with $H_0$. Even though this result is not expected to compete with current and
future EM observations (Amendola et al. 2018b; Scolnic et al. 2018), it is somehow
not surprising since the EMRIs used in Laghi et al. (2021) are distributed up to
$z \lesssim 0.7$. Future analysis, possibly including a larger number of events, thus including
low-SNR events which are typically coming from high redshifts, are expected to give
more informative results on cosmological parameters beyond $H_0$.

Preliminary investigations of the full ΛCDM cosmological model where the
curvature term Ω<sub>k</sub> is not fixed to 0, seem to indicate that, even for moderate redshift
sources such as the loudest EMRIs (SNR > 100) at z<sim 0.7, LISA will provide some
simultaneous constraints on all cosmological parameters. In the fiducial scenario of
Laghi et al. (2021), Ω<sub>k</sub> can be constrained with an accuracy of ~30%, while Ω<sub>m</sub> can
be constrained with an accuracy of ~55%, while retaining an accuracy on H<sub>0</sub> of
~2% (all 68% CL) (Laghi 2021). These results will need to be further investigated
in the future, however they are suggestive that a possible simultaneous inference of
different cosmological parameters with LISA standard sirens will indeed be possible.
In particular the full LISA cosmological analysis with results obtained from the
combination of all classes of LISA standard sirens, namely SOBBHs, EMRIs and
MBBHs, should substantially increase the accuracy of the estimates on Ω<sub>m</sub> and Ω<sub>k</sub>
above, thanks especially to the combination of cosmological datasets from different
redshift ranges which might break some of the degeneracy between the cosmological
parameters. A future combined investigation will thus be needed to thoroughly assess
the ability of LISA to constrain parameters beyond H<sub>0</sub>.

2.3.4 Tests of ΛCDM at high-redshift

The most interesting class of LISA standard sirens are MBBHs, not only because
they are expected to produce detectable EM counterparts (bright sirens) but also
because they will be detected at high redshift and thus can be employed to test the
cosmic expansion at early epochs largely unexplored by current EM cosmological
surveys. As we will show in Sect. 2.4 and in Sect. 4, MBBHs will have the potential
to test different cosmological models alternative to ΛCDM. Here we briefly mention
how well LISA can test ΛCDM itself at high z, taking into account possible general
deviations (to be discussed shortly) and by comparing with EM probes at similar
redshifts. Note also that in analogy to EM distance observations, LISA MBBHs can
as well be employed to probe the fundamental assumptions of ΛCDM, e.g. the
cosmological principle (Cai et al. 2018).

Let us first of all recall how well MBBHs can test ΛCDM. As shown in the
previous section, in the most optimistic cases Ω_m can be tested at the ~10% level,
which reflects the fact that at z > 1 the universe is expected to be matter dominated
and thus to provide information mainly on Ω_m. Any deviation from the standard
matter-dominated evolution at redshift 1≤z≤4 however will be constrained by LISA,
at a level not attained by current EM observations.

To put in context the potential of LISA we can compare its constraining power
with other EM measurements of the cosmic expansion at high z. LISA will in fact
provide independent and complementary data which will not only deliver useful and
accurate information on deviations from ΛCDM but will also be used to check and
cross-validate EM measurements, expected to be sparse and inaccurate at such high-
redshift. As a clear example, LISA can successfully compete with quasar cosmological observations at $z \sim 2$ and above (Speri et al. 2021). Current quasar distance measurements indicate possible issues in the Hubble diagram at $z \gtrsim 2$, where deviations from $\Lambda$CDM have already been claimed (Risaliti and Lusso 2019) (see also Velten and Gomes 2020; Yang et al. 2020b; Banerjee et al. 2021). To understand if such apparent deviations are due to systematic effects or are indeed due to new physics beyond $\Lambda$CDM, complementary measurements at the same redshift range will be needed. As shown in Speri et al. (2021), four MBBH standard sirens detected by LISA are on average enough to unequivocally confirm or rule out the apparent deviation claimed in Risaliti and Lusso (2019), and thus to reveal if indeed this is due to systematics in the quasar Hubble diagram or to new physics. Standard siren observations will play a crucial role to test deviations from the $\Lambda$CDM.

GW observations provide a more reliable distance luminosity measurement than other EM observations, such as quasars. Standard sirens rely in fact on fundamental predictions of GR and not on phenomenological relations between observed quantities. Although the quasar Hubble diagram will become more precise and accurate in the time between now and when LISA will fly, for example with observations by the eROSITA (Merloni et al. 2012) and Euclid (Amendola et al. 2018b; Barnett et al. 2019) missions, intrinsic systematics on the quasar cosmological measurement might still be unresolved. LISA will thus provide a unique complementary test of the cosmic expansion at $z \gtrsim 2$, which not only will yield accurate measurement of deviations from $\Lambda$CDM but it will also be used to cross-validate any results obtained by other EM observations at the same redshift range, notably with quasars.

## 2.4 Probing dark energy

This subsection discusses how LISA can probe the fundamental nature of DE through GW standard sirens. Here only simple alternative DE models are considered, and in particular we assume that the underlying gravitational theory remains GR. Models of DE based on modification of GR are discussed in details in Sect. 4.

### 2.4.1 Equation of state of dark energy: $w_0$ and $w_a$

Deviations from the standard cosmological model include the presence of a DE fluid with effective equation of state (EoS) given by: $p_{DE} = w_{DE}(a)\rho_{DE}$, where $p_{DE}$ and $\rho_{DE}$ are the effective pressure and density of the DE fluid respectively. Although the functional dependence of $w_{DE}(a)$ may be very non-trivial (see e.g. self-accelerating cosmologies), for practical purposes we consider here only a simple phenomenological parametrization introduced by Chevalier-Polarski-Linder (CPL) (Chevallier and Polarski 2001; Linder 2003):

$$w_{DE}(a) = w_0 + w_a(1 - a), \quad (4)$$

where $w_0$ and $w_a$ are constants and indicate, respectively, the value and the time
derivative of $w_{DE}$ today. We refer to Sect. 4 about modified gravity for more complicated, model dependent forms.

In all the scenarios studied by Belgacem et al. (2019c) (which depend on the BH seeds and on the assumptions about the error on redshift), the GW sources considered are expected not to contribute relevantly to improve the knowledge on $w_0$, with respect to what is already known from current cosmological observations. More precisely, the $1\sigma$ error $\Delta w_0$ only goes from $\Delta w_0 = 0.045$ (using CMB, type Ia supernovae and BAO data) to $\Delta w_0 = 0.044$ adding MBBH standard sirens to the datasets, even in the best scenario where LISA alone can only reach a 20% relative $1\sigma$ uncertainty on $w_0$ (Tamanini 2017). It is important to remark, however, that these outcomes are based on a mission duration of 4 years and sources are limited to MBBHs with EM counterparts. Significant improvements are expected by extending the data taking time or by combining with information from other sources, notably EMRIs.

Indeed, recent investigations (Laghi et al. 2021) suggest that EMRIs will deliver constraints on $w_0$ of the order of $\sim 10\%$, when inferred simultaneously with $w_a$. When assuming prior knowledge of $H_0$ and $\Omega_m$, constraints on $w_0$ are estimated at the $\sim 7\%$ level in a realistic EMRI scenario, reaching $\sim 5\%$ in the best case (all 90\% C.L.). These results are better than what expected with dark sirens analysis for 3 G ground-based detectors (Belgacem et al. 2018b), and compatible with what they can achieve with bright sirens (Belgacem et al. 2019b). While relevant information on $w_0$ can be obtained with moderately low-redshift events, $w_a$ is expected to be measurable only with higher-redshift events, which in the joint cosmological inference of $w_0$ and $w_a$ of Laghi et al. (2021) are not considered, and thus no relevant measurement of $w_a$ is obtained.

Although from these estimates it seems that LISA standard sirens will not be competitive with future EM observations, we stress that they will anyway provide independent and complementary measurements which will increase our confidence on any insight on the nature of DE. This is strikingly important for modified gravity models of DE, where GWs can indeed provide orthogonal information with respect to EM observations; see Sect. 4.

## 2.4.2 Alternative dark energy models

The CPL phenomenological parametrization discussed in the previous subsection constitutes only a simple approach for detecting deviations from the $\Lambda$CDM. The theoretical description of DE may require more sophisticated models. Here, we briefly discuss such a possibility focusing on scenarios that change the background cosmological expansion with respect to $\Lambda$CDM. We assume the standard evolution equation for GWs. Possible modifications to GW propagation through cosmological spacetimes—motivated by modified graviton dispersion relations, or non-minimal couplings of the dark-energy sector with curvature—are described in Sect. 4.

Models for DE can include quite a large number of parameters, leading to rich dynamics for the DE sector as a function of redshift. Some investigations aim to dynamically explain the puzzling small value for the present-day acceleration rate, leading to a time-dependent evolution of the DE density. These include quintessence
scenarios (Wetterich 1988; Zlatev et al. 1999), which can be generalised to
kinetically-driven (Armendariz-Picon et al. 2000) and kinetic-braiding models
(Deffayet et al. 2010) without modifying the propagation properties of GWs. Other
scenarios aim to alleviate the coincidence problem relating DE with DM at
intermediate redshifts, for example in the DE-DM interacting models (Wetterich
1995; Amendola 2000), or in Chaplygin gas cosmology (Kamenshchik et al. 2001).
More exotic possibilities include holographic DE, associating the present-day
acceleration of the universe with the size of the particle horizon (Li 2004). See e.g.
Copeland et al. (2006b), Li et al. (2011), Ishak (2019), Huterer and Shafer (2018) for
reviews, including analysis of cosmological implications and observational prospects
of DE scenarios. A recent resurgence of interest on DE model building, based on the
previous approaches, has been motivated by the $H_0$ tension discussed in Sect. 2.3.1;
see Di Valentino et al. (2021) for a review. Among many examples, such tension can
be alleviated in scenarios with DM-DE interactions or with features at small or
intermediate redshifts (Di Valentino et al. 2017, 2016; Keeley et al. 2019; Raveri
2020), in ranges that might be probed with GW sirens. See e.g. Knox and Millea
(2020) for a comprehensive discussion on this topic, including comparison between
theoretical ideas and existing cosmological constraints.

These theoretical models suggest that distinctive DE effects can occur at different
redshifts, from very small to relatively large values of z. The capability of LISA to
probe cosmological expansion in a large range of redshifts, as discussed in the
previous sections, provides invaluable opportunities for building independent
cosmological tests (see Sect. 2.3.4) and thus to probe different DE scenarios, in a
complementary way with respect to EM probes. A clear example is given by the
investigations of Cai et al. (2017), Caprini and Tamanini (2016) where LISA
forecasts for testing cosmological models allowing for DE-DM interactions or for
early DE have been produced using MBBH as bright sirens. Similar analyses using
LISA dark sirens are still missing in the literature and will constitute material for
future explorations. By considering the results obtained with both SOBBHs and
EMRIs for standard cosmological models (see Sect. 2.3), it would be very interesting
to further develop these studies by designing an efficient, unified method to
reconstruct the redshift dependence of DE with GWs, similar for example to what
already done with EM observations; see e.g. Sahni and Starobinsky (2006).
Assumptions going beyond the linear CPL parameterisation discussed in the previous
section might be better suited to test specific alternative scenarios, or for improving
the DE reconstruction in a wider redshift interval, e.g. through polynomial fitting
(Alam et al. 2003). Alternatively, methods based on a principal component analysis
of a binned parametrization of the signal as a function of redshift, as proposed in
Huterer and Starkman (2003), might be adapted and applied to GW observations with
LISA.

Besides the alternative DE models already considered in the literature, there are
plenty of others that can still be tested by LISA and for which detailed analyses will
be needed in the forthcoming years in order to understand how well LISA will
constrain them and thus to better assess and expand the science case of the mission.
We conclude by mentioning again that the nature of DE can be further investigated if
this is connected to an underlying gravitational theory beyond GR. In this case new
observational signatures might appear, giving rise to a richer phenomenology and to more promising LISA results. Such models and analyses will be discussed in Sect. 4 in the context of modified theories of gravity (Baker et al. 2021).

## 2.5 Synergy with other cosmological measurements

GW observations by LISA will provide a wealth of cosmological information as we have seen in the previous sections. But LISA will not be alone in the quest of understanding the cosmic history. It is thus of great importance to assess how LISA could complement with other facilities. In what comes next, we discuss LISA synergies with EM observatories and other GW detectors.

### 2.5.1 Integration with standard electromagnetic observations

To date, most studies constraining cosmological parameters with GWs have focused on the Hubble parameter, see Sect. 2.3.2. This is because the cosmological information in GW amplitudes is primarily held in the luminosity distance of the source, $D_L$, which for low redshifts ($z \ll 1$) can be approximated to $D_L \simeq cz/H_0$. The majority of LISA sources will exist at higher redshifts where this approximation breaks down. Instead the full expression for the luminosity distance must be used (see Eq. (1)) and so GW amplitudes in principle have sensitivity to further cosmological parameters traditionally measured electromagnetically, such as the fractional densities $\Omega_m$ and $\Omega_\Lambda$, and the DE EoS parameters $w_0$ and $w_a$, see Sects. 2.3.3 and 2.4.

A key outstanding question is whether specific combinations of EM and GW probes have the ability to improve constraints on these parameters and break degeneracies between them. Details about the galaxy surveys co-temporal with LISA are not available at present, but as a conservative approach we can consider ‘Stage IV’ experiments planned over the next decade, such as DESI, Euclid, the Vera C. Rubin Observatory, and the Roman Space Telescope (Aghamousa et al. 2016; Laureijs et al. 2011; Schmidt et al. 2020; Spergel et al. 2015). The corresponding CMB data will come from the LiteBIRD mission (Hazumi et al. 2012). Direct cross-correlation of GW sources with galaxy catalogues will be considered in Sect. 2.6. Here we focus instead on probe combination, though we note that there is a lack of comprehensive studies on this topic in the current literature.

One advantage leveraged by Stage IV galaxy surveys is the ability to combine multiple probes from the same instrument. Most commonly the main probes are shear power spectra from galaxy weak lensing, BAO measurements from galaxy clustering, supernovae, and, additionally, strong gravitational lenses in some analyses. As examples of the expected constraints, Zhan and Tyson (2018) provides forecasts for standard cosmological parameters using weak lensing, BAO and supernovae data from the Vera C. Rubin Observatory, combined with current Planck CMB data. The resulting 68% confidence intervals on {$w_0, w_a$} and $\omega_m = \Omega_m h^2$ are {$7.05 \times 10^{-2}, 1.86 \times 10^{-1}$} and $7.73 \times 10^{-4}$ respectively, with $h$ being the normalised Hubble parameter. Similar results are expected from the Euclid space
mission, which should yield 68% confidence intervals around
$\{3 \times 10^{-2}, 1.3 \times 10^{-1}\}$ for $\{w_0, w_a\}$. Comparing these values to the forecasts using
LISA EMRI detections presented in Laghi et al. (2021), it seems unlikely that LISA
will be able to offer competitive constraints on $\Omega_m$. However, in mildly optimistic
scenarios they may offer moderate constraints on the DE parameters, with (Laghi
et al. 2021) forecasting 95% confidence intervals on $\{w_0, w_a\}$ of approximately
$\{6.5 \times 10^{-2}, 6.5 \times 10^{-1}\}$ for a four-year LISA mission and optimistic MBH
population models.

An alternative strategy is to find probe combinations where GW data can break
degeneracies existing between EM probes. One such example is put forwards in Qi
et al. (2021), which determines forecasts for the combination BNS data from the
DECihertz Interferometer Gravitational wave Observatory (DECIGO) with mea-
surements of redshift drift (the Sandage–Loeb effect, Loeb 1998) from the SKA and
the European-Extremely Large Telescope (E-ELT). Redshift drift measurements use
high-resolution spectroscopy of the HI emission line (SKA) or Lyman-α absorption
lines in quasar spectra (E-ELT) over long time frames (10 years+) to directly
measure tiny shifts the line frequencies. Although experimentally challenging, this
constitutes a direct measurement of $H(z)$, as opposed to the integrated effect of $H(z)$
probed by luminosity distances; see Eq. (1).

GWs can likewise offer a direct measurement of $H(z)$ through the dipole of the
luminosity distance. Equation (1) gives the luminosity distance-redshift relation in a
perfectly homogeneous and isotropic universe; in fact small corrections to this are
induced by gravitational lensing and peculiar velocities. As shown in Bonvin et al.
(2006a), Bonvin et al. (2006b) for the EM case (see e.g. Nishizawa et al. 2011 for the
GW case), the peculiar motion of our Galaxy with respect to the CMB frame induces
a dipole mode in the luminosity distances. This dipole moment is given by

$$
D_L^{(1)}(z) = \frac{|\mathbf{v}_0| (1+z)^2}{H(z)}, \qquad (5)
$$

where $\mathbf{v}_0$ is the dipole anisotropy in the CMB, estimated to be approximately
$369.1 \pm 0.9$ km/s. Combining measurements of this dipole from DECIGO BNS
sources with the redshift drift measurements above, Qi et al. (2021) finds substantial
breaking of degeneracies in the $H_0 - \Omega_m$ plane, leading to $1-\sigma$ bounds on $\{H_0, \Omega_m\}$
of $\{0.78 \text{ km s}^{-1}\text{Mpc}^{-1}, 0.006\}$, competitive with current bounds. Mild improve-
ments on the constant DE EoS were also obtained ($\sigma_{w_0} \sim 0.03$), though no mean-
ingful constraint on $w_a$ was possible. Further analyses are needed to properly
understand if such method can equally be applied to LISA.

On cosmological scales, a bias prescription is used to model how GW events trace
the large-scale DM distribution. Analogously to galaxy bias, the GW bias is
modelled as scale-independent on large scales. Using the parameterisation
$b_{\text{GW}}(z) = b_{\text{GW}}^0 (1+z)^\alpha$, Mukherjee et al. (2021d) finds the parameters $b_{\text{GW}}^0$ and $\alpha$
to be uncorrelated with $H_0$ and $\Omega_m$ for a $\Lambda$CDM model ($w_0 = -1$ fixed). This
suggests cautious optimism that combined GW and EM constraints on cosmological
parameters should not be strongly sensitive to the GW bias prescription.
## 2.5.2 Complementarity with other gravitational-wave observatories

In addition to the synergies with other cosmological surveys, LISA will also build new synergies with other GW detectors. In particular, LISA will be able to detect the early inspiral phase of compact binary coalescences that later merge within the frequency bands of Earth-based interferometers. Some example signals are displayed in the left panel of Fig. 3. Because the same signal is detected across different frequencies, these events are known as multi-band. Different populations of BBHs could become multi-band sources. Most notably, the SOBBHs that current LIGO/Virgo detectors are observing could have been seen if LISA was online a few years before these detections (Sesana 2016). Nonetheless, LISA high frequency sensitivity limits their number (Moore et al. 2019). If present in nature, IMBBHs would be more promising candidates (Amaro-Seoane and Santamaria 2010; Jani et al. 2019; Sedda et al. 2020), and possibly yield interesting cosmological results (Muttoni et al. 2022). There is however a limit to their masses, because if they are too heavy they will merge before reaching the frequencies of ground-based detectors. (This limitation could be removed with a deci-Hertz observatory (Sedda et al. 2020).)

The right panel of Fig. 3 displays the fraction of multi-band events, defined as the subset of LISA detections merging within 10 years and being detected by a ground-based instrument (Ezquiaga and Holz 2021). For concreteness we consider Advanced LIGO (aLIGO), its possible upgrade (A+), Voyager and the third-generation detectors Cosmic Explorer (CE) and Einstein Telescope (ET). Interestingly, the multi-band fraction peaks where the upper end of the PISN mass gap is expected to be found, implying that far-side binaries could be promising multi-band sources (Ezquiaga and Holz 2021). As noted in Gerosa et al. (2019), for $M_{\text{tot}} \lesssim 100 M_{\odot}$, there is no difference for the multi-band ratio between 2 G and 3 G detectors because the fraction is limited by LISA high-frequency range. On the contrary, for $M_{\text{tot}} > 200 M_{\odot}$ the difference among ground-based detectors are sizeable and depends mostly on their low frequency sensitivity.

Besides individual sources, LISA and ground-based detectors could share SGWBs. For example the background of unresolved SOBBHs could be within the

<figure><img src="image_7.png" /><figcaption>**Fig. 3** On the left, GW signals as seen by LISA and different ground-based detectors. The times indicate the time to merge which depends on the initial frequency. On the right, fraction of total events observed by LISA and that will merge within 10 years, being detected by a ground-based detector. Image [right panel] reproduced with permission from Ezquiaga and Holz (2021), copyright by AAS</figcaption></figure>
reach of LISA detectability (Sesana 2016). Similarly, binaries above the PISN gap could leave an additional background (Mangiagli et al. 2019).

Finally, LISA could also complement other space-based detectors. In particular, there are several proposals such as Taiji (Ruan et al. 2020a) and TianQin Wang et al. (2020a) that could potentially fly at the same time as LISA. These additional detectors would help improve the cosmological inference (Wang et al. 2022a, b; Baral et al. 2020; Zhao et al. 2020). Moreover, any further detector in the deci-Hertz frequency band would be an excellent addition to LISA science (Sedda et al. 2020).

Altogether, LISA and the other GW observatories will all contribute to a standard siren measurement of $H(z)$ across a broad redshift range.

## 2.6 Cross-correlation and interaction with large scale structure

GW maps of resolved events or SGWB can be cross-correlated with other large-scale-structure (LSS) tracers to perform a variety of astrophysical and cosmological measurements. Here we focus on the use of the SGWB measurements while technicalities about this signal are postponed to Sec. 5.

### 2.6.1 Cross-correlations with resolved events

The analysis of the cross-correlations between galaxy surveys and resolved GW events from compact object mergers has a rich scientific potential (Laguna et al. 2010). From binary mergers with an EM counterpart, such an analysis can constrain DE and modified gravity models (Camera and Nishizawa 2013). A similar analysis for sources as the first ones detected by LIGO/Virgo, yields constraints on the distance-redshift relation, the Hubble constant and other cosmological parameters (Oguri 2016). PBH scenarios as well as different astrophysical models can also be tested (Raccanelli et al. 2016; Raccanelli 2017; Scelfo et al. 2018, 2020).

The cross-correlation of the resolved GW sources with galaxies is a further promising avenue. It enables the measurement of the redshift of the GW sources that do not have EM counterpart. This technique makes it possible to measure the value of Hubble constant, the matter density, the DE EoS and its redshift evolution from the GW sources. Along with probing the expansion history of the universe, it also provides evidence on the relation between the GW sources and DM distribution through the redshift-dependent galaxy bias parameter $b(z)$ (Oguri 2016; Raccanelli et al. 2016; Scelfo et al. 2018; Mukherjee and Wandelt 2018; Mukherjee et al. 2021d).

The observable considered for such studies is the (3D) angular power spectrum:

$$C_{\ell}^{\text{GW, LSS}}(z_1, z_2) = \frac{2}{\pi} \int dk k^2 P(k) \Psi_{\ell}^{\text{GW}}(k) \Psi_{\ell}^{\text{LSS}}(k). \quad (6)$$

The two $\Psi_{\ell}^{X}(z)$ kernels encapsulate the physical processes in place:
$$ \Psi_{\ell}^{X}(k) = \int N^{X}(z)b^{X}(z)D(z)W(z)\Delta_{\ell}^{X}(k,z)dz, \qquad (7) $$

where X stands for either GW or LSS when considering one or the other observable, and W(z) are the observational window functions related to the experiment specifications. The Δℓ(k, z) terms are the observed overdensities including effects from the intrinsic clustering, velocity effects (redshift-space distortions and Doppler), lensing and gravitational potentials, respectively, and they include information on galaxy clustering, gravity, and of course cosmological and astrophysical parameters:

$$ \Delta_{\ell}(k, z) = \Delta_{\ell}^{\text{den}}(k, z) + \Delta_{\ell}^{\text{vel}}(k, z) + \Delta_{\ell}^{\text{len}}(k, z) + \Delta_{\ell}^{\text{gr}}(k, z); \qquad (8) $$

for the complete expression of those terms, see e.g. Bonvin and Durrer (2011), Scelfo et al. (2018). In particular, this correlation can constrain DE and modified gravity models (Camera and Nishizawa 2013; Raccanelli 2017). In this case, the modifications due to different models of gravity and DE enter in the Ψℓ kernels, and the advantage of using the GW-LSS cross-correlation comes from the complementarity with other measurements as well as the potentially higher redshift range for the GW bin.

The fact that GW events trace the LSS allows us to also test and constrain astrophysical models (Scelfo et al. 2020; Libanore et al. 2021). In this case, the change in the model Cℓ will be in the merger rate, the redshift distribution and the bias of the compact objects' hosts.

Merging of BHs that are the endpoint of stellar evolution happens almost exclusively in galaxies that had sufficient star formation, and therefore in halos with a relatively large galaxy bias ($b_g > 1$). Conversely, PBHs preferentially merge in lower biased objects, and thus have a lower cross-correlation with luminous galaxies (Bird et al. 2016). Therefore, the cross-correlation of GW maps with galaxies, which can measure the bias of the BBH hosts, provides information on the abundance of PBHs (Raccanelli et al. 2016; Scelfo et al. 2018) (for more details on PBHs, see Sect. 10).

Further correlations contain other useful information. The tomographic shear maps and the number density distribution of GW sources, combined with shear-shear and GW-GW auto-correlations, also constrain the cosmological parameters (Osato 2018), though these require a large number of GW detections which LISA is not expected to observe. Nevertheless thanks to GW-LSS cross-correlations with LISA and future galaxy survey data, the detection of GW weak lensing seems viable (Mukherjee et al. 2020b). The cross-correlation between GW weak-lensing and CMB-lensing might also allow to test fundamental predictions of GR (Mukherjee et al. 2020a).

## 2.6.2 Cross-correlations with the stochastic gravitational-wave background

A complementary study to the cross-correlation between GW resolved sources and LSS is the cross-correlation with the astrophysical SGWB. There are at least two astrophysical backgrounds that LISA will detect: one contribution generated by the galactic binary (GB) mergers, which is expected to dominate at low frequencies (up
to ~$10^{-3}$ Hz), and one coming from extragalactic BBH inspirals, expected to be relevant at larger frequencies (~$10^{-3}$–$10^{-2}$ Hz). Several phenomena in the early universe also source a stochastic signal whose strength is poorly predicted (see Sect. 5 and subsequent sections). The SGWB is then expected to be a combination of an astrophysical and a cosmological component, and a priori any of the two can dominate the SGWB signal in the LISA band. In order to be sensitive to both of them, it is of fundamental importance to find ways to disentangle the two signals. The spectral shape of each contribution is a standard tool to try to disentangle the components, however, due to the richness of sources expected in the LISA band, it is worth finding other ways to characterise the SGWB contributions. A promising approach consists in cross-correlating the SGWB with matter distribution at late times. Since, as we will see, the GW energy density depends on cosmological perturbations (besides astrophysical dependencies), it correlates with other cosmological probes. Some analyses and forecasts of the cross-correlation signals between GW observatories and future galaxy surveys, as e.g. Euclid and SKA, are presented in Contaldi (2017), Cusin et al. (2017), Cusin et al. (2018b), Cusin et al. (2018a), Jenkins et al. (2019a), Jenkins and Sakellariadou (2018), Jenkins et al. (2018), Cusin et al. (2019b), Cusin et al. (2020), Jenkins and Sakellariadou (2019), Jenkins et al. (2019b), Bertacca et al. (2020), Pitrou et al. (2020), Mukherjee and Silk (2020), Alonso et al. (2020b). These cross-correlations not only can be useful to disentangle the origin of the SGWB, but represent completely new observables to infer cosmological information. For instance, along with the spatial fluctuation of the astrophysical SGWB, its temporal fluctuation provides a measurement of the high redshift merger rate of the astrophysical sources contributing to the SGWB (Mukherjee and Silk 2020).

As for the case of resolved events, the observable quantity is the angular cross power spectrum, between the galaxies overdensity and the energy density of the astrophysical SGWB, hereafter labelled AGWB:

$$C_{\ell}^{\text{AGWB} \times \text{LSS}} = \frac{2}{\pi} \int k^2 dk P_{\zeta}(k) \Delta_{\ell}^{\text{AGWB}}(k) \Delta_{\ell}^{\text{LSS}}(k), \quad (9)$$

where $P_{\zeta}(k) = A_s(k/k_{\text{pivot}})^{n_s-1}$ is the scale-invariant curvature power spectrum, with $A_s$ and $n_s$ the amplitude and tilt respectively, while $k_{\text{pivot}}$ is the pivot scale. The two transfer functions, $\Delta_{\ell}^{\text{AGWB}}(k)$ and $\Delta_{\ell}^{\text{LSS}}(k)$, contain astrophysical and cosmological information (see the previous section). The astrophysical information can be included focusing on the anisotropies of the astrophysical SGWB energy density. The total GW energy density per logarithmic frequency $f_0$ and solid angle $\Omega_o$ along the line-of-sight $\mathbf{n}$ of a SGWB is (Cusin et al. 2017; Bertacca et al. 2020)

$$\Omega_{\text{GW}}^{\text{TOT}}(f_0, \mathbf{n}) = \frac{f_0}{\rho_c} \frac{d\rho_{\text{GW}}^{\text{tot}}}{d\rho_c d\Omega_o}. \quad (10)$$

It contains both a background (monopole) contribution in the observed frame $\bar{\Omega}_{\text{GW}}(f_0)$, which is homogeneous and isotropic, and a direction-dependent contribution $\Omega_{\text{GW}}(f_0, \mathbf{n})$. Starting from these, we can define the total relative fluctuation
as $\Delta_{\text{AGWB}}^{\text{TOT}}(f_o, \hat{\mathbf{n}}) = (\Omega_{\text{AGWB}}^{\text{TOT}} - \bar{\Omega}_{\text{AGWB}}^{\text{TOT}})/\bar{\Omega}_{\text{AGWB}}^{\text{TOT}}$. The contributions to the astro-physical SGWB energy density fluctuation are computed in Contaldi (2017), Cusin et al. (2017), Cusin et al. (2018b), Cusin et al. (2018a), Cusin et al. (2019a), Jenkins et al. (2019a), Jenkins and Sakellariadou (2018), Cusin et al. (2019b), Cusin et al. (2020), Jenkins and Sakellariadou (2019), Jenkins et al. (2019b), Bertacca et al. (2020), Pitrou et al. (2020), Alonso et al. (2020b). Here, following Bertacca et al. (2020), we report its expression in the Poisson gauge

$$ \Delta\Omega_{\text{AGWB}}(f_o, \hat{\mathbf{n}}, \theta) = \sum_{[i]} \int dz W(z) \mathcal{F}^{[i]}(f_o, z, \theta) \left\{ b^{[i]} D \right. \\ \left. + \left(b_{\text{evo}}^{[i]} - 2 - \frac{\mathcal{H}'}{\mathcal{H}^2}\right) \hat{\mathbf{n}} \cdot \mathbf{V} - \frac{1}{\mathcal{H}} \hat{\mathcal{D}}_{\parallel}(\hat{\mathbf{n}} \cdot \mathbf{V}) - (b_{\text{evo}}^{[i]} - 3) \mathcal{H} V \right. \\ \left. + \left(3 - b_{\text{evo}}^{[i]} + \frac{\mathcal{H}'}{\mathcal{H}^2}\right) \Psi + \frac{1}{\mathcal{H}} \Phi' + \left(2 - b_{\text{evo}}^{[i]} + \frac{\mathcal{H}'}{\mathcal{H}^2}\right) \int_0^{\chi(z)} d\bar{\chi}(\Phi' + \Psi') \right. \\ \left. + \left(b_{\text{evo}}^{[i]} - 2 - \frac{\mathcal{H}'}{\mathcal{H}^2}\right) \left(\Psi_o - \mathcal{H}_0 \int_0^{\tau_0} d\tau \frac{\Psi(\tau)}{1+z(\tau)}\bigg|_o - (\hat{\mathbf{n}} \cdot \mathbf{V})_o\right) \right\}, \quad (11) $$

where the *density*, *velocity*, *gravity* and *observer* terms, in the first, second, third and fourth line, respectively contain all the cosmological information. On the other hand, the function $\mathcal{F}^{[i]}(f_o, z, \theta)$ contains all the astrophysical dependencies: e.g. the mass and spin distribution of the binary, the emitted GW energy spectrum, the clustering properties of GW events and the details of the GW detectors. The functions $b^{[i]}(z, \theta)$ and $b_{\text{evo}}^{[i]}(z, \theta)$ are the bias and the evolution bias of the i-th type of GW source, which specify the clustering properties of GW sources and characterise the formation of sources.

The cross-correlation analysis of the astrophysical SGWB (from sources at all redshifts along the line of sight) with galaxy number counts at a given redshift leads to a tomographic reconstruction of the redshift distribution of the sources (Cusin et al. 2018a; Mukherjee and Silk 2020; Cañas Herrera et al. 2020; Alonso et al. 2020b; Cusin et al. 2019b; Yang et al. 2020a). Subtleties about the noise and other characteristics of the detector play an important role, so that cross-checks on the detector performances are possible by means of this analysis (Bertacca et al. 2020; Alonso et al. 2020b). Overall, the cross-correlation analysis shows that the combination of galaxy surveys with the astrophysical SGWB can be a powerful probe for GW physics and multi-messenger cosmology.

### 2.6.3 Large-scale structure effects on gravitational-wave luminosity distance estimates

Here we discuss the effect of cosmological perturbations and inhomogeneities on estimates of the luminosity distance of compact object binary mergers through GWs.
It is important to account for such effects on GW propagation to obtain robust measurements for precision cosmography.

The main attempts to investigate perturbation effects on GWs involve the integrated Sachs–Wolfe effect (Laguna et al. 2010), peculiar velocities or accelerations (Kocsis et al. 2006; Bonvin et al. 2017; Tamanini et al. 2020; Mukherjee et al. 2021c; Nicolaou et al. 2020), lensing (Pyne et al. 1996; Wang et al. 1996; Takahashi and Nakamura 2003; Sathyaprakash et al. 2010; Dai et al. 2017; Takahashi 2017; Baker and Trodden 2017; Dai and Venumadhav 2017; Haris et al. 2018; Dai et al. 2018; Oguri 2018; Contigiani 2020; Ezquiaga and Zumalacárregui 2020; Ezquiaga et al. 2020b) and environmental effects (Barausse et al. 2014). Coherent peculiar velocity of the binaries at low redshift and weak gravitational lensing by intervening inhomogeneities can affect the identification of the hosts’ redshift. This introduces changes of typically a few percent (but occasionally much larger) in the flux, while not significantly affecting the redshift, and thus provides a source of noise in the $D_L(z) - z$ relation (Sathyaprakash et al. 2010; Hirata et al. 2010). Using the local wave zone approximation to define the tetrads at source position (Maggiore 2007), the corrections to the luminosity distance read (Bonvin et al. 2017; Bertacca et al. 2018)<fnref n="3" />

$$ \frac{\Delta D_L}{D_L} = \left(1 - \frac{1}{H\bar{\chi}}\right) v_{\parallel} - \int_0^{\bar{\chi}} d\tilde{\chi} \frac{(\bar{\chi} - \tilde{\chi})}{\tilde{\chi}\tilde{\chi}} \Delta_{\Omega}\Phi + \frac{1}{H\bar{\chi}}\Phi - 2\left(1 - \frac{1}{H\bar{\chi}}\right) \int_0^{\bar{\chi}} d\tilde{\chi}\, \Phi' - 2\Phi + \frac{2}{\bar{\chi}}\int_0^{\bar{\chi}} d\tilde{\chi}\, \Phi, \quad (12) $$

where prime denotes the derivative with respect to $\eta$ and $H = a'/a$, $\Delta\Omega \equiv \bar{\chi}^2\bar{\nabla}_{\perp}^2 = \bar{\chi}^2(\bar{\nabla}^2 - \bar{\partial}_{\parallel}^2 - 2\bar{\chi}^{-1}\bar{\partial}_{\parallel}) = (\cot\partial_{\theta} + \bar{\partial}_{\theta}^2 + \bar{\partial}_{\varphi}/\sin^2\theta)$. We can recognise in Eq. (12) the presence of a velocity term (the first term), followed by a lensing contribution, and the final four terms account for the Sachs–Wolfe, Integrated Sachs–Wolfe, volume and Shapiro time-delay effects.

The additional $D_L$ uncertainty due to the inclusion of perturbations is expected to peak at low-$z$ due to velocity contributions; however, velocity effects rapidly decrease and lensing takes over (Bertacca et al. 2018). Those results indicate that the amplitude of the corrections could be important for future interferometers such as LISA.

In presence of DE or modifications of gravity, the GW luminosity distance might differ from the EM signals also for large-scale fluctuations (Garoffolo et al. 2020; Dalang et al. 2020; Jana et al. 2021; Dalang et al. 2021). In particular, linearised fluctuations of the GW luminosity distance contain contributions directly proportional to the clustering of the DE field (Garoffolo et al. 2020, 2021) and by combining luminosity distance measurements from GW and supernovae sources, it is possible to uncover field inhomogeneities detecting them directly. see e.g. Van Den Broeck et al. (2010) for an analysis of weak lensing effects in the measurement of the DE EoS with LISA.

<footnote n="3">³ For simplicity, we have dropped all contributions evaluated at the observer, assuming concordance background model, and work in Poisson Gauge.</footnote>
# 3 Gravitational lensing of gravitational-wave signals

Section coordinators: D. Bacon, M. Zumalacarregui. Contributors: D. Bacon, G. Congedo, G. Cusin, J.M. Ezquiaga, S. Mukherjee, M. Zumalacarregui.

## 3.1 Introduction

Similarly to EM radiation, GWs feel the gravitational potential of both massive objects and also the LSS while traveling across the universe. This opens the possibility of both probing the distribution of structure in the cosmos and the fundamentals of the underlying gravitational interactions. LISA will offer a unique perspective since it will detect high-redshift GWs in a lower frequency band compared to ground-based detectors, increasing the lensing probabilities and the detectability of diffraction effects.

GW lensing phenomena are characterized by two properties: the convergence $κ^4$ (governing whether there is weak or strong gravitational lensing) and a dimensionless frequency $w$ (governing whether wave optics or geometric optics is relevant). The convergence, $κ$, is defined as the integral along the line of sight of the redshift-weighted second derivative of the Newtonian potential. Depending on the strength of the gravitational potential, two gravitational lensing regimes exist: weak ($κ ≪ 1$) or strong ($κ ≈ 1$). In the first case, the main observable effect is magnification or demagnification of the observed flux (or equivalently a change in the inferred luminosity distance). In the second case the effect can be both a magnification/ demagnification and also production of multiple signals with a time delay that is a function of the lens properties. We will examine both regimes in Sects. 3.2 and 3.3 below.

The low frequency and phase coherence of GWs allows the observation of wave effects. It is convenient to define the dimensionless frequency

$$w = 8\pi GM_L f = 4\pi \left( \frac{M_L}{10^8 M_{\odot}} \right) \left( \frac{f}{\text{mHz}} \right), \quad (13)$$

as well as a magnification factor $F(w) = h(w)_{\text{lensed}}/h(w)_{\text{unlensed}}$, where $M_L$ is the redshifted lens mass. GW lensing is accurately described by geometric optics only if $w ≫ 1$ (more precisely, $w\Delta\tilde{t}_j ≫ 1$ for all images). Then, the amplification factor (the ratio between the lensed and unlensed waveform) is a sum over multiple images $j$:

$$F(w) = \sum_j |\mu_j|^{1/2} \exp(iw\Delta\tilde{t}_j - i\pi n_j), \quad (14)$$

where $\mu_j$ is the magnification, $\Delta\tilde{t}_j$ is the time delay (in units of $4GM_L$) evaluated on the $j$-th image position. Here $n_j$ is the Morse phase of the image, respectively 0, 1/2, 1 for minima, saddle point and maxima of the time delay function (Schneider et al.

<footnote n="4">⁴ Because GW sources are effectively point-like, no image distortions are observable and shear influences GW observations only through its effect on the magnification. This is very different from lensing of galaxies, where shear distortions are directly observable.</footnote>
1992; Takahashi and Nakamura 2003). These image types are also known as type I, II and III respectively.

For low ($w$ considerably less than 1) or intermediate ($w \sim 1$) frequencies it is necessary to consider the wave optics regime, where the amplification factor is

$$F(w) = \frac{w}{2\pi i} \int dx^2 \exp(iw\Delta\tilde{t}(\vec{x}, \vec{b})), \qquad (15)$$

where the time delay $\Delta\tilde{t}(\vec{x}, \vec{b})$ is now a function of the (normalized) lens plane coordinate $\vec{x}$ and source location $\vec{b}$. (see Takahashi and Nakamura 2003 for details). The low frequency limit $w \ll 1$ corresponds to $F \sim 1$, or no magnification. This can be understood as a wave not being sensitive to an object whose effective size is smaller than its wavelength. The high frequency limit corresponds to the geometric optics result Eq. (14), as can be obtained from a Gaussian expansion around the images, which correspond to the extrema of the time delay $\delta\tilde{t}$ (sub-dominant corrections to geometric optics can be computed (Takahashi 2004)). We will discuss wave effects in gravitational lensing and LISA opportunities in Sect. 3.4.

## 3.2 Weak lensing

The first regime that GWs will undergo is weak lensing. Because of the relatively deep observations LISA will be able to achieve (median $z \sim 2$, reaching out much deeper depending on source types) most of the observed events (if not all) will be subjected to lensing by the large scale structure. In recent years weak lensing has become one of the primary tools in cosmology to study the distribution of DM and the nature of DE (through the evolution with redshift). Galaxy lensing surveys such as KiDS (Heymans et al. 2021) and DES (Troxel et al. 2018) are now providing tighter constraints on cosmological parameters, such as geometry parameters ($H_0$, $\Omega_m$, and $\Omega_\Lambda$) and matter clustering parameters ($\sigma_8$ and $n_s$) to a few percent level. Soon Euclid (Blanchard et al. 2020) and the Vera C. Rubin Observatory (Abell et al. 2009) will further constrain those parameters, extending the analysis to DE ($w_0$ and $w_\Lambda$) and various models beyond $\Lambda$CDM, down to the percent level. Most likely, the study of systematic errors affecting those measurements will be the topic of the next decade or so. GWs have the potential to revolutionize the field with virtually bias-free measurements of the luminosity distance. In fact, LISA will give us the first deep luminosity distance measurements that will be both accurate (as this relies on assuming GR is the correct theory) and also precise with typical errors up to a few percent depending on source type and position in the sky. At the same time, the typical root-mean-squared error due to lensing is 0.02 for $z < 1$ and ramping up quite rapidly with redshift (Cusin and Tamanini 2021). This will be a source of unwanted bias and extra noise in the Hubble diagram inference, which must be accounted for in the analysis. However, if modelling of the lensing is introduced, then this systematic error can become a new piece of useful information. For these reasons LISA will also establish the groundwork for planned second generation space-based detectors.
### 3.2.1 As a source of noise and bias for standard sirens

Magnified sources are easier to detect than de-magnified sources. This simple fact biases the distribution of lensing magnifications of an observed source sample. Lensing selection effects are often neglected when estimating the lensing-induced uncertainty on the cosmological distance measurement from high-redshift GW sources. However, when selection effects are included, flux conservation is no longer enforced, and the mean of the magnification distribution is shifted from 1 for sufficiently high-redshift sources. This introduces an irreducible (multiplicative) bias on the distance reconstruction, independent of the sample size. in Cusin and Tamanini (2021) the effect of selection bias on a population of MBBHs is examined, for different scenarios of MBH formation. It is shown that, while the effect of the bias on the distance estimator for sources with an EM counterpart is typically below the variance threshold, it becomes relevant for high redshift sources when the statistics of detectable events is large (e.g. when MBHs form at high redshift from remnants of Population III stars), and it should be taken into account in population studies.

### 3.2.2 As a probe complementary to standard sirens

Thanks to the variety and large number of sources, LISA will be able to probe a wide range of redshifts with sufficient statistical significance to constrain cosmological parameters. LISA will be able to measure luminosity distances to within a few percent with optimistic estimates pushing this limit down to one percent. Hence, weak lensing by the LSS will be expected to dominate the error budget of any luminosity distance measure for most sources with EM counterparts (e.g. MBBHs), whose median redshift is ~ 2 (Caprini and Tamanini 2016). This has to be taken into account in geometric analyses, such as the luminosity distance-redshift relationship as illustrated in the previous sections. Here we explore an avenue to extract additional cosmological information: it is well known that the weak lensing signal for background objects is correlated with foreground large scale structure. Foreground over- (under-) densities will magnify (demagnify) the observed luminosity distance to background sources by a small factor. This can be expressed in terms of the convergence via

$$D'_{L} = (1 - \kappa)D_{L}, \qquad (16)$$

where $D_L$ is the (unknown) unlensed distance, $\kappa$ is the projected 2D convergence, and $D'_L$ is the (observed) lensed distance. The observed luminosity distance now incorporates information about the convergence field, which in turn is cosmology dependent. In fact its power spectrum has cosmological information imprinted into it that can be extracted from data. It was long ago speculated that this signal could be exploited if large enough statistics (very high number densities and very precise distance measurements) were available in the future with second generation space detectors such as BBO (Cutler and Holz 2009). This approach was then applied to third-generation ground based detectors and deci-Hz space detectors to prove that the method would in principle be applicable to DE and modified gravity in a
tomographic analysis (Camera and Nishizawa 2013). More recently, it was shown that a joint analysis of lensing and distance-redshift on a combined dataset of LISA and third generation ground detectors would allow us to jointly constrain geometry parameters (such as $H_0$ and density) and clustering parameters ($\sigma_8$ and $n_s$) at the same time with percent precision (ΛCDM with curvature, nominal distance errors, and nominal number density) (Congedo and Taylor 2019; Mpetha et al. 2023). This is because a geometry probe such as distance-redshift can act as a strong prior to break degeneracies that are usually present in lensing. LISA will have good enough statistics to prove that such analysis is indeed possible, bringing in independent constraints on additional cosmological parameters such as $\sigma_8$. LISA will be able to perform a combined analysis of this type with other third-generation ground-based detectors such as ET, as well as proposed space-borne detectors operating in other bands (One such proposal is a deci-Hz detector currently under evaluation by ESA (Baker et al. 2021)). In this way LISA will provide an independent measurement of weak lensing, and expand the field of statistical cosmology to GWs in the mHz range.

## 3.3 Strong lensing

The probability that a GW is lensed depends on the redshift of the source and the distribution of lenses. The merger rate at high redshift depends on the astrophysical formation channels which can affect the event rate of well detected events and the sub-threshold lensed events observed by LIGO/Virgo (Mukherjee et al. 2021a) or, in future, by CE and ET (Piórkowska et al. 2013; Biesiada et al. 2014; Ding et al. 2015). One of the robust ways to predict the lensing event rate is by using the detection or an upper bound on the amplitude of the SGWB (Mukherjee et al. 2021b; Buscicchio et al. 2020). For high-redshift GWs, the intervening matter between the source and the observer magnifies the GW signal, thereby introducing a systematic error in the luminosity distance determination. Correcting for this weak lensing effect is relevant when inferring cosmological parameters from GW standard sirens (Hughes and Holz 2003; Hirata et al. 2010); see Sect. 3.2. But a fraction of these GW events pass close enough to the lens so that lensing produces multiple images of the same event. LISA could detect a few strongly-lensed massive BBH during its mission (Sereno et al. 2010). The detection of multiple lensing events could provide new means for cosmography with LISA (Sereno et al. 2011). Owing to its better sky localization and parameter estimation, LISA may be less susceptible to strong lensing false alarms (unlensed events that mimic lensing due to coincidental overlap of parameters) compared to the LIGO/Virgo detector network.

In the strong lensing regime the time delay between the images increases linearly with the mass of the lens, reaching delays of months for a $10^{12} M_\odot$ galaxy lens. Each of the lensed images has a different magnification, but they also acquire a fixed phase shift depending on how many times the image has crossed a lens caustic (Schneider et al. 1992); recall the $i\pi n_j$ term in Eq. (14). The origin of this phase shift can be also understood from the folding of the wavefront, which in strong-lensing configurations
produces a phase of the waveform associated with the properties of the lensing potential (Schneider et al. 1992; Dai and Venumadhav 2017).

In the case of minima of the lensing potential (type I images), there is no phase shift and the lensed waveform reproduces exactly (modulo magnification) the shape of the emitted signal. For maxima (type III), there is an overall change of the sign, $e^{-i\pi} = -1$, but the morphology of the signal remains intact. However, for saddle points (type II), the associated phase shift can distort the waveform, since the frequency-independent phase shift will introduce a frequency-dependent time delay that will differentially transform the GW signal (Ezquiaga et al. 2021). For non-precessing, quasi-circular orbits of equal mass binaries the emitted GWs are dominated by a single quadrupole radiation mode (22). In that case, the lensing phase shift accounts for a $\pi/4$ shift of the coalescence phase (Dai and Venumadhav 2017; Ezquiaga et al. 2021). Nonetheless, whenever higher modes, precession or eccentricity are relevant, the lensed waveform is modified so that it does not conform with expectation from (unlensed) templates in GR (Ezquiaga et al. 2021). This sets up the possibility that type II images near the detection threshold might be missed with standard template bank searches, but also offers the opportunity of identifying strongly lensed events with a single type II image if the SNR is high enough. Explorations of the capabilities of third generation ground-based detectors to identify type II images are performed in Wang et al. (2021b).

Precise measurement of the phase of long duration signals by LISA could be key in distinguishing different types of lensed images (Ezquiaga et al. 2021). The observation of very high SNR MBBHs and sources with very asymmetric masses, such as IMRIs and EMRIs, could give LISA a unique opportunity to identify strongly lensed events, even when detecting a single image, i.e. one term in Eq. (14), if the other images are received outside of the observing window or are too faint to be recovered.

Another very strong lensing effect can exist for LISA MBBH sources. If we observe these in EM during the inspiral, we expect self-lensing flares to occur when the two BHs are aligned with the line-of-sight (to within about an Einstein radius). In a recent pair of papers (Davelaar and Haiman 2022a, b), a binary emission model was ray-traced, and a distinct feature was found in the light curve imprinted by the BH shadow from the lensed BH. A dip occurs in the light curve when the foreground BH lenses the shadowed part of the background BH. This could make it possible to extract BH shadows that are spatially unresolvable by high-resolution very-long baseline interferometry. These shadows are another probe of the metric around the MBH.

## 3.4 Wave effects

GWs can be emitted at low frequencies ($w \lesssim 1$), allowing the observation of wave diffractive phenomena. For typical LISA sources, wave optics as in Eq. (15) needs to be considered for lenses with masses $M_L \sim 10^6 - 10^9 M_\odot$, cf. Eq. (13). In contrast, wave effects in gravitational lensing are not observed for EM sources, due to the
finite size of these sources and their high EM frequencies ($w \gg 1$), cf. Figs. 12, 13 of Oguri (2019).

Wave effects produce a frequency-dependent amplification factor $F$. At very low frequencies $w \ll 1$ no magnification occurs, $F \sim 1$, as a wave is not affected by an object smaller than its wavelength. As the frequency increases, an oscillatory pattern emerges, caused by the superposition of different images. This pattern persists for GWs even in the geometric optics limit $w \gg 1$ through the exponents in Eq. (14). We do not see this pattern for optical EM sources due to incoherent emission and lack of time resolution. Note that diffraction prevents the magnification for GWs from diverging, even for perfectly aligned, symmetric lenses with normalized impact parameter $|\vec{b}|/R_E \to 0$. An inspiralling binary sweeps through a range of frequencies before coalescence; a gravitational lens imprints the frequency-dependent amplification factor into the waveform. These effects are shown in Fig. 4.

Several studies address the detectability of wave effects by LISA. Sereno et al. (2010) estimates that a handful of strongly lensed events ($\lesssim 5$) will be observed in a 5 year mission. Because the lensing probabilities are dominated by galactic halos with $M_L \sim 10^{12} M_\odot$, wave effects are likely to be unimportant for most strongly lensed LISA sources. However, frequency-dependent diffraction effects in the weak lensing regime may be detectable for a small fraction of sources: (Gao et al. 2022) find that 0.1–1.6% of LISA’s MBBH with total mass $10^5–10^{6.5} M_\odot$ and redshift 4–10 may present detectable wave-optics. LISA provides new opportunities to observe these frequency-dependent effects and use them to probe the large-scale structure of the universe at the scale of sub-halos.

Wave effects can help identify GW signals (or groups of them) as lensed events and characterize gravitational lenses. We discuss both possibilities below.

### 3.4.1 Identification of lensed signals

Diffraction-induced features in the gravitational waveform could be used to identify the event as lensed. The main target for ground detectors is the first diffraction peak, appearing at low frequencies (Dai et al. 2018). These diffraction effects can be

<figure><img src="image_8.png" /><figcaption>**Fig. 4** GW lensing and wave effects. Left panel: Amplification $|F(w)|$ caused by a point lens in the wave optics regime. Low frequency signals $w \ll 1$ undergo no magnification. An oscillatory pattern emerges at higher frequencies. Right panel: Imprint of strong lensing on a MBBH source, as observed by LISA. The frequency dependence (onset of magnification, diffraction pattern) carries information about the gravitational lens</figcaption></figure>
detected even for impact parameters comparable to the Einstein radius $b \gtrsim R_E$, for which there is only one image in the geometric optics limit. Orbital parameters are expected to have some degeneracies with diffraction, but it is in principle possible to distinguish them. Spin precession creates a modulation that affects the amplitude more than the phase (while diffraction affects both equally). Orbital eccentricity can be distinguished by higher harmonics in GWs, which are not induced by gravitational lensing.

Wave effects may produce systematic errors in different analyses. Diffraction-induced distortions make the GW phase appear to arrive earlier than an unlensed EM signal (Takahashi 2017; Morita and Soda 2019), although this is only an apparent superluminality (Suyama 2020; Ezquiaga et al. 2020b). This effect has to be taken into account (Ezquiaga et al. 2020b) when, for example, inferring constraints on the speed of gravity from the possible pre-merger modulated EM brightness of a MBB (Haiman 2017).

The interplay between wave optics and gravitational polarizations may provide alternative means to characterize lensed GWs. This has been investigated using low-frequency corrections to geometric optics. These corrections cause the polarization plane defined in geometric optics to become blurred due to diffraction effects, which leads to the rise of apparent vector and scalar polarization modes (Cusin and Lagos 2020). Moreover, GWs in curved space propagate on different geodesics depending on their circular polarization (an effect known as the gravitational spin hall effect) (Yamamoto 2018; Andersson et al. 2021). These effects might provide means to characterize lensed GWs. As they are suppressed by factors of the GW frequency, these phenomena are more important for LISA than for ground-based detectors.

### 3.4.2 Reconstruction of large scale structure and lens properties

GW diffraction is sensitive to the power spectrum of matter density fluctuations on very small physical scales (Oguri and Takahashi 2020). Ground detectors are sensitive to this power spectrum at wavenumbers where baryonic effects are smallest, $k \sim 10^6 h \text{Mpc}^{-1}$ for $f \sim 0.1 - 1\text{Hz}$. Individual high SNR sources at high redshift, such as LISA observed MBBHs, may be sensitive to frequency-dependent modulations of the GW phase and amplitude (SNR $\gtrsim 500$ at mHz frequency, $z_S \sim 3$, cf. Fig. 12 of Oguri and Takahashi (2020)). Lensing induced diffraction can be enhanced substantially in non-standard scenarios that affect the small-scale power spectrum (e.g. PBHs).

Wave effects can then improve the chances of detection and provide additional information about the lens. Inspiralling binaries detectable by LISA sweep through a range of frequencies that might include the onset of diffraction for $w \sim 1$ and/or the characteristic diffraction pattern for $w > 1$; see Eq. (13). For the case of LISA MBBHs lensed by point masses and singular isothermal spheres, diffraction allows for a determination of lens mass and source position at the level of 10% to $\sim 0.1\%$ for $M_L \lesssim 10^7 M_\odot$ and $M_L \gtrsim 10^8 M_\odot$ for SNR $\sim 1000$ (Takahashi and Nakamura 2003). Point lenses can be identified for large impact parameters $b/R_E \sim \mathcal{O}(10)$ (rather than $\sim 1$), increasing significantly the likelihood of detection $\propto b^2$.
Diffraction is also sensitive to cuspy (or singular) matter distributions, producing an image with magnification $\mu \sim (GMf)^{\alpha-3}$ for a profile with central density $\rho \sim r^{-\alpha}$ ($0 < \alpha \le 2$) (Takahashi 2004). Frequency-dependent wave effects can be used to individually detect DM subhalos with $M \lesssim 10^7 M_{\odot}$ and even measure their mass profiles, although LISA might not provide enough events with sufficiently high SNR (Choi et al. 2021).

Investigations of microlensing (the effect of small-scale structure, such as stars, on the images produced by a large lens) have focused on LIGO-Virgo sources. Microlenses have masses such that wave diffraction effects are important at LISA frequencies, and would normally not leave an observable imprint (magnification $\sim 1$). However, they can be observed if they are present in regions of high magnification of the macrolens, such as caustics (Diego et al. 2019). Nevertheless, in most situations the macrolens model will be a sufficient description (Cheung et al. 2021). GW microlensing is also sensitive to the stellar mass function, particularly in regions of high magnification (Mishra et al. 2021).

# 4 Constraints on modified gravity theories

Section coordinator: D. Vernieri. Contributors: T. Baker, E. Belgacem, G. Calcagni, M. Crisostomi, N. Frusciante, K. Koyama, L. Lombriser, M. Maggiore, S. Mukherjee, M. Sakellariadou, I. D. Saltas, D. Vernieri, F. Vernizzi, M. Zumalacarregui.

## 4.1 Models and theories

GW astronomy offers a new window to test gravity complementing cosmological and astrophysical tests. GWs detectors can allow to measure the cosmological expansion, thus providing independent measurements of both the equation of state of dark energy and Hubble parameter. In particular, in the latter case, it can be used to help resolving the Hubble tension. Furthermore there are keys features in the propagation of GWs in theories beyond GR, such as the modification in the friction term, difference in the speed of propagation, additional polarizations, etc, which will provide unique ways for testing signatures of new physics beyond GR.

To test GR using the propagation of GWs, it is useful to compare it to other models that predict observational deviations. The simplest case that one can consider is the presence of an additional scalar interaction, e.g. a so-called scalar-tensor theory of gravity. Over the years, increasingly more general theories devoid of Ostrogradsky instabilities (Ostrogradksi 1850) have been introduced with the aim to account for the most general modification in this setup. The advent of the Galileon field theory (Nicolis et al. 2009) quickly led to its covariant generalization (Deffayet et al. 2011) and so to the reappraisal of the Horndeski theory (Horndeski 1974), i.e. the most general scalar-tensor theory with second order field equations. Relaxing this last assumption in a suitable way, has resulted in the beyond-Horndeski (Zumalacárregui and Garcia-Bellido 2014) or Gleyzes–Langlois–Piazza–Vernizzi model (Gleyzes et al. 2015a, b), and subsequently degenerate higher-order scalar-tensor theories (Langlois and Noui 2016; Crisostomi et al. 2016; Ben Achour et al. 2016).
At low energy, those relevant for cosmology, the Lagrangian of these theories schematically reads (Pirtskhalava et al. 2015)

$$ \mathcal{L} = \sum_{n=0}^{\infty} \sum_{m=0}^{3} c_{n,m}(\phi) \Lambda_2^4 \left( \frac{(\partial\phi)^2}{\Lambda_2^4} \right)^n \left( \frac{\Box\phi}{\Lambda_3^3} \right)^m + \dots, \quad (17) $$

where $c_{n,m}(\phi)$ are dimensionless coefficients that depend mildly on $\phi$, $\Lambda_2$ corresponds to the spontaneous Lorentz breaking scale, $\Lambda_3$ is related to the ultraviolet (UV) cut-off of these scalar-field theories while the ellipses refer to suitable gravitational terms necessary to keep the equations of motion of the propagating degrees of freedom of second order. Higher derivative terms with $m > 0$ have been introduced in cosmology to explain the accelerated expansion without a cosmological constant and simultaneously pass Solar System tests (Luty et al. 2003; Nicolis et al. 2009) using the so-called Vainshtein screening mechanism (Vainshtein 1972; Babichev and Deffayet 2013). In this case, the two scales must be chosen as $\Lambda_2 \sim (M_{Pl}H_0)^{1/2}$ and $\Lambda_3 \sim (M_{Pl}H_0^2)^{1/3}$, where $M_{Pl}$ is the Planck mass and $H_0$ the value of the current Hubble constant.

The above models represent low energy effective field theory (EFT) whose cut-off scale can be chosen according to the certain type of phenomenology they aim to describe, e.g. DE. An alternative path to follow is to work out cosmological phenomena in a top-down way (Calcagni 2020) starting from a fundamental theory of gravitational and matter interactions, such as quantum gravity.

Quantum gravity is a generic name applied to any theory where the gravitational interaction is consistent with the laws of quantum mechanics. Exploring the cosmological implications of these theories is important to assess how LISA can contribute to our knowledge of fundamental physics.

The corrections to cosmological dynamics and evolution of the universe arising in quantum-gravity scenarios are expected to be very small since they dominate at Planck scales $\ell_{Pl}$ or high curvature, compared to horizon scales $H_0^{-1} \sim 10^{60} \ell_{Pl}$ or low curvature. However, this may hold true for perturbative curvature corrections to GR, while other mechanisms may enhance deviations from classical gravity either by inflation of metric fluctuations to cosmic scales at early times or by non-perturbative cumulative effects on the propagation of GWs. The former mechanism will be explored in Sect. 9.3, while the latter will be developed in this section in parallel with the other models. In particular, the dimension of spacetime in quantum gravity changes with the scale probed and this can affect both the dispersion relation of GWs and the luminosity distance of GW astrophysical sources (Calcagni et al. 2019a, b).

Quantum effects can also manifest themselves as non-local terms in the quantum effective action, and these can affect the long-distance dynamics of the theory. In particular, assuming a non-local term that corresponds to a dynamical generation of a mass for the conformal mode leads to a predictive model that passes all current cosmological constraints (Maggiore 2014; Belgacem et al. 2018c, 2020) and predicts modified GW propagation. In this model the effect can be quite large, leading, at the redshifts accessible to LISA, to deviations from GR that can be as large as 80% (Belgacem et al. 2019a, 2020).
The gravitational interaction can be modified by considering a Lagrangian constructed with a general function of the scalar associated with non-metricity Q, i.e. f(Q) (Nester and Yo 1999; Beltrán Jiménez et al. 2020b; Bajardi et al. 2020). In this theory at least two additional scalar modes are introduced with respect to GR. The equivalent of GR is obtained when the Lagrangian coincides with the non-metricity scalar (Beltrán Jiménez et al. 2019). The f(Q) theory modifies the propagation equation of GWs due to an additional scale-independent friction term which is related to a time-dependent effective Planck mass (Beltrán Jiménez et al. 2018). As a consequence, the GW luminosity distance, $D_L^{GW}$, differs from the standard EM luminosity distance, $D_L^{em}$ (Frusciante 2021). For a specific model it has been shown (Frusciante 2021) that at the redshifts relevant for LISA, the deviations with respect to the standard cosmological model can be relevant depending on the value of the free parameter and that GW detectors such as LISA show a strong power in constraining it.

## 4.2 General expression for modified cosmological gravitational-wave propagation

### 4.2.1 Modified cosmological gravitational-wave propagation

The GW propagation over cosmological distances for generalised theories of gravity and other extensions of ΛCDM can be parametrised by the wave equation (Saltas et al. 2014; Gleyzes et al. 2015c; Lombriser and Taylor 2016; Nishizawa 2018; Belgacem et al. 2018a; Ezquiaga and Zumalacárregui 2018; Lombriser 2018; Belgacem et al. 2019c)

$$h_A'' + 2[1 - \delta(\eta, k)]\mathcal{H}h_A' + [c_T^2(\eta, k)k^2 + m_T^2(\eta, k)]h_A = \Pi_A, \quad (18)$$

where the primes indicate derivatives with respect to conformal time $\eta$. Here, we assume that the GW propagates far away from the source in the form of a plane wave. The modified friction ($\delta$), speed ($c_T^2$), mass ($m_T^2$) and the source ($\Pi_A$) are in general assumed to be functions of both time and wave-number $k$.<fnref n="5" /> Furthermore, parity-violating theories could also introduce a dependence of these quantities on the polarization index $A$.

As regards the friction term, within theories beyond GR next to the standard Hubble friction, one may encounter an extra contribution described through $\delta$. The latter can for instance parametrise an effective Planck mass evolution rate or the impact of extra dimensions, or the effect of other modifications of GR at the cosmological scale. The wave propagation may further be modified by a deviation in its speed $c_T$ and the mass of the graviton $m_T$. Finally, in standard GR the source term is due to the anisotropic stress tensor (a typical mechanism that can generate it is neutrino free streaming, see e.g. Section 19.5.3 of Maggiore (2018)), and in modified gravity can in general get further contributions (as is the case for instance in bimetric theories (Hassan and Rosen 2012)). In several modified gravity theories, the only

<footnote n="5">⁵ The function $\delta$ is also indicated with other names in the literature, such as $v$ and $\alpha_M$, related to $\delta$ by $v = \alpha_M = -2\delta$. Here we follow the notation of Belgacem et al. (2019c).</footnote>
modification is in the friction term and is wavenumber-independent, so that Eq. (18) reads

$$h_A'' + 2\mathcal{H}[1 - \delta(\eta)]h_A' + k^2 h_A = 0. \quad (19)$$

In that case, the effect of the modification is such that, from the measurement of the waveform of the inspiral of a coalescing binary, we do not get the standard luminosity distance obtained from EM signals (that, in this context, we will denote by $D_L^{\text{em}}$), but rather a GW luminosity distance, $D_L^{\text{GW}} := 1/h_A$ defined as the inverse of the strain amplitude, up to a constant coefficient. In the case of Eq. (19), this quantity is related to $D_L^{\text{em}}$ by (Belgacem et al. 2018a)

$$D_L^{\text{GW}}(z) = D_L^{\text{em}}(z) \exp \left[ - \int_0^z \frac{dz'}{1+z'} \delta(z') \right]. \quad (20)$$

For the purpose of comparison with the data, it is convenient to have a simple parametrization of this effect in terms of a few parameters, rather than a full function of redshift (similarly to the ($w_0, w_a$) parametrization of the DE EoS). In practice, just as for the DE EoS, a parametrization is useful only if it is quite economical, with at most two parameters, otherwise it will be difficult to extract them from the data. For modified GW propagation, a convenient choice is (Belgacem et al. 2018b)

$$\frac{D_L^{\text{GW}}(z)}{D_L^{\text{em}}(z)} = \Xi_0 + \frac{1 - \Xi_0}{(1+z)^n}, \quad (21)$$

which depends on the two parameters ($\Xi_0, n$). This parametrization reproduces the fact that, as $z \to 0$, the ratio $D_L^{\text{GW}}/D_L^{\text{em}}$ must go to one since as the redshift of the source goes to zero, there can be no effect from modified propagation. In the opposite limit of large redshifts, in Eq. (21) $D_L^{\text{GW}}/D_L^{\text{em}}$ goes to a constant value $\Xi_0$. This reproduces what happens in typical DE models, where deviations from GR only appear in the recent cosmological epoch, so that $\delta(z)$ goes to zero at large redshift and $D_L^{\text{GW}}(z)/D_L^{\text{em}}(z)$ saturates to a constant. Equation (21) indeed fits the explicit results from a large class of modified gravity models (Belgacem et al. 2019c). The modification in the luminosity distance shown in Eq. (21) can be tested from observations by combining with EM observational probes such as the BAO and CMB. By combining the BAO angular scale $\theta_{\text{BAO}}(z)$ with the measured luminosity distance from the GW sources and the sound horizon scale from CMB $r_s$, one can write a unique relation to measure the parameters $\Xi_0$ and $n$ (Mukherjee et al. 2021e):

$$D_L^{\text{GW}}(z) \theta_{\text{BAO}}(z) = \left[ \Xi_0 + \frac{1 - \Xi_0}{(1+z)^n} \right] (1+z) r_s. \quad (22)$$

If GR is the correct theory of gravity, the above equation indicates that the product of $\theta_{\text{BAO}}(z)$ and $D_L^{\text{GW}}(z)$ should scale as $(1+z)r_s$ as a function of redshift. Otherwise, it will exhibit a departure from $(1+z)r_s$ which can be explored to reconstruct the parameters $\Xi_0$ and $n$ for both bright and dark standard sirens as a function of redshift (Mukherjee et al. 2021e).
The ($Ξ_0$, *n*) parametrization is the one in terms of which have been written explicit codes and Markov Chain Monte Carlo (MCMC), in order to assess the sensitivity of LISA to modified GW propagation, including the degeneracies of $Ξ_0$ with other parameters such as $w_0$, see below. Other two-parameter parametrizations can be found in Belgacem et al. (2019c), Calcagni et al. (2019b).

### 4.2.2 The damping term (standard sirens)

Late-time constraints on the modified damping term have been forecasted in Lombriser and Taylor (2016), Belgacem et al. (2018a), Amendola et al. (2018a), Belgacem et al. (2019c), Baker and Harrison (2021), Mukherjee et al. (2021e), Frusciante (2021) based on standard siren tests (Schutz 1986; Holz and Hughes 2005), and early-time modifications can be constrained by CMB B-modes (Amendola et al. 2014). However, the dark sirens detected by LIGO and Virgo already impose some actual limits on $Ξ_0$ (Arai and Nishizawa 2018; Belgacem et al. 2018b; Lagos et al. 2019; Finke et al. 2021). The modified friction term can also be constrained with GW data only analyzing the mass distribution of compact binaries (Ezquiaga 2021), similarly to the methods developed to bound $H_0$ and other cosmological parameters (Farr et al. 2019; Ezquiaga and Holz 2021; Mastrogiovanni et al. 2021). So far there has not been much exploration of a possible frequency dependence of $\delta$. (See, however, Belgacem et al. (2019c) for forecasts on oscillations in the GW amplitude.)

In specific models of modified gravity, the parameters effecting the GW luminosity distance also impact EM observables. For the Horndeski family of theories, constraints of this type can be found in Alonso et al. (2017), Spurio Mancini et al. (2019), Bonilla et al. (2020), Allahyari et al. (2022), and the combination of EM+GW constraints in Noller and Nicola (2019), Baker and Harrison (2021).

In certain theories GWs can decay via interactions. This effect is induced by non-linear couplings from $m = 2$ in Eq. (17), as the spontaneous violation of Lorentz symmetry allows GWs to decay into two scalar excitations perturbatively (Creminelli et al. 2018) or with resonant enhancement (Creminelli et al. 2019) for which, importantly, LISA will probe a different and complementary region of the parameter space already constrained by LIGO/Virgo.

### 4.2.3 The gravitational-wave speed

The GW speed $c_T$ is constrained by a variety of measurements. The detection of ultra high energy cosmic rays implies a strong constraint on gravitational Cherenkov radiation from subluminal propagation of the waves, as otherwise the radiation would decay away at a rate proportional to the square of their energy $O(10^{11}$ GeV) before reaching us (Moore and Nelson 2001; Caves 1980; Kimura and Yamamoto 2012). For galactic $O(10$ kpc) or cosmological $O(1$ Gpc) origin, the relative deviation in $c_T$ is constrained to be smaller than $O(10^{-15})$ or $O(10^{-19})$, respectively. This bound, however, only applies for subluminal propagation, redshifts of $z \lesssim 0.1$, and modifications in the high-energy regime.
Another constraint on $c_T$ at the subpercent level can be inferred from the energy loss in binary pulsar systems (Beltran Jimenez et al. 2016; Brax et al. 2016). Constraints on $c_T$ have also been discussed as forecasts for potential arrival time comparisons with nearby supernovae emissions (Nishizawa and Nakamura 2014; Saltas et al. 2014; Lombriser and Taylor 2016; Brax et al. 2016), which are however very rare, for LISA eclipsing binary systems (Bettoni et al. 2017), or the weak bounds that can be inferred without counterpart emissions from BBH mergers (Cornish et al. 2017), from the comparison of the GW arrival times between the terrestrial detectors (Blas et al. 2016), or the CMB B-mode power spectrum (Raveri et al. 2015) for early-time modifications. A stringent and prominent direct constraint on deviations of $c_T/c = 1$ of $\lesssim \mathcal{O}(10^{-15})$ was obtained from the arrival times of the GWs from the LIGO/Virgo event GW170817 (Abbott et al. 2017d, c) and its EM counterparts. As anticipated, the measurement left a strong impact across a wide range of cosmic acceleration models (Lombriser and Taylor 2016; McManus et al. 2016; Lombriser and Lima 2017; Creminelli and Vernizzi 2017; Sakstein and Jain 2017; Ezquiaga and Zumalacárregui 2017; Baker et al. 2017). This can be observed from Eq. (17): for $\dot{\phi}^2 \sim \Lambda_2^4$ and $H\dot{\phi} \sim \Lambda_3^3$, a $c_T$ very close to c implies a fine tuning of the coefficients $c_{n,2}$ and $c_{n,3}$ for all n.

Importantly, however, the constraint only applies to low redshifts of $z \lesssim 0.01$ and the LIGO/Virgo frequency range (Battye et al. 2018; de Rham and Melville 2018). In particular, it was argued that UV completion terms for modified gravity theories naturally recover a luminal speed of gravity in the high-energy limit tested by GW170817 while allowing deviations at lower energies relevant to modifications that could drive cosmic acceleration (de Rham and Melville 2018). LISA will provide a threefold improvement over a GW170817-like bound. With the detection of GWs from massive BBHs up to $z \sim 10$ and their EM counterparts $c_T$ will be tested across much larger distances, tightening the current constraint. Additionally, a more robust measurement of arrival-time delays can be achieved (Mangiagli et al. 2020; Haiman 2017; Tang et al. 2018). Furthermore, LISA tests the frequency range below the expected UV transition of modified gravity models (de Rham and Melville 2018). Besides providing a new test of GR at different energy scales, a measurement of the speed of GWs with LISA is thus of particular relevance to cosmic acceleration models. Another way to relax the constraints on $c_T$ is by considering larger values of $\Lambda_3$ or smaller values of $\dot{\phi}^2$. For instance, it was noted that by raising $\Lambda_3$ in Eq. (17) by three orders of magnitudes is enough to agree with most of the constraints discussed in this section (Noller et al. 2020). Although this means that cosmological effects of these theories become irrelevant, their consequences could become manifest in compact astrophysical objects, potentially probed by LISA.

### 4.2.4 Dispersion

An energy $E$ dependent GW velocity is also introduced with non-vanishing graviton mass $m_g = \mu h/c^2 \neq 0$, where the group velocity becomes $(v_g/c) \simeq 1 - (m_g c^2/E)^2$. The combination of current LIGO/Virgo sources yields a constraint of $m_g \leq 1.76 \times$
$10^{-32}$ GeV/$c^2$ (Abbott et al. 2021f). More distant and more massive sources are generally more effective in constraining the dispersion relation and the lower energy range of LISA will be favourable for tightening the bounds on $m_g$. Currently, the strongest bound on the mass of the graviton with $m_g \le 7 \times 10^{-41}$ GeV/$c^2$ is inferred from weak gravitational lensing (Choudhury et al. 2004).

Frequency dependent modifications in the velocity term of Eq. (18) can also more generally be parametrised with the group velocity $(v_g/c) \approx 1 + (\alpha - 1)A_\alpha E^{\alpha-2} + \mathcal{O}(A_\alpha^2)$, where $\alpha$ and $A_\alpha$ may parametrise quantum-gravity effects, extra dimensions, or Lorentz invariance violations (Mirshekari et al. 2012). The parametrisation recovers for instance the effect of a massive graviton for $\alpha = 0$ with $A_0 = m_g^2 c^4$. In quantum gravity, the parameter $\alpha$ depends on the dimension of spacetime and of momentum space (Calcagni et al. 2019b).

### 4.2.5 Interaction terms and gravitational-wave oscillations

Finally, the presence of a source term $\Pi \ne 0$ further modifies the GW amplitude with an oscillatory correction of order $\Pi \gamma_{ij} / (Hk)$ (Nishizawa 2018). This effective source term appears naturally in theories with multiple tensor fields. These tensor fields could have a fundamental origin like in bigravity (de Rham et al. 2011; Hassan and Rosen 2012) or they could be an effective field arising from a combination of multiple vector fields as in Yang-Mills (Cervero and Jacobs 1978; Galtsov and Volkov 1991; Darian and Kunzle 1997), Abelian multi gauge fields in a gaugid configuration (Piazza et al. 2017) and multi Proca fields (Armendariz-Picon 2004; Hull et al. 2015; Allys et al. 2016; Beltran Jimenez and Heisenberg 2017). (see Beltrán Jiménez et al. 2020a for a survey of the theory landscape.) The induced GW oscillations have been studied in concrete examples such as (Max et al. 2017, 2018) and gauge field DE (Caldwell et al. 2016; Caldwell and Devulder 2019). The presence of additional tensor modes can leave an imprint in the luminosity distance, introduce waveform modulations and chiral effects (Beltrán Jiménez et al. 2020a).

The effect of GW oscillations on the luminosity distance can be constrained with LISA standard sirens (Belgacem et al. 2019c) (see Sec. 3.3). In the coherent regime, GW oscillations appear as frequency dependent modulations of the waveform and can thus be tested without the need of redshift information, cf. Fig. 11 in Belgacem et al. (2019c). GW oscillations in bigravity depend on a mixing angle (defining the interaction between the two spin-2 fields) and are suppressed by the graviton mass as $\tilde{m}_g^2/f^2$. The lower frequency range and high SNR will allow LISA to improve constraints on $m_g$ to $m_g \lesssim 2 \cdot 10^{-25}$ eV, an improvement by 3 orders of magnitude with respect to current LVK bounds ($m_g \lesssim 10^{-22}$) eV (Max et al. 2017).

### 4.2.6 Polarisation effects, gravitational-wave lensing and triple systems

GW propagation beyond GR can lead to effects on standard $(+, \times)$ or novel GW polarizations (e.g. scalar and vector waves). Due to Lorentz symmetry and parity, on the homogeneous FLRW background $h_\times, h_+$ interacts only with additional tensor
polarizations, either fundamental (bigravity, multigravity) or composite (multiple vector fields), leading to the effects discussed in Sect. 4.2.5. Theories including parity-violating interactions can lead to propagation effects on left/right polarized GWs (Alexander et al. 2008; Okounkova et al. 2022).

Interactions between standard and novel polarizations in parity-preserving theories occur beyond the homogeneous background, leading to new phenomena in GW lensing beyond GR (Ezquiaga and Zumalacárregui 2020). Locally, the propagating field is a superposition of the $+, \times$ and novel polarizations, known as propagation eigenstates. Propagation eigenstates evolve independently with a well defined speed that can depend on the position and direction. In general, each propagation eigenstate has a different speed, leading to birefringence (polarization dependent arrival time and deflection angles). Birefringence causes very distinct effects when the (total) time delay between different polarizations is shorter than the duration of the signal in the observable band. Long signals at high redshift (e.g. MBBHs) will allow LISA to probe birefringence in ways complementary to ground detectors.

GW lensing beyond GR can be particularly powerful in configurations in which the source and the lens are in close proximity. Interesting targets for LISA are stellar mass binaries in the vicinity of a MBH (a scenario suggested by the possible EM counterpart to GW190521, Graham et al. 2020). GW birefringence effects can be strongly enhanced by the strong gravity, potentially allowing tests of Horndeski theories at even higher precision than GW170817 (cf. Sect. 4.2.3 and Fig. 16 in Ezquiaga and Zumalacárregui 2020). GW emission by a binary near a MBH can be affected by the scalar hair of the central body (Brax et al. 2020). A stellar-mass binary orbiting a MBH undergoes a characteristic Doppler shift and (possibly) a strong lensing pattern that can be used to identify the triple system (Toubiana et al. 2021). Those systems can excite the quasi-normal modes of the central BH during the inspiral, with stellar-mass binaries detected by LISA sensitive to central masses $M \sim 10^6 M_\odot$ (Cardoso et al. 2021). These configurations will allow LISA to perform novel tests of GR, either through the confirmation of the triple-nature of the system or by direct detection of quasi-normal modes.

### 4.2.7 Relation to ppE framework

The modifications in Eq. (18) can also be cast into the parameterized post-Einsteinian (ppE) framework (Yunes and Pretorius 2009), introduced to parametrise the effects of GR departures in the dynamical strong-field regime on the gravitational waveforms from the binary coalescence of compact objects,

$$h(f) = \left( 1 + \sum_j \alpha_j u^j \right) e^{i \sum_k \beta_k u^k} h_{\text{GR}}(f), \qquad (23)$$

where $u \equiv (\pi Mf)^{1/3}$ with chirp mass $M$ and frequency $f$. The GR waveform $h_{\text{GR}}(f)$ is reproduced for vanishing $\alpha_j$ and $\beta_k$. The series in $\alpha_j$ can be expressed as an integral over $\delta$ whereas the series in $\beta_k$ may be expressed as an integral involving $c_T$ and $m_T$ (Nishizawa 2018). In the ppE subclass of the generalised inspiral-merger-ringdown
waveform model $h(f) = e^{i\delta\Phi_{\text{GIMR}}} h_{\text{GR}}(f)$ (Li et al. 2012) (see e.g. (Nishizawa 2018) for the details of this parameterization) the modification is hence independent of $\delta$. The above parametrization has been further extended to include parity violating effects and then used to put constraints on parity violation in gravity by using the LIGO/Virgo O1/O2 catalog (Yamada and Tanaka 2020).

## 4.3 Open problems

An interesting open problem is the mixture of propagation and source effects for extended theories of gravity and the role of screening mechanisms. A post-Newtonian expansion for the source emission in screened regimes can be performed using a scaling relation (McManus et al. 2017; Renevey et al. 2020). For some GR extensions the relevant modification in $v$ are in fact determined by the screened environments of emitter and observer rather than the cosmological background whereas for $c_T$ screening effects may safely be neglected.

In addition to the effects on the propagation of GWs, modified gravity theories can also affect the generation of GWs. Shift-symmetric scalar-tensor theories are particularly suited for describing DE since they mediate a long range (massless) interaction. However, in these theories, there is a no-hair theorem that proves that the only BH solution is locally isometric to Schwarzschild (Hui and Nicolis 2013; Creminelli et al. 2020a). One way to avoid the no-hair theorem is to consider a time dependent scalar field solution and a number of hairy BH solutions have been found in Horndeski, beyond-Horndeski as well as degenerate higher-order scalar-tensor theories (Babichev et al. 2016). Another possibility for hairy BH solutions is to involve the coupling with the Gauss-Bonnet invariant (Sotiriou and Zhou 2014). A no-hair theorem also exists for neutron stars (Barausse and Yagi 2015; Lehébel et al. 2017). Again the time dependence of the scalar field can evade this no-hair theorem (Sakstein et al. 2017; Ogawa et al. 2020). It is still an open question whether these hairy solutions can be formed dynamically and how they modify the generation of GWs such as BH ringdown (Noller et al. 2020). We note that even if modified gravity theories predict the same BH solutions in static environments, this is not necessarily the case in dynamical situations, and also their perturbations are different in different theories (Barausse and Sotiriou 2008; Tattersall et al. 2018).

An interesting possibility to evade such no-hair theorems, at least for neutron stars, is to allow for a conformal coupling with matter. The conformal coupling breaks the shift symmetry, but so long as this is Planck suppressed, this breaking will be soft. In this case non-trivial solutions featuring screening exist (Babichev et al. 2009) and, most importantly, pass the test of non-linear numerical evolution (ter Haar et al. 2021).

These modified gravity theories admit the existence of scalar GWs. A detection of scalar GWs is a smoking gun for the deviation from general relativity (Scharre and Will 2002). Screening mechanisms are expected to suppress the generation of scalar GWs (Dar et al. 2019), however, very recently, numerical relativity simulations have been performed, which point at a partial breakdown of the screening in BH collapse (ter Haar et al. 2021; Bezares et al. 2021) and in the late inspiral and merger of BNSs
(Bezares et al. 2022). In more detail, stellar collapse seems (quite surprisingly) to produce a very low frequency signal potentially detectable by LISA, while waveforms from BNSs seem to deviate from their GR counterparts at the quadrupole (but not dipole) multipole order.

Finally, there are various theoretical constraints on modified gravity models considered in this section. A notable one is the decay of GWs into scalar field perturbations (Creminelli et al. 2018, 2019). This is due to the low strong-coupling scale of the theory, which is close to the energy scale probed by LIGO/Virgo. Other constraints can be placed in terms with $m = 1$ and $m = 2$ in Eq. (17), and come from avoiding instabilities in the scalar field sector that can be induced by passing GWs (Creminelli et al. 2020b). If new states are present at a scale parametrically below the cut-off, the theory is of no use for the GWs predictions at LIGO/Virgo (de Rham and Melville 2018). LISA will detect GWs at much lower frequencies and will play an essential role in constraining these theories.

# 5 Stochastic gravitational-wave background as a probe of the early universe

Section coordinator: S. Kuroyanagi. Contributors: N. Bartolo, C. Caprini, G. Cusin, D.G. Figueroa, S. Kuroyanagi, M. Peloso, S. Renaux-Petel, A. Ricciardone.

## 5.1 Introduction

The SGWB contains information on the GW events that occurred across the whole history of the universe. Its measurement has an immense value for cosmology. It gives access to stages of the early universe that cannot be directly probed by either EM or neutrino observations. The universe before the CMB epoch was indeed not transparent to photons, and neutrinos produced in the very early universe are too weak to be detected at current or forthcoming experiments. We now introduce some key concepts about the SGWB and its observables, while we focus on the physics of its cosmological sources in the subsequent sections.

### 5.1.1 What is a stochastic gravitational-wave background?

Early universe phenomena that emit GWs typically lead to the production of a background of GWs of stochastic nature. This means that the tensor perturbation $h_{ij}(x, t)$ that defines the background is a random variable with different realisations everywhere in space. SGWBs can therefore be characterised only at a statistical level, by means of ensemble averages. As there is only one observable universe, what is customary in cosmology is to invoke the *ergodic hypothesis*, stating that either spatial or temporal averages are equivalent to an ensemble average of the underlying statistical distribution. In other words, we interpret that, by observing today a large region of the universe (or for this matter a given region for long enough time), we have access to many realisations of the system. This holds under two conditions: 1) if the universe is almost homogeneous and isotropic at the time of the GW production,
so that the initial condition for the generation of GWs is the same (in a statistical sense) everywhere; 2) if the GW source respects causality. Under these circumstances, the properties of SGWBs from the early universe can be studied by means of the ergodic hypothesis.

Due to causality, a cosmological GW source acting at a given time in the early universe cannot produce a signal correlated at length and time scales larger than the cosmological horizon at the moment of GW generation. If we denote with a subscript $p$ the time of production, the (physical) correlation scale of the emitted GWs must satisfy $\ell_p \le H_p^{-1}$.<fnref n="6" /> Equivalently, the GW signal cannot be correlated at time scales larger than $\Delta t_p \le H_p^{-1}$. At the present time we have access to much larger length and time scales than today's redshifted scale associated with $H_p^{-1}$, so a SGWB signal in our detectors is perceived as the superposition of many signals uncorrelated in time and space. The number of independent signals in a given region today can be actually counted, knowing the evolution of the universe and the details of the GW generation mechanism.

We note that the above arguments remain valid also for causal sources continuously emitting GWs over a long period, say during several Hubble times. The paradigmatic example of this is the case of a cosmic string network, which emits GWs continuously, all the way from the moment of the phase transition that created it, till today. The SGWB signal in this case is mostly dominated (at least in the Nambu Goto (NG) approach to strings, see Sect. 7) by the superposition of the GWs emitted by sub-horizon string loops. This GW signal is perceived today as a background formed by the superposition of many GW emissions at different times of cosmic history, and from many different regions. Therefore the observation of this signal today cannot be resolved beyond its stochastic nature, exactly for the same reasons discussed above. The main difference in fact between a continuously sourced background and one arising from a source localised in time, is rather that the former extends over a long frequency range, precisely because its source has been emitting during many Hubble times.

In the case of inflation, the above arguments do not apply, as the causal horizon grows exponentially during the inflationary phase. Yet, it is well known that inflation produces a SGWB, as we discuss in Sects. 8 and 9. The inflationary background of GWs is actually considered to be stochastic because of the intrinsic quantum nature of the generating process. In particular, this background originates during inflation due to the quantum vacuum fluctuations of tensor metric perturbations, which become a random variable. The tensors become effectively classical during the inflationary accelerated expansion as their wavelengths are exponentially stretched to super-Hubble scales, where they acquire very large occupation numbers. (See the discussion e.g. in Caprini and Figueroa (2018).) This transition after Hubble-crossing renders the tensor perturbation a stochastic variable. When the tensor perturbations re-enter the Hubble radius during cosmic evolution after inflation, they form a GW signal that is intrinsically stochastic.

<footnote n="6">⁶ Let us note that we have used the inverse Hubble rate $H_p^{-1}$ as the cosmological horizon, as that is a good approximation for most of the cosmological evolution of the universe, except during inflation.</footnote>
Finally, let us notice that the superposition of astrophysical GWs also forms a SGWB, which will act as a foreground for cosmological SGWBs. We will briefly review possible astrophysical background contributions expected in the mHz band in Sect. 5.2.2.

### 5.1.2 How do we characterize a stochastic gravitational-wave background?

Since we speak about random tensor variables, the Fourier modes $h_r(\mathbf{k}, t)$ are also considered to be random variables. For statistically homogeneous and isotropic, unpolarised and Gaussian GW backgrounds, we define the tensor amplitude power spectrum as

$$ \langle h_r(\mathbf{k}, t) h_p^*(\mathbf{q}, t) \rangle = \frac{8\pi^5}{k^3} \delta^{(3)}(\mathbf{k} - \mathbf{q}) \delta_{rp} h_c^2(k, t), \quad (24) $$

with $h_c$ a dimensionless real function, depending only on the time $t$ and the comoving wave-number $k = |\mathbf{k}|$. The previous factor $8\pi^5$ has been chosen so that we can write

$$ \langle h_{ij}(\mathbf{x}, t) h_{ij}(\mathbf{x}, t) \rangle = 2 \int_{0}^{+\infty} \frac{dk}{k} h_c^2(k, t), \quad (25) $$

with the factor 2 in the RHS as a convention motivated by the fact that the LHS involves the contribution from two independent polarisations. It is then clear that $h_c(k, t)$ represents the characteristic tensor amplitude per logarithmic wave-number interval and per polarisation state, at a given time $t$. It is common to denote the combination $\frac{h_c(k,t)^2}{2} \equiv P_h(k,t)$ as power spectrum. This last quantity is used in Secs. 8 and 9.

Another relevant quantity to characterise a SGWB is the spectrum of GW energy density per logarithmic wave-number interval, $d\rho_{GW}/d\log k$, defined as

$$ \rho_{GW} = \frac{\langle \dot{h}_{ij}(\mathbf{x}, t) \dot{h}_{ij}(\mathbf{x}, t) \rangle}{32\pi G} = \frac{\langle h'_{ij}(\mathbf{x}, \eta) h'_{ij}(\mathbf{x}, \eta) \rangle}{32\pi G a^2(\eta)} = \int_0^{+\infty} \frac{dk}{k} \frac{d\rho_{GW}}{d\log k}, \quad (26) $$

where in the second equality we have converted the derivatives with respect to the physical time $t$ into derivatives with respect to the conformal time $\eta$. An expression for the GW energy density power spectrum $d\rho_{GW}/d\log k$ valid for free waves inside the Hubble radius, can be found using the simple relation $h_c'(k, \eta) \approx k^2 h_c^2(k, \eta)$, where the prime symbol ($'$) represents the derivative with respect to the conformal time, and $h_c'(k, \eta)$ characterises the expectation value $\langle h_{ij}'(\mathbf{x}, \eta) h_{ij}'(\mathbf{x}, \eta) \rangle = 2 \int_0^{+\infty} \frac{dk}{k} h_c'^2(k, \eta)$. Using this, we arrive at (Caprini and Figueroa 2018)

$$ \frac{d\rho_{GW}}{d\log k} = \frac{k^2 h_c^2(k, \eta)}{16\pi G a^2(\eta)}. \quad (27) $$

As $h_c(k, \eta) \propto 1/a(\eta)$ for sub-Hubble modes, the GW energy density is diluted as
radiation with the expansion of the universe, $\rho_{GW} \propto a^{-4}$, as expected for massless degrees of freedom.

Finally, in order to connect the above expressions with experimental observables, we need to re-express the GW spectra today in terms of the present-day physical frequency, $f = k/(2\pi a_0)$, associated with the comoving wave-number $k$. In order to do this, we first write the characteristic tensor amplitude per logarithmic frequency interval today as $h_c(f) = h_c(k, t_0)$. We then define the one-sided spectral density of a SGWB as

$$S_h(f) = \frac{h_c^2(f)}{2f}, \qquad (28)$$

which has dimensions Hz$^{-1}$. The reason to define this quantity is that it is directly comparable to the noise in a detector, parametrised by $S_n(f)$. On the other hand, it is convenient to normalise the spectrum of the GW energy density per logarithmic frequency interval like

$$\Omega_{GW}(f) = \frac{1}{\rho_c} \frac{d\rho_{GW}}{d\log f}, \qquad (29)$$

where $\rho_c = 3H^2/(8\pi G)$ is the critical energy density at time $t$. The quantity traditionally considered by cosmologists is then

$$\Omega_{GW}^{(0)}(f) = \frac{4\pi^2}{3H_0^2} f^3 S_h(f), \qquad (30)$$

which corresponds to the normalised energy density spectrum today. In terms of the dimensionless amplitude $h_c = \sqrt{2f} S_h$, we can also write the following relations

$$S_h(f) = 7.98 \times 10^{-37} \left(\frac{\text{Hz}}{f}\right)^3 h^2 \Omega_{GW}(f) \frac{1}{\text{Hz}}, \qquad (31)$$

$$h_c(f) = 1.26 \times 10^{-18} \left(\frac{\text{Hz}}{f}\right) \sqrt{h^2 \Omega_{GW}(f)}. \qquad (32)$$

### 5.1.3 Current constraints by other observations

There are several ways to place constraints on the SGWB energy density. We briefly review different types of observations which provide constraints on the SGWB at different frequencies.

BBN places a bound on a primordial (i.e. prior to BBN) SGWB as the latter contributes to the total energy density of extra relativistic species and affects the expansion rate of the universe during BBN. In order not to spoil BBN by changing the resulting light-element abundances, the energy density of the SGWB should satisfy (Caprini and Figueroa 2018; Allen 1996)
$$ \int d(\ln f) \Omega_{\text{GW}}(f) h^2 \le 5.6 \times 10^{-6} (N_{\text{eff}}^{\text{(upper)}} - 3.046), \quad (33) $$

where $N_{\text{eff}}^{\text{(upper)}}$ is the upper bound on the effective number of relativistic degrees of freedom $N_{\text{eff}}$. Applying the $2\sigma$ upper limit for $N_{\text{eff}}$ from BBN (obtained by observations of $^4$He and D) as $N_{\text{eff}}^{\text{(upper)}} = 3.41$ (Cyburt et al. 2016), we obtain $\Omega_{\text{GW}}h^2 < 2.3 \times 10^{-6}$ for a logarithmic frequency bin. This constraint is valid for GWs generated before the BBN epoch $T \sim 1$ MeV.

The CMB is also a powerful tool to constrain the primordial SGWB. CMB anisotropies on large scales, both in temperature and in polarisation (E and B modes), can be induced by the tensor metric perturbations. Non-detection of such anisotropies provides a bound on the tensor-to-scalar ratio $r$ and the current best bound is given by the combined Planck and BICEP2/Keck data up to 2018: $r < 0.036$ at $k = 0.05$Mpc$^{-1}$ (Akrami et al. 2020; Ade et al. 2018; Tristram et al. 2021; Ade et al. 2021), which corresponds to $\Omega_{\text{GW}}h^2 < 1.5 \times 10^{-16}$ at $f = 7.7 \times 10^{-17}$ Hz. On small scales $f > 10^{-15}$ Hz, GW modes are well inside the cosmological horizon and behave like massless neutrinos, contributing to $N_{\text{eff}}$, thus affecting the growth of density perturbations as well as the expansion rate at recombination (Smith et al. 2006; Sendra and Smith 2012; Clarke et al. 2020). The updated constraint from the temperature anisotropy is $\Omega_{\text{GW}}h^2 < 1.7 \times 10^{-6}$ for adiabatic initial conditions and $\Omega_{\text{GW}}h^2 < 2.9 \times 10^{-7}$ for homogeneous (non-adiabatic) initial conditions at the 95% confidence level (Clarke et al. 2020). In the case of adiabatic initial conditions, GW perturbations evolve in the same way as neutrino perturbations. Most known sources of a SGWB produce an unperturbed background and we should impose homogeneous initial conditions, which assume no initial density perturbation.

Ground-based interferometer experiments, which have sensitivity to the frequency range of $10-100$Hz, are rapidly improving the upper limit on a SGWB. The latest constraint by the O3 run of LIGO/Virgo gives $\Omega_{\text{GW}} < 5.8 \times 10^{-9}$ (with $h = 0.679$) at the 95% confidence level for a flat (frequency-independent) SGWB (Abbott et al. 2021g). The limit is obtained by combining data from the earlier O1 and O2 runs and 99% of the sensitivity comes from the frequency band $20-76.6$ Hz.

Pulsar timing arrays probe low-frequency GWs at a frequency range of $10^{-9}$ and $10^{-6}$ Hz. Millisecond pulsars are known to have an extremely stable pulse frequency. GWs affect the pulse propagation and change the pulse times of arrival, thus their presence can be tested by regularly monitoring pulse frequencies. Radio telescope projects, such as NANOGrav (Alam et al. 2021), EPTA (Desvignes et al. 2016), and PPTA (Kerr et al. 2020), have been updating the upper limit on the GW amplitude, and recently reported the possible detection of a SGWB (Arzoumanian et al. 2020; Chen et al. 2021b; Goncharov et al. 2021). The strain amplitude for a $f^{-2/3}$ power-law spectrum inferred by the two experiments has central value at $h_c = 1.92 \times 10^{-15}$ and $h_c = 2.95 \times 10^{-15}$ at $f = 1$ yr$^{-1}$, respectively, which correspond to $\Omega_{\text{GW}}h^2 \simeq 2-5 \times 10^{-9}$ at $f = 3.2 \times 10^{-8}$ Hz. A confirmation of quadrupolar spatial correlations in the signal is needed to establish the detection.
A SGWB at µHz frequencies can in principle be observed in binary pulsars as well as in Earth–Moon-like or Earth–Satellites-like systems. The binary system’s orbits indeed exhibit resonant interactions with the SGWB in some regimes. Current estimates show that such binaries permit to rule out the existence of a SGWB with $\Omega_{GW} \gtrsim 6 \times 10^{-6}$ at frequencies $f \sim 1$ µHz, with improvement to $\Omega_{GW} \gtrsim 5 \times 10^{-9}$ by late 2030s (Blas and Jenkins 2022a, b).

## 5.2 Generation mechanisms

### 5.2.1 Cosmological

Among the targets of the LISA mission there is the cosmological SGWB, which includes different sources active in the early universe. Such a background is characterised by the spectral energy density $\Omega_{GW}(f)$ defined in Eq. (29) and by some peculiar features, which will be described in next sections, that can be useful in the process of characterisation and disentanglement from the astrophysical background, like chirality, non-Gaussianity, anisotropies and so on. A very well known example of a cosmological SGWB is the irreducible GW background due to quantum vacuum tensor fluctuations produced during inflation, which spans a large range of frequencies with an almost scale-invariant spectrum. For the simplest realisation of inflation, it has an amplitude that is too small to be detected by LISA. Besides this, there are several primordial mechanisms which can lead to a non-flat cosmological-SGWB frequency profile at the scales probed by GW interferometers (from $10^{-5}$ to $10^2$ Hz): from models where the inflaton is coupled with extra (gauge) fields (Sect. 8.2) (Barnaby and Peloso 2011; Cook and Sorbo 2012; Sorbo 2011; Barnaby et al. 2012; Dimastrogiovanni et al. 2017; Peloso et al. 2016; Domcke et al. 2016) to models with features in the scalar power spectrum (Sect. 8.3) (Flauger et al. 2010; Braglia et al. 2020; Fumagalli et al. 2021a), models where spacetime symmetries are broken during inflation (Sect. 8.4) (Endlich et al. 2014; Koh et al. 2013; Cannone et al. 2015a, b; Bartolo et al. 2016a; Ricciardone and Tasinato 2017; Bartolo et al. 2016a; Cannone et al. 2015b; Akhshik et al. 2014; Akhshik 2015), or scenarios where non-attractor phases characterise the universe evolution (Leach et al. 2001; Namjoo et al. 2013; Mylova et al. 2018). GWs sourced by second-order scalar fluctuations (Sect. 10.3) can further be associated with PBH formation (Sect. 10). These models are characterised by an amplitude which, still respecting the CMB bounds, have a large amplitude and a peculiar frequency shape which may enable detection by LISA. A dedicated analysis for the potential of the LISA space-based interferometer to detect the SGWB produced from different inflationary models has been performed in Bartolo et al. (2016b). This analysis has shown how LISA will be able to probe inflationary scenarios, in a complementary way to CMB experiments. Besides these, there are some post-inflationary mechanisms which can also generate GWs with a large amplitude at LISA scales: for instance several setups beyond the standard model of particle physics (BSM) exhibit a first-order phase transition (FOPT) around the TeV energy scale that peaks in the LISA frequency window (see Sect. 6). A dedicated analysis for the detection of a cosmological SGWB from
FOPTs have been done in Caprini et al. (2016), Caprini et al. (2020). Also cosmic defects can generate a cosmological SGWB which crosses the frequency window of the LISA detector. More precisely the GW signal from cosmic defects can be detected if the energy of the phase transition that created the defects is at the right scale (see Sect. 7). A recent analysis to probe the ability of LISA to measure this background, considering leading models of the string networks has been done in Auclair et al. (2020). In the most optimistic case, LISA might be able to probe cosmic strings with tensions $G\mu \gtrsim \mathcal{O}(10^{-17})$. It has been recently pointed out (Boileau et al. 2022) that, depending on different assumptions on the astrophysical background and the galactic foreground, LISA will be able to probe cosmic strings with tensions $G\mu \gtrsim \mathcal{O}(10^{-16} - 10^{-15})$.

The detection of any of these SGWBs from the early universe, would allow us to test high energy scales beyond the reach of particle colliders, like the Large Hadron Collider (LHC).

In Fig. 5 we collect GW cosmological signals expected to peak in the LISA frequency band and we compare them with the sensitivity of present and future GW detectors. To this end we use the PLS curves (Thrane and Romano 2013) designed to graphically assess the ability of a given detector to probe SGWBs: a power-law frequency spectrum crossing these curves is detected, provided it is not masked by a greater foreground.

<figure><img src="image_9.png" /><figcaption>Fig. 5 SGWB energy density $h^2\Omega_{\text{GW}}$ for different cosmological sources compared to the sensitivity of different GW detectors. As cosmological signals we have the vacuum GW contribution coming from inflation (grey dashed line) with $r = 0.044$ and $n_T = -r/8$, the signal expected in axion inflation models (cyan), the signal generated by cosmic string networks with $G\mu = 10^{-10}$ (brown), the signal generated by a FOPT with $v_w = 0.9$, $\alpha = 0.1$, $\beta/H_* = 50$, $g_* = 100$, $T_* = 200$GeV (pink) and the signal generated at second-order by the formation mechanism of PBHs with $f_{\text{PBH}} = 1$, $\sigma = 0.5$, $k_* = k_{\text{LISA}}$ (orange). For GW detectors we report the sensitivity of Planck (darker green), LITEBird (green), EPTA (blue), SKA (darker blue), LISA (red), DECIGO (purple), LIGO Design (black) and ET (darker black)</figcaption></figure>
### 5.2.2 Astrophysical

The astrophysical SGWB results from the incoherent superposition of signals emitted by numerous unresolved astrophysical sources from the onset of stellar activity until today. As any other background of radiation, the astrophysical SGWB is quantifiable through its isotropic energy density level and through the spatial angular power spectrum encoding its anisotropy. Many different astrophysical sources contribute to the astrophysical SGWB, including SOBBHs and BNSs (Abbott et al. 2016a; Regimbau et al. 2017; Mandic et al. 2016; Dvorkin et al. 2016a; Nakazato et al. 2016; Dvorkin et al. 2016b; Evangelista and Araujo 2014), merging MBBHs (Kelley et al. 2017), rotating neutron stars (Surace et al. 2016; Talukder et al. 2014; Lasky et al. 2013), stellar core collapse (Crocker et al. 2017, 2015) and population III binaries (Kowalska et al. 2012). The astrophysical information that can be extracted from the intensity and polarisation maps of the astrophysical SGWB, are the collective properties of a given population of astrophysical sources (redshift and mass distribution, local properties of galactic environment,...).

The astrophysical SGWB from BBHs is expected to be dominant in the LISA band (Dvorkin et al. 2016a) and below, and may become a source of confusion noise for other sources and cosmological background emissions. Observations with LISA will allow for the study of some aspects of BBH populations that are difficult to observe with ground-based interferometers. For example, at the mHz frequencies accessible to LISA, some of the binaries may not be fully circularised, and their residual eccentricities may provide an indication to their formation channel. In particular, binaries formed through dynamical processes in dense stellar clusters can have measurable eccentricities. These can be constrained for the subset of resolved merger, and in addition the distribution of eccentricities of the entire population may also affect the resulting astrophysical SGWB.

The detection of the BNS merger by the LIGO/Virgo network (Abbott et al. 2017d, e, 2020a) and the estimated rate of mergers in the local universe of $R = 13 - 1900 \text{ Gpc}^{-3} \text{ yr}^{-1}$ (Abbott et al. 2023) led to the conclusion that in the Hz band these sources may have a comparable contribution to the astrophysical SGWB relative to BBHs (Abbott et al. 2018b, 2021g). We may therefore expect that their contribution to the anisotropies of the astrophysical SGWB will also be important also for LISA. While it will be difficult to disentangle the relative contributions of BBHs and BNSs to the overall astrophysical SGWB, especially in view of the large modelling uncertainty in the BNS merger rates (Chruslinska et al. 2018; Giacobbo and Mapelli 2019), it is interesting to note that their host galaxies are expected to have different properties. In the isolated BH formation scenario discussed e.g. in Cusin et al. (2019b), BH masses are heavily influenced by the metallicity of their progenitor stars. Specifically, metal-poor stars retain most of their mass throughout their evolution and collapse to form heavier BHs. As a consequence, these BHs form preferentially in high-redshift and/or low-mass galaxies (Lamberts et al. 2016; Cao et al. 2018; Mapelli et al. 2018; Artale et al. 2019). In contrast, NSs can also form in metal-rich environments. In view of the different clustering properties of the host galaxy populations, BBHs and BNSs can in principle give rise to very different
anisotropic components of the astrophysical SGWB (Cusin et al. 2019b). Note that, to put this into practice, we need to understand better the time delays between formation and merger, which can change the host galaxy property if delays are long.

Finally, LISA will also allow one to study the astrophysical SGWB from other types of sources such as close white dwarf binaries (see e.g. Vecchio 2002), which may also produce anisotropies in the galactic plane (Ungarelli and Vecchio 2001; Kudoh and Taruya 2005). We refer the reader to Amaro-Seoane et al. (2023) for an in-depth discussion on the astrophysical populations leading to this and the aforementioned astrophysical SGWBs.

## 5.3 Characteristics of the stochastic gravitational-wave background

### 5.3.1 Frequency profile of the stochastic gravitational-wave background

As described above, there are many different mechanisms to generate a SGWB both of cosmological and astrophysical origin. The question that arises is, how to determine the origin of the SGWB, once it is detected. Identifying the SGWB source is a challenging but crucial task, essential to extract physical information from the SGWB detection. One of the most useful properties at this aim is the spectral shape of the SGWB. More specifically, one aims at reconstructing $\Omega_{GW}$ as a function of frequency, since the frequency profile of the SGWB depends on the generation mechanisms, thus providing a way to disentangle different SGWB sources. For example, the SGWB formed by the superposition of extragalactic BBHs and BNSs has a characteristic frequency profile of $\Omega_{GW} \propto f^{2/3}$, reflecting the frequency dependence of the inspiral phase of each event, thus enabling us to infer the origin. Also, see Kuroyanagi et al. (2018) for a list of astrophysical and cosmological sources, and examples of their frequency profiles. Reconstructing the frequency profile in detail also allows to explore the possibility that several sourcing components contribute to the SGWB within the same frequency band, leading to a complicated spectral shape arising from their superposition.

While the spectral shape of the primordial signal depends very much on the detail of the source properties, some fairly general considerations on the expected frequency profile of the cosmological SGWB are still possible, based essentially on causality (Caprini et al. 2009a). Sources active on a (conformal) time-scale $\Delta t \lesssim \mathcal{H}_p^{-1}$, where $\mathcal{H}_p^{-1}$ denotes the Hubble time when the source starts operating (we assume here that this occurs during the radiation dominated era), typically lead to peaked SGWB signals in terms of the variable $\Omega_{GW}(k)$. On (comoving) wavenumbers $k \lesssim 1/\Delta t$, in fact, the SGWB source is expected to be uncorrelated both in time and in space, since $k \lesssim 1/\Delta t \sim v/L < 1/L$ for $v$ and $L$ the typical speed and size of the anisotropic stresses. These latter are the part of the source energy momentum tensor leading to the GW generation. Since $\Omega_{GW}(k)$ is the GW energy density per logarithmic wave-number, white-noise anisotropic stresses (flat in $k$) lead to $\Omega_{GW}(k)$ increasing as $k^3$. This increase is expected to change at around $k_p \sim 1/\Delta t \sim v/L$, and the subsequent behaviour depends on the detail of the time and space structure of the source (Caprini et al. 2009a). The SGWB spectrum either peaks around $k_p$, and
decreases for $k > k_p$ with a slope depending on the source characteristics; or it takes
a shallower, but positive slope, in which case the peak occurs at a higher wave-
number. This latter can correspond to an inverse length $L^{-1}$, though other time and
length scales can also be relevant and show up in the $\Omega_{GW}(k)$ spectrum: for example,
the characteristic time-scale over which the source is coherent $\tau < \Delta t$ (as opposed to
the source duration $\Delta t$), or different length-scales, relevant in the space distribution of
the anisotropic stresses, other than $L$.

Notable peaks in the spectrum are then expected for SGWB sources characterised
by a finite, short duration $\Delta t \lesssim H_p^{-1}$. In particular, one example of sources giving rise
to peaked spectra that are relevant for LISA (since they are typically peaking in the
LISA band), are FOPTs related to the electroweak symmetry breaking, as described
in Sect. 6. On the other hand, for generation mechanisms that turn on at a given time
$H_p^{-1}$ in the radiation era, but continue to source GWs throughout the universe
evolution $\Delta t \gg H_p^{-1}$, the region of $k^3$ increase is less and less relevant, being pushed
towards the horizon today. One expects a wide frequency region over which the
signal features a slower increase with $k$ or, in some cases, is constant. The most
noteworthy example for LISA is the SGWB produced by topological defects in the
scaling regime, such as a NG cosmic strings network, which is exactly scale-invariant
in the LISA band for wide regions of the model parameter space (see e.g. Figueroa
et al. 2013 and the discussion in Sect. 7).

Another almost scale-invariant SGWB is the one generated by slow roll inflation.
This constitutes, however, an exception with respect to the cases described above.
The SGWB is in fact generated as the tensor perturbations reenter the horizon during
the radiation and matter eras, and therefore the spectral shape does not depend on the
causal evolution properties of some source anisotropic stresses, but on the
amplification of vacuum tensor perturbations during inflation. In some scenarios in
which the inflaton is coupled to an external field, the SGWB is indeed produced by
the field anisotropic stresses, and one can obtain blue-tilted spectra whose tilt
depends on the model. GWs actively generated at second order in perturbation theory
from large scalar fluctuations could also have peculiar features depending on the
inflaton dynamics. The scenarios pertaining to these categories that are relevant for
LISA are presented in detail elsewhere in this paper (cf. Sects. 8, 10.4, 10.3, and for a
review see Caprini and Figueroa (2018) and references therein).

5.3.2 Anisotropies and propagation effects

Angular anisotropies in the energy density of the SGWB can be an efficient way to
characterise its physical origin and properties. They provide a further tool to help in
disentangling a SGWB of cosmological origin from an astrophysical one, besides the
exploitation of their different frequency dependence. Angular anisotropies can be
imprinted both at the epoch of the SGWB generation and at later times, during its
propagation across cosmological perturbations. As such the anisotropies in the
SGWB can provide a new way to characterise and distinguish various generation
mechanisms of primordial SGWB and they allow one to probe the evolution of
cosmological perturbations. Because the universe is transparent to GWs for sub-
Planckian energies, the case of a cosmological SGWB represents a privileged observable to probe the physics of the early universe, and its anisotropies can preserve the memory of the initial conditions of the universe right after inflation.

We are interested in anisotropies and inhomogeneities in the energy density of the SGWB, therefore we allow the monopole Eq. (30) to be dependent on space and direction of observation

$$ \Omega_{\text{GW}} = \frac{1}{4\pi} \int d^2 \hat{n} \, \omega_{\text{GW}}(\eta, \vec{x}, q, \hat{n}), \qquad (34) $$

thus defining the energy density contrast as

$$ \delta_{\text{GW}} \equiv \frac{\delta\omega_{\text{GW}}(\eta, \vec{x}, q, \hat{n})}{\bar{\omega}_{\text{GW}}(\eta, q)} = \frac{\omega_{\text{GW}}(\eta, \vec{x}, q, \hat{n}) - \bar{\Omega}_{\text{GW}}(q, \eta)}{\bar{\Omega}_{\text{GW}}(q, \eta)}, \qquad (35) $$

where $q = 2\pi f$. Various approaches have been adopted to compute the angular anisotropies and their statistics (such as the angular power spectrum) both for the cosmological and astrophysical SGWB (Alba and Maldacena 2016; Contaldi 2017; Geller et al. 2018; Bartolo et al. 2019a, 2020b; Cusin et al. 2017, 2018b, 2019a; Pitrou et al. 2020; Bertacca et al. 2020; Cusin et al. 2018a; Jenkins et al. 2018, 2019a; Cusin et al. 2019b, 2020; Bartolo et al. 2020a; Valbusa Dall'Armi et al. 2021). A way to compute the SGWB anisotropies is to adopt a Boltzmann equation approach, similarly to CMB anisotropies (Contaldi 2017; Bartolo et al. 2019a, 2020b; Cusin et al. 2019a; Pitrou et al. 2020). In such an approach one considers the generation of high-frequency GW modes and their propagation across a background of lower frequency (large-scale) cosmological perturbations (which can be either scalar or tensor in nature). As for CMB photons, therefore, the propagating GWs become the cosmological carrier of the underlying cosmic inhomogeneities. Such an approach allows one to put in evidence at least two distinguishing features for a cosmological SGWB (Bartolo et al. 2019a, 2020b): first the anisotropies imprinted at the production epoch can be characterised by a strong frequency dependent contribution; secondly, if primordial non-Gaussianity are present in the background large-scale cosmological perturbations, then they will be left imprinted into the SGWB anisotropies. The bispectrum of the angular anisotropies of $\delta_{\text{GW}}$ turns out therefore to be a new probe of primordial non-Gaussianity, potentially measurable at interferometers, beyond the CMB and large-scale structure measurements (Bartolo et al. 2019a, 2020b). For these reasons, besides the information they provide for a SGWB from inflation, anisotropies can be a new probe for a whole series of phenomena. They can be produced at the epoch of generation of GWs from a phase transition (Geller et al. 2018; Kumar et al. 2021), and they can characterise also the SGWB which is unavoidably produced by second-order curvature perturbations in PBH formation scenarios (Bartolo et al. 2020a). Specific imprints in the SGWB anisotropies can be also generated by decoupled relativistic particles in the early universe, thus reinforcing the SGWB as a new window into the particle physics content of the universe (Valbusa Dall'Armi et al. 2021).

For a SGWB of astrophysical origin, the analytic derivation of energy density anisotropies can be found in Contaldi (2017), Cusin et al. (2017), Cusin et al.
(2018b), Cusin et al. (2019a), Pitrou et al. (2020), Bertacca et al. (2020). When adopting a Bolzmann-like description, one needs to add an emissivity term to the Vlasov equation for the graviton distribution function, accounting for the generation process at galactic scales (Cusin et al. 2019a; Pitrou et al. 2020). For extragalactic background components, the primary contribution to the energy density anisotropy comes from clustering (sources are embedded in the cosmic web), while a secondary source of anisotropy is due to line of sight effects (e.g. lensing, kinematic and volume distortion effects). Predictions for the energy density angular power spectrum have been presented in Cusin et al. (2018a), Jenkins et al. (2018), Jenkins et al. (2019a), Cusin et al. (2019b), Cusin et al. (2020), Bertacca et al. (2020) in the Hz band and in Cusin et al. (2020) in the mHz band. Anisotropies are typically suppressed by a factor $10^{-1}-10^{-2}$ with respect to the monopole, the range of variability depending on the underlying astrophysical model for star formation and collapse, and the angular power spectrum scales as $l^{-1}$ on large scales. Different physical choices for the process of BH collapse and mass distribution lead to differences up to 50% on the angular power spectrum in the mHz band, non degenerate with a global scaling (Cusin et al. 2020). With LISA it may be possible to constrain the dipole and quadrupole components of the angular power spectrum, for sufficiently high SNR detection (i.e. sufficiently high monopole) (Alonso et al. 2020a; Contaldi et al. 2020).

As shown in Sect. 2.6, where we have described how to compute the angular cross-spectrum between GW and other cosmological probes, the study of the cross-correlation of the SGWB energy density fluctuation with the LSS (e.g. galaxy distribution) is an interesting subject to examine to distinguish the origin (cosmological versus astrophysical) of a given background component. Unlike a cosmological SGWB, the extragalactic astrophysical background is expected to be highly correlated with the large-scale structure (see e.g. Cusin et al. 2017, 2018b, a; Jenkins et al. 2019a; Jenkins and Sakellariadou 2018; Jenkins et al. 2018; Cusin et al. 2019b, 2020; Jenkins and Sakellariadou 2019; Jenkins et al. 2019b; Bertacca et al. 2020; Pitrou et al. 2020; Mukherjee and Silk 2020; Alonso et al. 2020b; Adshead et al. 2021; Ricciardone et al. 2021; Braglia and Kuroyanagi 2021; Valbusa Dall'Armi et al. 2022). Ways to exploit this feature are discussed in Sect. 3.

### 5.3.3 Polarisations

As any background of radiation, a SGWB is fully characterised in terms of Stokes parameters, intensity (proportional to the background energy density), and *Q*, *U*, *V* parameters describing polarisation. Classical diffusion of GW radiation from massive objects can generate a net polarisation out of an unpolarised flux, playing a role analogue to Thomson scattering for CMB photons (Cusin et al. 2019a). The amount of polarisation that can be generated depends on the GW frequency, and it is more effective for large wavelength modes, for which wave effects are expected to be more important in an astrophysical context. An order of magnitude estimate of the effect gives that, in the mHz band, the net amount of polarisation generated by diffusion is suppressed by several orders of magnitude with respect to anisotropies in the intensity (Cusin et al. 2019a). As polarisation cannot be effectively generated during
propagation and astrophysical background components are expected to be statisti-
cally unpolarised at emission, the detection of a highly polarised background
component is a smoking gun of its cosmological origin.

As we review in Sect. 8.7.1, several inflationary mechanisms have been proposed that could produce a net circularly polarised SWGB, which is characterised by Stokes V parameter. A net chiral polarisation can be measured with a network of ground-based (Seto and Taruya 2007, 2008; Crowder et al. 2013; Smith and Caldwell 2017) or space-based (Orlando et al. 2021) interferometers.<fnref n="7" />

The measurement is more problematic in the case of a single planar instrument such as LISA. In this case, a left-handed GW with wave-vector $\vec{k}$ produces the same effects as a left-handed GW with wave-vector $\vec{k}_p$, where $\vec{k}_p$ has been obtained from $\vec{k}$ with a reflection on the plane of the detector. Therefore a difference between the two polarisations cannot be detected in the case of an isotropic SGWB.

A net polarisation can however be detected also by a planar instrument if the SGWB is not isotropic. As discussed in the previous subsection, It is natural to expect that a SGWB of cosmological origin has a dominant monopole component, with large-scale anisotropies of magnitude comparable to that of the CMB ones. This statement is, however, frame-dependent, and the most natural expectation is that the SGWB is isotropic in the CMB rest-frame. As seen in the CMB, the motion of the Solar System in this frame, with a velocity $v \simeq 10^{-3}$, produces a dipole anisotropy, with an amplitude suppressed by a factor $v$ with respect to that of the monopole. Seto (2006), Domcke et al. (2020b) studied how the dipole signal might allow one to measure a net chirality with LISA. This can be done through the cross-correlation between the A and E channels, which vanishes both in the case of isotropic and of unpolarised SGWB. As estimated in Domcke et al. (2020b), the SNR associated with this measurement is

$$ \text{SNR} \approx \frac{v}{10^{-3}} \left| \frac{\sum_{\lambda} \lambda \Omega_{\text{GW}}^{\lambda} h^2}{1.4 \cdot 10^{-11}} \right| \sqrt{\frac{T}{3 \text{ years}}} \quad (36) $$

where $\lambda = \pm 1$ refers to the right and left chirality, respectively, and where T is the observation time.

### 5.3.4 Non-Gaussianity

There are two types of non-Gaussianity discussed in the context of GW observations. One is the non-Gaussianity of inhomogeneities, which is defined in position or momentum space (see e.g. Sect. 8.7.3). GWs generated at sub-horizon scales cannot produce correlation across the horizon due to causality, thus the SGWB is Gaussian. Non-Gaussianity typically appears in GWs generated in the context of inflation, which could produce non-trivial spatial correlations stretched over the horizon (see Sects. 8 and 10.2.3). See Bartolo et al. (2019a, 2020b) for a detailed derivation of the non-Gaussianity expected in the SGWB, which is simply generated by the evolution

<footnote n="7">⁷ The detection of a circularly polarised SWGB at CMB scales was studied in Gluscevic and Kamionkowski (2010), Smith and Caldwell (2017), Gerbino et al. (2016), Thorne et al. (2018).</footnote>
through the background large-scale underlying inhomogeneities, similarly to what happens for CMB photons. This is computed through the angular bispectra (i.e. the three point function) of the graviton energy density.

The other type of non-Gaussianity is the one in the time signal (sometimes referred to as a SGWB in the “popcorn” or “shot noise” regime), which could be a useful statistical measure for a SGWB formed by overlapped short-duration events, such as the astrophysical background. If GW events are not frequent enough to overlap in time, the observed strain has a non-Gaussian distribution. Among the cosmological sources, the SGWB from cosmic strings could show this non-Gaussian feature (Regimbau et al. 2012) (see Method II in Sect. 7.3). For non-Gaussianity of astrophysical sources, see Sect. 12.2.

# 6 First-order phase transitions

Section coordinators: J. Kozaczuk, M. Lewicki. Contributors: M. Besancon, C. Caprini, D. Croon, D. Cutting, G. Dorsch, O. Gould, R. Jinno, T. Konstandin, J. Kozaczuk, M. Lewicki, E. Madge, G. Nardini, J.M. No, A. Roper Pol, P. Schwaller, G. Servant, P. Simakachorn.

## 6.1 Introduction

Cosmological FOPTs are one of the most attractive sources of GWs in the early universe (Witten 1984; Hogan 1986). A FOPT can occur when the Higgs or any other scalar fields are trapped in a metastable vacuum in the early universe. As the universe cools down, thermal or quantum fluctuations drive the field over or through the potential barrier, resulting in bubbles of the stable phase nucleating in the sea of metastable phase. These bubbles then expand and collide with each other to complete the transition. The collision of the bubbles and the fluid motion around them produce a SGWB. Since the generated GWs propagate to the present nearly without interaction, they maintain information about the processes that produced them at the time of generation. Their detection may therefore reveal some of the properties of the high-energy universe.

In the standard model of particle physics (SM) there are in principle two phase transitions, at two energy scales: the scale of the electroweak gauge symmetry breaking, and that of the chiral symmetry breaking in quantum chromodynamics (QCD). Given the coupling constants and the mass of the Higgs boson in the SM, it is known that the electroweak symmetry breaking is a crossover (Kajantie et al. 1996b; Csikor et al. 1999). Moreover, at almost zero quark chemical potential, the QCD phase transition cannot be first-order (Stephanov 2006). However, FOPTs are predicted in a number of extensions of the SM aimed at addressing open questions such as the origin of the observed baryon asymmetry (Kuzmin et al. 1985; Rubakov and Shaposhnikov 1996), the nature of DM, and the hierarchy problem.

The millhertz-band GWs that LISA will measure are in the best frequency band to explore FOPTs occurring between the electroweak and the multi-TeV scales. Therefore, it is important to understand how FOPTs source GWs, how such
transitions arise from microphysics, and consequently to appreciate the implications that the detection of GWs from FOPTs has for particle physics. We provide an up-to-date understanding on these topics in the following.

## 6.2 Determination of the relevant parameters

The main features of the SGWB spectrum can be determined from knowledge of four parameters related to the FOPT: the temperature $T_*$ of the plasma when the bubbles percolate, the transition strength $\alpha$, the inverse duration of the phase transition $\beta$, and the wall velocity $v_w$. These macroscopic parameters crucially depend on the underlying microscopic model, linking particle physics and cosmology. We dedicate this section to a discussion of how to compute the parameters entering the SGWB spectrum from the underlying microphysics.

The temperature at which the FOPT occurs is determined by the interplay between the nucleation of bubbles of the true-vacuum phase and the expansion of the universe. The nucleation rate per unit volume is given by $\Gamma = A(t) \exp[-S_c(t)]$, where $S_c$ is the critical bubble action. For thermally-induced nucleation, we have $A \sim T^4$ and $S_c = S_3/T$ with $S_3$ being the three-dimensional Euclidean action of the $O(3)$ symmetric bounce solution.<fnref n="8" /> The onset of the FOPT then occurs at the nucleation temperature $T_n$, at which on average one bubble is nucleated per horizon volume. It is roughly given by the temperature at which $S_c \sim 140$. The characteristic temperature for GW production is however the temperature $T_*$ at which the bubbles percolate, when approximately one third of the comoving volume has transitioned to the true vacuum. For moderately strong transitions, these temperatures are typically sufficiently close to take $T_* \sim T_n$. More precise formulas for the determination of the nucleation and percolation temperatures in from the critical action can be found in Caprini et al. (2020).

The transition strength parameter $\alpha$ describes the amount of energy released as a fraction of the radiation energy $\rho_{rad}$. Different definitions are used throughout the literature, corresponding to different ways of mapping the particle physics model to the bag EoS commonly used for determining the efficiency of converting the released energy into bulk motion of the fluid (Espinosa et al. 2010). A simple but reasonably accurate definition of $\alpha$ in terms of a given particle physics model is obtained from the difference of the trace $\theta$ of the energy-momentum tensor between the phases,

$$ \alpha = \frac{\Delta\theta}{\rho_{rad}}, \qquad \Delta\theta = \Delta V(T) - \frac{T}{4} \frac{\partial\Delta V(T)}{\partial T}, \qquad (37) $$

where $V(T)$ is the thermal effective potential. See Giese et al. (2020, 2021) for a comparison of other conventions regarding the transition strength as well as generalisations beyond the bag model. The efficiency factors for the conversion of the energy released in the transition into bulk motion can then be calculated from the transition strength $\alpha$ and the wall velocity (Espinosa et al. 2010; Giese et al. 2021).

<footnote n="8">⁸ For quantum tunnelling, $S_c = S_4$ is the Euclidean action of the $O(4)$ symmetric bounce in four dimensions.</footnote>