Probing anisotropies of the Stochastic Gravitational Wave Background with LISA

Introduction

One of the main targets of the LISA gravitational wave (GW) detector is the detection of a stochastic gravitational wave background (SGWB), which can shed light on the physics of the early universe and on astrophysical population properties not accessible with resolved sources. There are many possible astrophysical and cosmological sources which contribute to the stochastic background (see e.g., [1][2][3] for recent reviews), and up to now we have only upper bounds on its amplitude in [4], and on parameters characterising its directional properties [5,6], by the LIGO/Virgo collaboration. On the other hand we have a recent JCAP11(2022)009 claim of a possible detection of a SGWB signal in the nano-Hertz regime by the NANOGrav collaboration [7]. Typically it is expected that each source is characterized by a specific spectral shape [8,9], however, given the plethora of sources (both resolved and unresolved) which are present in the LISA band (i.e., milli-Hertz regime), it is important to study other features which can allow for a better characterization and detection of this signal.

Interesting properties can be extracted by measuring the anisotropies in the SGWB. The first attempts for extracting information on the anisotropies of the SGWB have been first done in [10] for the case of ground-based interferometers, in [11] for space-based interferometers, and in [12,13] for pulsar timing arrays.

The aim of this work, developed within the LISA Cosmology Working Group, is to analyze the capabilities of LISA [14] to detect anisotropies of the SGWB in the milli-Hertz band, making use of current instrument specifications, as well as of the latest theoretical characterizations of sources of SGWB anisotropies. The work is developed in two main parts: the first part is more theoretical in nature, and reviews our present understanding of cosmological and astrophysical sources for the SGWB and the properties of its anisotropies; the second part contains new results on the characterisation of angular response functions for LISA, accompanied by forecasts of the detectability of an anisotropic SGWB with LISA.

The theory part of our work starts with a review of a Boltzmann equation approach for analyzing anisotropies of the SGWB [15][16][17][18][19], similarly to what is commonly done for the Cosmic Microwave Background (CMB). This method is convenient for distinguishing e ects on anisotropies sourced at the moment of GW production, from anisotropies developed as GWs propagate through our inhomogeneous universe. We then discuss early universe sources of the SGWB, and we describe SGWB anisotropies produced from inflationary mechanisms and from the formation of primordial black holes (PBH). Similarly to CMB photons, gravitons are also a ected by the Sachs-Wolfe and Integrated Sachs-Wolfe e ects, both related to the propagation of GW through a perturbed universe. Besides these contributions, we discuss the intrinsic SGWB anisotropy generated at the moment of production, whose frequency-dependence represents a peculiar signature of GW. We then discuss a case where GW anisotropies are induced by primordial non-Gaussianity, showing that in certain scenarios such a contribution can be relatively large. Then we study the anisotropies expected in some post-inflationary mechanisms, like preheating in a scale invariant model [20,21] -even if the GW background in this model is typically peaked at larger frequencies [22], beyond the LISA frequency band. We finally review another two main GW sources that are characterized by anisotropies: phase transitions and topological defects. For phase transitions, if only cosmological adiabatic perturbations are considered, the fluctuations in any causally produced GW background will follow those in the CMB, and hence they are expected to be small [23,24]. For GW sourced by topological defects, anisotropies induced by a network of Nambu-Goto cosmic string loops have been computed in [25][26][27]. It has been shown that while the angular power spectrum C the quantity characterising the multipole decomposition of the SGWB spectrum -depends on the model of the loop network, the anisotropies are driven by local Poisson fluctuations in the number of loops, and the resulting angular power spectrum is spectrally white (i.e., C ¸= constant with respect to ¸), regardless of the particular loop distribution [25].

We then present the case of anisotropies generated from astrophysical sources of GWs. LISA will be sensitive to several astrophysical sources such as Super Massive Binary Black Holes (SMBBHs) with masses ≥ 10 4 -10 7 M § , stellar origin Binary Black Holes (SOBBHs), Extreme Mass Ratio Inspirals (EMRIs) and Galactic white dwarf Binaries (GBs). Beyond these resolvable sources, measurements by LISA will also be a ected by a huge number of JCAP11(2022)009 unresolvable events which will sum up incoherently, forming a SGWB [1,28,29]. At least two SGWB components are guaranteed to be present in the LISA band: a contribution due mostly due to GB inspirals in the low-frequency band of (up to ≥ 10 ≠3 Hz), and a contribution from extra-galactic BBH mergers expected at slightly higher frequencies (≥ 10 ≠3 -10 ≠2 Hz). The analytic derivation of the energy density anisotropies for an SGWB has been well studied in the literature [15,18,19,[30][31][32]. Predictions for the energy density angular power spectrum have been presented in [32][33][34][35][36][37] in the Hz band and in [37] in the mHz band (see [38] for a recent numerical code to estimate the angular spectrum of the anisotropies of the astrophysical GWB). Anisotropies show a range of variability depending on the underlying astrophysical model for star formation, mass distribution and collapse, and on the considered cosmological perturbation e ects. Due to its stochastic nature, we characterise the anisotropies in terms of their angular power spectrum taking into account all the cosmological and astrophysical dependencies.

We then move to the second part of this paper containing original results on prospects of detection of anisotropies of the SGWB with LISA, given the current instrument specifications. The characterisation of the angular resolution of space-based detectors as LISA has been pioneered in [39][40][41], and previous studies on LISA capabilities in detecting and characterising SGWB anisotropies include [42][43][44][45][46]. We start our analysis computing the angular response functions of LISA to the di erent multipoles for a statistically isotropic SGWB. We work in the A, E, T Time-Delay-Interferometry basis (see [47][48][49][50][51][52], as well as the comprehensive review [53]) and we compute the angular response as a function of frequency for the autocorrelation channels (AA, EE, and TT) and cross-correlated ones (i.e. AE, AT). We also give their analytic expression in the low frequency limit. We develop an estimator for the angular power spectrum C ¸, giving a simple analytic tool to estimate the total sensitivity of LISA to an anisotropic signal. With these tools we estimate the minimal amplitude of GW energy density needed for detecting higher multipoles. As a concrete example, we analyse the case of the kinematic dipole and quadrupole generated by boosting with a factor -© v/c an isotropic SGWB. We find that for one year of observation, -GW ≥ 2 ◊ 10 ≠11 is required to observe a dipolar signal with LISA. We also find that the detector response to the quadrupole has a factor ≥ 10 3 -relative to that of the dipole.

We then perform a Fisher matrix analysis aimed at forecasting the amplitude required on the lowest multipoles of the SGWB angular power spectrum for being detectable with LISA, given the current information on LISA strain and angular resolution sensitivity. We consider a power-law SGWB spectrum peaking at some multiple ¸characterised by a fiducial amplitude and spectral tilt.

The peak in sensitivity for ¸= 1 occurs at higher frequencies than that for ¸= 0, 2. Therefore, if we choose the pivot scale of the power-law signal to coincide with the peak sensitivity frequency of the ¸= 0, 2 multipoles, so that their detectability is only weakly sensitive to the spectral index, we then find a greater sensitivity for ¸= 1 in the case of a positive spectral tilt.

Finally we apply the maximum likelihood map-making method for stochastic backgrounds proposed in [54] to the LISA detector, in order to provide estimates for the noise angular power spectrum N ¸. We simulate and map the noise directly in the sky domain, and we take into account the specific scan strategy of LISA, which describes how the sky signal is sampled as a function of time.

The structure of this paper is as follows: in section 2 and 3 we review the main cosmological and astrophysical GW sources and their predicted angular power spectra; in JCAP11(2022)009 section 4 we present the LISA angular response function to di erent multipoles and the Signal-to-Noise (SNR) estimator for anisotropic signals. In section 5 we perform a Fisher matrix analysis for the amplitude and spectral tilt of a SGWB signal characterized by a power-law behaviour. Finally in section 6 we compute the noise angular power spectrum of LISA for di erent multipoles using a map-making approach. A conclusion and some technical appendices conclude the work.

Theoretical framework

The SGWB energy is controlled by the energy density spectrum GW defined as

with dfl GW being the energy density in GW contained in the comoving momentum interval q to q + dq, and fl c,0 corresponding to the critical energy density of the present-day universe. As we are going to discuss, we expect that the quantity GW is characterized by an averaged isotropic component plus a direction-dependent component. Both the isotropic and the anisotropic contributions are two key observables that can be targeted by the GW LISA detector. Several cosmological sources can produce a monopole GW energy density within the reach of the LISA sensitivity: inflationary models beyond vanilla single-field scenarios, where the inflaton is coupled with extra (gauge) fields [55][56][57][58][59][60][61] to models with features in the scalar power spectrum [62][63][64], or models where space-time symmetries are broken during inflation [65][66][67][68][69][70][71][72], or scenarios where non-attractor phases characterize the Universe evolution, [73][74][75], or second-order scalar induced GWs which are also responsible for PBH formation [76][77][78][79][80][81][82][83][84][85][86][87][88][89][90][91][92][93][94][95]. Also post-inflationary mechanisms can generate GW signals within the reach of the LISA detector: expected signals come from first order phase transitions beyond the Standard Model of particle physics, and from the subsequent generation of topological defects, including the irreducible SGWB from any network of cosmic defects. Forecasts about the detection of the isotropic monopole contribution have been performed in previous publications of the LISA Cosmology Working group: for inflationary scenarios in [96], for phase transitions in [97,98], and for cosmic strings in [99].

All such backgrounds are also expected to display anisotropies (direction dependence) in the GW energy density GW (f, n), which can be generated either at the time of their production [20,21,23,[100][101][102][103] or during their propagation in our perturbed universe [15][16][17]104]. For this reason anisotropies in the SGWB energy density can be considered as a new tool to characterize and distinguish various generation mechanisms of primordial GW. At the same time, they can be considered as a probe of the evolution of cosmological perturbations.

As shown in [15][16][17]105], SGWB anisotropies show strong analogies with those of the Cosmic Microwave Background (CMB), at least in the geometrical optics limit [106][107][108][109]. For this reason they can be treated using the Boltzmann equation approach, i.e., computing and evolving the equation for the gravitons distribution function f in a perturbed FLRW background, analogously to what is done for CMB photons. At zeroth order in the perturbations, the isotropy and homogeneity of the background imply that the graviton distribution depends only on time and on their frequency. The gravitons propagate freely, and their physical momentum redshifts during the propagation, as CMB photons. There is however a marked di erence between the graviton and the photon distribution, namely the JCAP11(2022)009 initial population of gravitons is not expected to be thermal, as we have in mind specific production mechanisms, such as inflation [55,56], phase transitions [23], or enhanced density perturbations leading to primordial black holes (PBH) [84,85,101], occurring at energy densities much smaller to what be required for the thermalization of the produced gravitons. This induces a sort of 'memory' of the initial state in the distribution.

The production mechanism could occur inhomogeneously in the observed universe, in a way that correlates to the large scale perturbations. This would result in an anisotropic signal arriving on Earth. Besides this initial condition, an additional anisotropic contribution is induced by the GW propagation in our perturbed universe. Working at the linearized level in a regime of a large hierarchy q ∫ k between the GW (comoving) momentum q and the (comoving) momentum k of the large scale perturbations, the graviton propagation is a ected by a Sachs-Wolfe (SW) e ect, which is dominating on large scales, and by an Integrated Sachs-Wolfe (ISW), similarly to CMB photons. An important di erence with respect to the CMB photons is associated with the 'decoupling' time of the two species: while the CMB temperature anisotropies are generated only at the last scattering surface, or afterward, the universe is instead transparent to GWs at all energies below the Planck scale. For this reason, the SGWB provides a snapshot of the universe right after inflation, and its anisotropies retain precious information about the primordial cosmological evolution.

The Boltzmann equation for the graviton distribution function f (x µ , p µ ), with x µ the graviton position and p µ = dx µ /d⁄ its momentum, is given by

where L © d/d⁄ is the Liouville operator, while C and I account, respectively, for the collision of GWs along their path, and for their emissivity from cosmological and astrophysical sources [15]. In the case of a cosmological SGWB, the emissivity term can be treated as an initial condition on the GW distribution, while, as we will see in section 3, in the case of an astrophysical background it is related to the astrophysical process that generate the GW signal at various redshifts, such as the black hole merging. On the other hand, we disregard the GW collision term since it a ects the distribution at higher orders in an expansion series in the gravitational strength 1/M Pl , where M Pl is the Planck mass. We assume that our universe is well described by a perturbed FLRW metric

where a(÷) is the scale factor as a function of the conformal time ÷, and scalar fluctuations, and h ij the transverse-traceless tensor fluctuations. We can then solve the Boltzmann equation (2.2), at background and linear levels. The background Boltzmann equation simply reads ˆf /ˆ÷ = 0, and it is solved by any distribution that is function only of the comoving momentum q, namely f = f (q). This implies that the physical momentum of the individual gravitons redshifts proportionally to 1/a. At linearized level, the evolution equation for f becomes [15][16][17] ˆf ˆ÷

where ni = qi is the direction of motion of the gravitons. The distribution function f is related to the GW energy density by

JCAP11(2022)009

where we use the spectral energy density GW introduced in eq. (2.1), which depends also on the position x where the energy density is evaluated. The su x 0 indicates a quantity evaluated today. We can account for a possibly anisotropic dependence by defining the quantity Ê GW through

and then the bar quantity ¯ GW (÷ 0 , q) is defined as spatial average (over the evaluation point x) of the above quantity GW (÷ 0 , x, q). With these ingredients we can introduce the density contrast

where the homogeneous and isotropic fractional energy density is obtained from the zeroth order distributions function f . We decompose, as for the CMB, the density contrast in spherical harmonics,

and, under the assumption of statistical isotropy, we define the multipole coe cients through

As shown in [15][16][17], it is useful to re-define the graviton distribution function as "f © ≠q ˆf ˆq (÷, x, q, n) , to simplify the first order Boltzmann equation, that now in Fourier space reads 1

where the terms on the right hand side (r.h.s.) define the so-called source function S(÷, k, n), prime denotes a derivative with respect to conformal time, and µ is the cosine of the angle between k and n. The GW density contrast is related to the and to the background energy density fractional contribution ¯ GW [16,17],

where we recall that q = qn is the graviton comoving momentum. Many of the cosmological GW scenarios mentioned above have a GW spectrum well described by a simple power law in frequency (i.e., ¯ GW Ã q n T ). In these cases the previous relation reduces to

, where n T is the tensor spectral index.

The solution of the eq. (2.10) can be decomposed as

where I, S, and T stand for Initial, Scalar and Tensor sourced terms respectively. The scalar and tensor terms correspond to the induced anisotropies arising from the propagation of 1

In the CMB case CMB = " T/T .

JCAP11(2022)009

considering perturbations around a homogeneous and isotropic solution, it becomes clear that tensor fluctuations in the gravity sector, i.e. gravitational waves, are a universal prediction of inflation.

The existence of a cosmological stochastic gravitational wave background (SGWB) can be tested across a wide range of scales, from its e ects on the CMB B-mode polarisation, all the way to direct detection via laser interferometers. In what follows, we shall focus on the latter possibility and clarify how anisotropies in the GW energy density, imprinted at the epoch of the SGWB generation, may directly probe inflationary dynamics.

These anisotropies, encoded in the first contribution I (÷ in , k, q) in eq. (2.8), carry the imprints from the initial conditions because the Universe is essentially transparent to GWs. This is to be compared to CMB photons for which anisotropies at the initial epoch "÷ in " are erased by the multiple collisions photons su er prior to the recombination epoch. We stress that, interestingly, anisotropies due to initial condition are strongly model dependent and thus provide the opportunity to test and distinguish among di erent inflationary models. To give one example, in the case of single-field adiabatic initial conditions (and for scale-invariant primordial gravitational waves) one would get, in the language of eq. (2.8):

where (÷ in , k) is the gravitational potential perturbation (in Poisson gauge), see [16,17,105,110].

Anisotropies provide a precious handle on the particle content of the very early Universe. We would like now to single out the two necessary conditions underlying the e ectiveness of anisotropies specifically as a probe of inflationary interactions: (i) naturally, a primordial GW spectrum amplitude at small scales that is well-above the sensitivity curve of laser interferometers such as LISA; (iia) a su ciently sizeable long-short mode coupling (i.e. squeezed primordial non-Gaussianity) [16,17,101,111,112], or (iib) an anisotropic background tout court [16,17].

Each of the property in (i) and (ii) are unlikely to characterize single-field slow-roll (SFSR) models of inflation. Indeed, the typical frequency profile of SFSR realisations is that of a slightly red-tilted GW spectrum, with a signal below the LISA sensitivity threshold. 2 Non-Gaussianities associated to the same SFSR paradigm are also small. Remarkably, there is a growing literature on multi-field inflationary realisations that comply with both requirements. Interesting examples of anisotropies induced by primordial non-Gaussianities include those occurring in models with light spin-2 field(s) during inflation [114,115] and set-ups with a nonstandard symmetry breaking patterns (see e.g. [65,116,117]). For examples of anisotropies engendered by an anisotropic background we refer the reader to [16,17], where the case of GWs sourced by gauge fields in axion inflation is discussed. This set-up leads to anisotropies with a significant frequency dependence, in contradistinction to what happens for CMB photons.

A general treatment of anisotropies from initial (i.e. inflationary) conditions is made possible by the Boltzmann equation and the theoretical framework expounded in section 2.1. In the remainder of this subsection, we shall describe and highlight the importance of anisotropies as a probe of primordial non-Gaussianities in the sense of (iia) defined above. We will put aside (iib) as well as assume, and later quantify, a su ciently large primordial bispectrum so as to render the anisotropy via long-short mode coupling the leading contribution. It is 2 Noteworthy exceptions include models where an attractor phase is preceded by non-attractor evolution, see e.g. [113].

JCAP11(2022)009

convenient, before elaborating on the explicit form of non-Gaussianities-induced anisotropies, to make contact with the form they take in the context of the Boltzmann treatment. The e ect of a squeezed scalar-tensor-tensor (STT) primordial bispectrum on GW anisotropies is captured by the "initial conditions" term I in eq. (2.8) via:

where k L underscores the specific bispectrum configuration (squeezed) under scrutiny and F NL is a placeholder for primordial non-Gaussianity of the STT type. An analogous relation exists for anisotropies induced by TTT-type correlators, i.e. GW non-Gaussianities.

The anisotropies of the GW energy density induced by, respectively, squeezed STT and TTT non-Gaussianity, have the following form [102,118,119]:

and

where d = ÷ 0 ≠ ÷ in is the elapsed from horizon re-entry to the present for the mode q, and the non-linearity parameters have been defined as

and the bispectra B sq are understood as defined in standard fashion from the squeezed limit of the three-point function in Fourier space.

The bispectrum component that appears in eqs. (2.18)-(2.19) is the leading physical contribution to the three-point functions. It is often the case that those bispectrum diagrams that include interactions mediated by additional (w.r.t. the single-field slow-roll case) fields give the largest contribution in terms of non-Gaussianities, squeezed or otherwise.

In order to identify the regime where non-Gaussianities provide the leading contribution to anisotropies, it su ces to report here that, schematically:

where A S is the amplitude of the primordial scalar power spectrum and r is the tensor-to-scalar ratio. The regimes of interest are then, respectively, those where the conditions F STT, sq NL ∫ 1 and F TTT, sq NL Ô r ∫ 1 hold true. It is instructive to recall, for illustrative purposes, the analytical approximation to the angular power spectrum of STT-induced anisotropies:

which has been obtained under the simplifying assumptions of a direction-independent, scaleinvariant, F STT, sq NL as well as a scale-invariant P ' . Note that, in the large F sq NL limit, due

Preheating and phase transitions

In standard preheating scenarios, a daughter or 'preheat' field ‰ is coupled to an inflaton " via some interaction term involving the two fields. If the inflaton potential exhibits a monomial shape at the stages following inflation, the inflaton oscillates around the minimum of its potential after inflation, inducing a non-adiabatic time evolution in the interactive mass of the preheat field. This leads to an e cient resonant production of the daughter species [121][122][123][124][125], the e ciency of which depends on the inflaton-daughter coupling, as well as on the details of the inflaton potential (see e.g. [126,127] for more recent analysis). This particle production mechanism is known as parametric resonance, and it corresponds to a non-perturbative, non-linear, and out-of-equilibrium e ect [128]. We speak about broad resonance when the choice of interaction and inflationary model leads to an excitation of the ‰ field modes within broad band(s) of momenta. In this case, a significant production of gravitational waves (GWs) takes place [22,[129][130][131][132].

In large field inflationary models, the daughter field is typically 'heavy' during inflation, as the inflaton field takes super-Planckian amplitudes. It is possible however, to find some coupling values for which the daughter field is light during most of the inflationary era, but becomes heavy only towards the last e-foldings of inflation (when the inflaton rolls down its potential towards smaller values). In this case, after inflation ends, ‰ displays amplified perturbations on super-horizon scales, just as the inflaton. At the onset of preheating, subhorizon vacuum fluctuations serve as an initial condition for parametric resonance, but these are super-imposed over almost homogeneous values ‰ i of the daughter field. 3 This is precisely the crucial ingredient for the development of anisotropies in the GW background. The value of ‰ i changes at super-horizon scales according to a variance

N , where N is the number of e-folds for which ‰ is a light degree of freedom, and H inf is the inflationary Hubble scale. Initial quantum fluctuations of the daughter field ‰ at sub-horizon scales are exponentially stimulated via parametric resonance. When non linearities become relevant in JCAP11(2022)009 the system, i.e. when ‰ back-reacts on the inflaton ", the dynamics of the sub-horizon modes ‰ k are influenced by the value of ‰ i within each given patch. The spatial distribution of the field ‰, and hence of the source of the GWs, will be then di erent at causally disconnected regions. As a result, a di erent amount of GWs is produced at each super-horizon region, in correspondence with the di erent values of ‰ i .

The anisotropies in the GW energy density spectrum from preheating have been studied in detail in the scale invariant model V (") = 1 4 ⁄"foot_1 + 1 2 " 2 ‰ 2 [20,21], chosen because of its computational convenience. GW anisotropies should be however a relatively common phenomenon arising in other preheating scenarios, as long as the appropriate conditions are met. In the mentioned scenario, the lightness of ‰ before the last e-folds of inflation is guaranteed if the coupling constant is taken to be g 2 /⁄ ≥ O(1). The dynamics of preheating proceeds as usual, but the initial conditions at the onset of parametric resonance are such that at each super-horizon volume there are di erent values ‰ i , drawn from a Gaussian distribution with variance

N . In practice one just needs to run simulations with free values of ‰ i , simply restricted to

Employing the 'separate Universe' approach, refs. [20,21] compared the peaks of the GW energy density spectrum from simulations with di erent initial values of ‰ i , run for the choice g 2 /⁄ = 2. While the GW backgrounds were always peaked at the same frequency, as expected, the peak amplitudes of the GW spectra di ered significantly. For example, in the left panel of figure 1 we show two GW spectra obtained for slightly di erent values of ‰ i , and it is clearly appreciated that one amplitude is larger than the other by a factor ≥ 2-3. In other words, the actual value of ‰ i influences the evolution of the sub-horizon gradients of ‰, and hence the production of GWs. To be concrete, GW was observed to vary up to a factor ≥ 5 between slightly di erent values of ‰ i (the non-linear dynamics is actually chaotic [133], so small variations of ‰ i can lead to a large variation of sub-horizon dynamics of the modes ‰ k ). The level of anisotropy produced in the energy density of the resulting GW background is characterized by the angular power spectrum C GW ¸of the relative GW spectral energy-density fluctuation [cf. eq. (2.9)], which can be written as a function of the ‰ i values. A general formula applicable to all scenarios characterized by a light spectator field during inflation is [20,21] 4

where "‰ i © ‰ i ≠ ‰ i , with ‰ i the mean value over the currently observable universe. This implies that the angular power spectrum of the GW energy density anisotropy is scale invariant, i.e. characterised by a plateau at small multi-poles, ¸(¸+ 1) C GW ¸Ã const., analogous to the large angular scale Sachs-Wolfe plateau for the temperature anisotropies in the CMB. In the analysed preheating scenario, the relative amplitude of the GW energy density spectrum, for a reference value of ‰ i = 3.42 • 10 ≠7 M Pl (here M Pl ƒ 1.22 ◊ 10 19 GeV is the Planck mass), was found to have spatial fluctuations as Ò l (l + 1) C GW ¸= 0.017 ± 0.003. For other values of ‰ i the anisotropy amplitude is also similar, always at the O(1) % level, see right panel of figure 1. For comparison, recall that the relative amplitude of the CMB temperature fluctuations is of the order of O(0.001) %. The GW anisotropies obtained in this model are therefore very large.

JCAP11(2022)009

3.2×10 -7

3.4×10 -7

3.6×10 -7

3.8×10 The red dot shows the amplitude for the reference value ‰ i = 3.42 ◊ 10 ≠7 M Pl . Both plots are taken from ref. [21].

The details of the GW anisotropy, if ever observed, could provide a powerful way to di erentiate between di erent inflationary and preheating sectors. The GW background from scale invariant preheating just discussed, is however peaked at very large frequencies [22], way above the observational frequency window accessible to LISA (or to any other ground-or space-based planned detector for this matter). So the example mentioned only serves as a proof of principle, at least for what can be detected in the foreseeable future. On top of this, it is important to remember that a quartic potential model is also ruled out by CMB data

The mechanism just described corresponds to the imprint of intrinsic anisotropies in the energy density of the GW background from preheating. However, in general, other e ects causing anisotropy can be also present. As a matter of fact, any causally sourced GW background will exhibit, in general, anisotropies in the spatial distribution of its energy density at large scales. This is simply due to Doppler, Sachs-Wolfe, and Integrated Sachs-Wolfe e ects [15], similarly as the anisotropies arising in the photons of the CMB. This type of anisotropies concern actually not only preheating, but also GW backgrounds from cosmological phase transitions, and in general from any causally driven mechanism creating GWs after inflation (as well as inflationary GWs themselves). For cosmological adiabatic perturbations, the fluctuations in any causally produced GW background will simply follow those in the CMB, and hence they are expected to be very small [23,24], of the order of ≥ 10 ≠5 . If primordial fluctuations carry however an isocurvature component, this need no longer be true. Ref. [24] has recently shown that in non-minimal inflationary and reheating settings leading to large non-Gaussian perturbations, a non-Gaussian GW background is also expected, even when the rest of the cosmological fluids inherit predominantly Gaussian fluctuations. Primordial isocurvature perturbations can survive in the GW background say from a cosmological phase transition, exhibiting significant non-Gaussianity, while obeying observational bounds from the CMB or Large-Scale Structure surveys. Probing such inherited non-Gaussianity in causally generated GW backgrounds seems to be however a marginal possibility at LISA [24], and rather more futuristic proposals such as the detectors DECIGO or BBO are needed.

Topological defects

Gravitational wave sources with an inhomogeneous spatial distribution would lead to anisotropies in the SGWB, in addition to the anisotropies induced by the nature of spacetime along the line of propagation of the GWs. An inhomogeneous distribution of cosmic strings, formed generically [134] in the early Universe as a result of a phase transition, followed by a spontaneous symmetry breaking characterized by a vacuum manifold with non-contractible closed curves, will lead to anisotropies in the SGWB.

Several studies [25][26][27] in the literature have calculated the anisotropies induced by a network of Nambu-Goto cosmic string loops, addressing the question of whether the model for the loop distribution will a ect the angular power spectrum. It has been shown that while the amplitude of the resulting power spectrum C ¸depends on the model of the loop network, the anisotropies are driven by local Poisson fluctuations in the number of loops, and the resulting angular power spectrum is spectrally white (i.e., C ¸= constant with respect to ¸), regardless of the particular loop distribution [25].

We show in figure 2 the amplitude of the SGWB angular power spectrum as a function of the string tension Gµ for three cosmic string loop distributions, dubbed "Model 1, 2, 3". The first, Model 1, is the original one-scale model where all loops have the same size set by a free parameter -, chosen here to be -= 10 ≠12 . While this model is rather obsolete, we illustrate it here since it has been shown that it leads to significant anisotropies in the PTA frequency band [26]. Models 2 [135,136] and 3 [137,138] are based on di erent computer simulations and they di er on the way they model the production and cascade of loops from the super-horizon cosmic string network.

We find that, regardless of the adopted cosmic string loop model and the considered string tension, the predicted angular power spectrum C ¸is too small to be detected with LISA. Note that in both models 2 and 3, the monopole should be detectable by LISA for Gµ & 10 ≠17 [99] (though the presence of astrophysical foregrounds reduces the sensitivity somewhat to Gµ & 10 ≠16 [139]).

Aside from the extra-galactic population of cosmic string loops discussed above, several authors have studied the possibility of a population of loops being captured in the halo of the Milky Way [140][141][142]. These loops would then give rise to an anisotropic SGWB signal which would trace the density profile of the galactic halo. Using the results from ref. [142] we calculate here the corresponding C ¸spectrum, which is shown in figure 3. Again, this signal is too weak to be detected by LISA.

JCAP11(2022)009

the fractional GW energy density

in terms of the critical energy density of a spatially flat universe fl c = 3H 2 M 2 p and the curvature perturbation power spectrum P ' .

The predicted amount of angular anisotropies can be determined by the two-point correlation function of the density field fl GW in di erent angular directions. For a Gaussian curvature perturbation one expects those to be undetectable, given the capability of the LISA experiment to measure anisotropies between spatial points separated by non-negligible fractions of the present horizon [101]. Indeed, according to the Equivalence Principle, the anisotropies will be highly suppressed by a factor (k PBH |x ≠ y|) ≠2 π 1, since the characteristic scales of the scalar perturbations are much smaller than those spatial distances, k PBH |x ≠ y| ∫ 1, and the emission takes place near horizon crossing.

This conclusion does not hold in the presence of primordial non-Gaussianity correlating short (' PBH ) and long scales (' L ). Indeed, large scale modulation of the power spectrum may lead to anisotropies at large-scales imprinted at formation [101]. Assuming a local, scaleinvariant, shape of non-Gaussianity ' = ' g + 3 5 f NL ' 2 g and keeping into account propagation e ects (see section 2.3 for details), one can compute the two-point correlation function of the GW energy density contrast in spherical harmonics " GW,¸m as in eq. (2.9), obtaining

in terms of the power spectrum at large scales P ' L and the momentum dependent non-Gaussian parameter

Non-Gaussianity in the curvature perturbation is constrained to fall in the range ≠11.1 AE f NL AE 9.3 at 95% C.L. [151] by the Planck collaboration. It is important to stress that its presence would also generate a significant variation of the PBH abundance on large scales given the impact of long modes inducing a modulation of the power on small scales. As isocurvature modes in the DM density fluid are strongly constrained by CMB observations [152], one can put an upper bound on the fraction of the Dark Matter (DM) in our universe composed by PBHs formed in the presence of non-Gaussianities [153,154].

For a monochromatic and lognormal power spectra of curvature perturbations at small scales, peaked at the LISA maximum sensitivity frequency, the GW anisotropy are plotted in figure 4, where the coloured region identifies the range of parameters allowed by the Planck constraints and the dot-dashed lines identify the present epoch GWs abundance evaluated at the peak frequency. The non-linear parameter has been assumed to be f NL > ≠1/3 to avoid the inconsistency of the perturbative approach in the PBH abundance computation happening at larger negative values, see the related discussion in [155,156]. The main finding JCAP11(2022)009 is that a large fraction of DM composed by PBHs would imply a highly isotropic and Gaussian SGWB, up to propagation e ects. On the other hand, the detection of a sizeable amount of anisotropy in the signal related to the PBH scenario would imply that PBHs can account only for a small fraction of the DM in the universe [101].

Propagation e ects

Independently from the initial anisotropies in the SGWB of cosmological origin, we do expect a minimal level of anisotropies in all of the scenarios described above due to the propagation of GWs through (large-scale) cosmic inhomogeneities, while travelling from the generation surface till the observation point. Such anisotropies represent an unavoidable, irreducible contribution which indeed carries precious cosmological information, being sensitive to the evolution of cosmological perturbations and to the initial conditions from which cosmic structures formed. Employing the general formalism of Boltzmann equations explained above, the underlying cosmological perturbations leave specific imprints in the statistics of the SGWB anisotropies in terms of, e.g., angular power spectra.

From eqs. (2.15), we can infer some properties about the SGWB anisotropies due to their propagation through cosmological perturbations: similarly to CMB, gravitons are a ected by the Sachs-Wolfe contribution, which represents the energy lost by a graviton which escapes from a potential well, and by the Integrated Sachs-Wolfe (ISW) e ect, due to tensor and scalar perturbations, the latter producing an anisotropy which is roughly proportional to the total variation of the potentials + . An important point to stress here is the "initial" time ÷ i , which has an impact both on the SW and on the ISW contributions. The numerical evaluation of the angular power spectrum for the cosmological SGWB has been performed in [109] (see also [120]), modifying the publicly available code CLASS, usually employed for the computation of CMB anisotropies [157] and adapting it to the SGWB.

In figure 5 we report the angular power spectrum of the cosmological SGWB due to propagation e ects sourced by scalar perturbations and we compared it to the CMB one coming from temperature anisotropies. We can see that the SGWB spectrum shows a larger amplitude compared to the CMB. This can be explained considering the graviton "decoupling" time, which occurs earlier compared to CMB photons and so gravitons feel for longer time JCAP11(2022)009 the propagation e ects. In such a figure, we also report the contribution from the SW and the ISW separately, to show their behaviour at di erent angular scales.

From the left plot we can see that at large angular scales (i.e., low ¸), the SW contribution is dominating while moving to smaller scales (i.e., ¸& 100), the ISW contribution starts to be larger. On the other hand, from the right plot we can quantify the expected di erence among the CMB and SGWB anisotropies.

Interestingly enough, by measuring or constraining angular anisotropies of the SGWB, it is also possible to probe the level of primordial non-Gaussianity possibly present both in the scalar and tensor cosmological perturbations through which the SGWB propagates. Indeed such primordial non-Gaussianity will left be imprinted into the GWs passing through the background large-scale underlying inhomogeneities, similarly to what happens for CMB photons. This entails to go beyond the power spectra statistics and to compute higher-order correlation functions, such as the angular bispectra of the GW energy density fluctuations [16,17]

where b GW ¸¸Õ ¸ÕÕ is the so-called reduced bispectrum (see e.g. [158,159]). For example, as shown in [16,17], in the case of primordial local non-Gaussianity in the curvature perturbations

' g (x) being the linear Gaussian part of the perturbation, one finds

It is important to stress that similar results follow in the case of primordial non-Gaussianity in the large-scale tensor perturbations. Therefore the 3-point correlation function of GW energy density anisotropies provides for the first time a way to probe at interferometers primordial non-Gaussianity of large-scale tensor modes through the imprints that the latter leave in the spatial distribution of GW energy density as described by the second equation JCAP11(2022)009 in (2.15). Indeed, for a su ciently high SGWB, it might happen that primordial (scalar/tensor) non-Gaussianity can be measurable through the detection of the SGWB anisotropies at interferometers, even in cases where such primordial non-Gaussianity are not measurable at CMB scales.

As it is clear from the previous discussion, the seeds that give rise to anisotropies during the propagation, are the same for photons and gravitons. This naturally induces a cross-correlation among these two messengers. Recently, the cross-correlation between CMB and SGWB anisotropies induced during the propagation has been studied in [110,120], and focusing on the initial anisotropy in [102,103]. All these studies have shown that such a cross correlation signal will be within the reach of the LISA detector, and it will be extremely useful to extract information on cosmological parameters, pre-recombination physics and the non-linear parameter f NL to measure primordial non-Gaussianity.

Astrophysical sources of anisotropies

The astrophysical stochastic gravitational-wave background (AGWB) is generated by the incoherent superposition of signals emitted by a large number of resolved and unresolved astrophysical sources. In di erent frequency bands, several astrophysical sources can contribute to the AGWB, as merging stellar-mass black hole (SOBHB) or binary neutron stars (BNS) [160][161][162][163][164][165][166][167], super-massive black hole binaries (SMBHB) [168], rotating neutron stars [169][170][171], stellar core collapse [172,173] and population III binaries [174] (see, e.g., [1] for a review). As the cosmological GW background, also the AGWB is characterized by an isotropic energy density contribution and through the spatial angular power spectrum encoding its anisotropy.

Based on the recent observations of merging black holes and neutron star binaries by the LIGO and Virgo detectors [175][176][177][178][179][180], it is estimated that the stochastic background from unresolved stellar-mass compact binaries may be detected within a few years of operation of such a network [181]. Its anisotropic component is constrained by LIGO/Virgo observations up to ¸= 4 [182] which results in upper limits on the amplitude of the dimensionless energy density per units of logarithmic frequency in the range GW (f = 25Hz, ) < 0.64-2.47 ◊ 10 ≠8 sr ≠1 for a population of merging binary compact objects, where denotes the angular dependence. Methods to measure and map the AGWB in the LIGO and LISA frequency ranges are discussed in [5,10,11,53,54,[183][184][185][186][187].

Traditionally, the energy density of the AGWB has been modeled and parameterized under the assumption that both our universe and the distribution of sources are homogeneous and isotropic (see e.g. refs. [1,164]). This is a rather crude approximation: GW sources are located in galaxies embedded in the cosmic web; moreover, once a GW signal is emitted, it is deflected by the presence of massive structures, such as galaxies and compact objects. It follows that the energy flux from all astrophysical sources has a stochastic, directional dependence.

The first prediction of the AGWB angular power spectrum was presented in [33,34] following the framework introduced in refs. [30,31]. This framework is flexible and splits the cosmological large-scale structure and sub-galactic scales so that it can be applied to any source contributions and any frequency band. The astrophysical dependence of the angular power spectrum on the detail of the underlying astrophysical model has been studied in [25,[35][36][37][188][189][190][191] and di erent formal aspects of the derivation of anisotropies and their interpretation are discussed in [15,18,19,32,187]. The relative importance of cosmological and astrophysical e ects depends on the frequency band chosen, hence o ering the possibility JCAP11(2022)009 to distinguish di erent astrophysical processes. Due to their stochastic nature, anisotropies can be statistically characterized in terms of their angular power spectrum and they also correlate with other cosmological observables such as weak lensing, galaxy number counts and CMB anisotropies.

The study of the cross correlations with electromagnetic observables provides complementary information and might improve the signal to noise of the anisotropic searches [36,187,192,193]. Moreover, cross-correlating the background that collects contribution from sources at all redshifts along the line of sight, with EM observables (such as galaxy number counts) at a given redshift, allows one to get a tomographic reconstruction of the redshift distribution of sources [33,36,187,192].

GW energy density for astrophysical sources

The total present-day GW energy density per logarithmic frequency f o (where f o = q, see the previous section) and solid angle o along the line-of-sight ŝ (note that ŝ = ≠n) of the AGWB is defined as [149,194].

and it represents the fractional contribution of GWs to the critical energy density of the Universe today fl c,0 = 3H 2 0 /(8fiG); dfl tot GW is the total energy density of GWs in the frequency interval of today {f o , f o + df o }. See also the definition in eq. (2.6). Such a quantity contains both a background (monopole) contribution in the observed frame, which is homogeneous and isotropic, i.e. Êtot GW = ¯ tot GW , and a direction-dependent contribution

As usual, we consider the local wave zone approximation at the source position: in other words, the observer "at the emitted position" is defined in a region with a comoving distance to the source "su ciently large" such that the gravitational field is "weak enough" but still "local", i.e., its wavelength is small w.r.t. the comoving distance from the observer ‰ (see also [32,195]). Considering an observer that measures a GW signal in a fixed direction n, one expects that the total gravitational energy density in such a direction is given by the sum of all the (unresolved) astrophysical contributions along the line of sight contained in a given volume dV e (ŝ) and can be expressed as

from which we can build the total gravitational energy density tot GW =

where

and [i] is related to the summation over all unresolved astrophysical sources that produce the SGWB. Here dA o is the unit surface element at observer [31]. The vector ◊ = {M h , M ú , m, ◊ú } represents all the parameters which are: the halo mass M h , the mass of stars that give origin to the sources M ú ; m indicates the masses of the compact objects and ◊ ú includes the astrophysical parameters related to the model (like spin, orbital parameters, star formation rate). In eq. (3.3), n

h is the (physical) number of halos at a given mass M h , within the JCAP11(2022)009

physical volume dV e , weighted with the parameters ◊ of the sources at x µ e . The letter "e" stands for "evaluated at the emission (source)" while "o" for "evaluated at the observer". Using the energy conservation we have

where we have redefined

GWo . Here T is the proper time of the observer and D A is the angular diameter distance. Now, defining the total GW density as

we can easily obtain the expression for the energy density

Here, x µ (‰) are the comoving coordinates in the real frame (the "physical frame"), where ‰ is the comoving distance from the source to the detector. The previous expression depends on the position, but we can also define a position independent, isotropic version of it by integrating over a spatial volume: we denote the corresponding quantity df tot GW /df o d o with a bar, as in section 2.1.

We use the observer frame where we perform observations (also called "cosmic GW laboratory" in [195]). This is the correct frame where we want to reconstruct 3D maps/catalogs of galaxies by using both EM and GW signals. Let us point out that if we use unperturbed coordinates, instead of the observer coordinates, we are not able to interpret correctly the correlation between for istance the AGWB and EM sources from observed galaxies. This could induce a wrong estimate of our results [32,195]. In particular, we consider coordinates that are flattened in our past gravitational wave cone so that the GW geodesic from the source can be defined with the following conformal space-time coordinates

Here, ÷ 0 is the conformal time today, ‰(z) is the comoving distance to the observed redshift and ŝ is the observed direction of the GW, i.e.

Using ‰ as an a ne parameter in the observer's frame, the total derivative along the past GW-cone is

Setting up a mapping between the observed frame and the "physical frame" in the following way x µ (‰) = xµ ( ‰) + x µ ( ‰), where x µ ( ‰) is a suitable linear perturbation that shifts the comoving four-coordinates from the real-space frame to the observed frame. Then using the decomposition of eq. (3.3), we obtain

⁄ N [i] (z, f e , ◊)

JCAP11(2022)009

with 3 the total comoving number density at a given redshift or ‰. Notice that, by construction, the quantity ¯ tot GW is isotropic. At linear order, we obtain the following AGWB energy density fluctuation

⁄ N [i] (z, f e , ◊ )

where a(‰)

and "f is the linear perturbation of the frequency of the GW due to the anisotropies. Finally, let us mention that when the integration along the line of sight is performed, one should also consider the normalized selection window function w(z), whose form depends, besides redshift, on the sensitivity/characteristics of the GW detector (see [32] and [38] for more details about the window function). So we finally have

and

Connection with the halo and stellar mass function and with the star formation rate. In general, the isotropic component of equation (3.10) is given by [196]

where dE

GWe /df e /d e is the energy spectrum per unit solid angle, p [i] ( ◊ ) is the probability distribution of the source parameters ◊ and R [i] is the observed comoving merger rate density of i-th unresolved type of source. In particular, the event rate (per unit of redshift) can be derived from the cosmic star formation rate. For instance, assuming for simplicity that the gravitational emission occurs shortly after the birth of the progenitor, it turns out that

where the (1 + z) term corrects the cosmic star formation rate (SFR) by the time dilation due to the cosmic expansion and dfl

ú /dT e is the (density) cosmic SFR in M § , Mpc 3 and yr ≠1 .

JCAP11(2022)009

Here ⁄ [i] (z, ◊) is a generic function which depends on the initial mass function M ú and, in general, on other parameters of the sources, as the halo mass M h . So then we have

GWe (z, f e , x µ e , ◊) df e dT e d e . (3.17)

Now, let us consider events with short emission (i.e. burst sources), as merging binary sources (BH-BH, NS-NS and/or NS-BH). Then we have

where dN

GWe /dT e dM ú is the merging rate of events for each halo and at a given stellar mass M ú , and the comoving density 3 of the halos can be rewritten as

i.e. the comoving density at a given M h . In order to give a very simple example, let us assume that nh and Nh are equal for all sources. In this case they do not depend only on M h and we can remove the index [i] and N h (M h , z) can be related to the fraction of mass F (M h , z) that is bound at the epoch z in halos of mass smaller than M h , i.e.

d Nh (M h , z)

where f(z) is the comoving background density (e.g., Press & Schechter (1974) [197], Sheth & Tormen (1999) [198] or Tinker (2008) [199] mass fraction). Following [200,201], it is possible to express the mass function in terms of the multiplicity function of halos g(M )

Physically, this quantity gives the fraction of mass that is bound in halos per unit logarithmic interval in mass. Finally let introduce the (mean) SFR that it is connected with

GWe in the following way

Note that s(M h , z) defined in [200,201] can be related to the SFR in the following way

JCAP11(2022)009

where we have defined

where P m (k) is the matter power spectrum today and S ¸are the source functions which include all the e ects described in eq. (3.24). The index [i] refers to the specific unresolved astrophysical source while the greek index stands for the various contributions to the GW energy density anisotropies.

To have some physical insight into the information encoded in the anisotropies of the AGWB, we can use the Limber approximation. The general expression of the angular power spectrum reduces to [30]

where ¸is the multipole in the spherical harmonic expansion, P (k) is the galaxy power spectrum, r is the (comoving) distance, related to the momentum k via the Limber constraint kr = ¸+ 1/2. Each astrophysical model predicts a functional dependence of the astrophysical kernel ˆr ¯ , defined as

It follows that the angular power spectrum depends on the astrophysical model chosen to describe sub-galactic physics. In particular, low ¸are sensitive to the low-redshift value of the kernel ˆr ¯ . The angular power spectrum of the anisotropies in the AGWB from merging stellar-mass binary BHs in the mHz band where LISA operates, has been computed in [37] using the astrophysical framework of [164]. The result for three di erent frequencies is shown in figure 6. It has been shown that AGWB anisotropies are very sensitive to sub-galactic astrophysical modeling. In particular, di erent descriptions of stellar evolution and black hole binary formation lead to fractional di erences in the angular power spectrum of anisotropies up to ≥ 50%, independently on the global normalization (monopole) [36,37].

Monopole and anisotropies contain complementary astrophysical information and studying the latter will allow one to break degeneracies between di erent astrophysical ingredients and potentially to constrain them separately.

Systematic e ects on the angular power spectrum

In estimating the anisotropies of the astrophysical gravitational-wave background, the finite number of the sources that contribute to the background at any given time and the very short time each of them spends in the frequency of the interferometer, induce a white noise (¸-independent) term, W, in the angular power spectrum C GW ¸:

where V is some direction-independent (due to statistical isotropy) function describing the variance, r H © 1/H 0 is the Hubble radius and C LSS ¸stands for the angular power of the . Angular power spectrum of anisotropies for three frequencies in the mHz band, for the reference astrophysical model of [37] and [36]. Multiplication by (¸+ 1) emphasizes the large-scale behaviour of eq. (3.30), while we multiplied the spectrum by the monopole amplitude to show the frequency scaling of anisotropies. The shaded region corresponds to the cosmic variance limit. Adapted from [37].

intrinsic, astrophysical anisotropy. The shot noise dominates over the true astrophysical power spectrum; the latter may be recovered with long enough observing runs and su cient removal of a large number of foreground sources [188].

To calculate the true, astrophysical angular power spectrum of a statistically-isotropic gravitational-wave background, a novel method, based on combining statistically-independent data segments, was proposed in [189]. The proposed estimator, constructed from the crosscorrelations between statistically-independent time intervals, reads (3.33) where N • © T /• (with T the total observing time) denotes the number of segments and

(3.34)

The estimator (3.33) is unbiased since

where the subscripts S, stand for performing first the cosmological and the shot noise average. In the limit of many data segments, N • ∫ 1, the estimator (3.33) has the lowest-variance

and, in this sense, it is the most e cient one. It is worth noting that the term W T in the above equation is the same as the one appearing in the mean of the standard estimator ÈC (std)

) is still a ected by the shot noise. However, now the shot noise only adds to the variance of the estimator and it does not a ect the angular power spectrum, as in the standard case.

¸-dependent response functions in the A, E, T channels

As shown by eq. (4.16), the anisotropic LISA response functions transform as spherical harmonics under rotations. One can therefore consider the ¸≠dependent response function

that is invariant under rotations, and therefore constant in time (it does not depend on the orientation of the LISA triangle). As we show in the next subsection, it provides an estimate for the response of LISA to a statistically isotropic SGWB, see subsection 4.4. From the properties in eq. ( 4. 19) we learn that

It is customary to consider linear combinations of the F i measurements considered so far

which we write more compactly as

These combinations (that we have normalized as in ref. [9], so that the rotation matrix associated with these transformations is orthogonal) diagonalize the noise covariance matrix, in the hypothesis that LISA is an equilateral triangle, with identical instruments at the vertices. In terms of the A, E, T channels the response function formally reads

We evaluate these linear combinations, accounting for the identities in eq. (4.19), namely

The resulting expressions acquire di erent forms for even and odd multipoles. Specifically, for odd ¸we find (4.13). In each term, we have kept up to the leading f -dependent term. We recall that the response functions are symmetric in the channels, that RȨ E = RA A , and that RȨ T = RA T .

and for even ŖA

We also note that the property in eq. (4.18) implies that the response fuction is symmetric in the channels, RO Õ O = RO O Õ . These expressions can be evaluated numerically, for arbitary frequency, or evaluated analytically in the small frequency regime. For the first few multipoles, we obtain the values in the table 1. In figures 7 and 8 we show instead a comparison between the full shape of the response functions and the small frequency expressions for these first multipoles.

Sensitivity to ¸-multipoles

Eq. (4.38) provides the SNR for the detection of a SGWB which is the sum of all possible multipoles contributions. Although we have not explicitly written it, the response functions R ¸m OO Õ (f ) also depend on time, as they are functions of the positions of the satellites. A full analysis of the separate contributions of the various multipoles would then require a component JCAP11(2022)009 separation, which is in practice the inversion of the time-dependent streams measured by the satellite to the multipole amplitudes p ¸m. We leave this discussion to section 6. Here we estimate the relative sensitivity of LISA to di erent ¸≠multipoles by assuming that only one multipole dominates the SGWB and that multipoles with the same ¸but di erent m are obtained from the same Gaussian statistics. This amounts to assuming a statistically isotropic SGWB, with correlators given by eq. (2.9).

Taking this into account, the expected SNR (4.38) can be written as a sum over the various multipoles, (4.39) where, for each multipole,

where we recall that the response function R OO Õ (f ) is the quantity defined in eq. (4.28) and rescaled as in eq. (4.12). In the following, we can work directly in terms of RO O Õ (f ) by rescaling the noise functions accordingly, see eqs. (B.14) and (B.15). Moreover, ad discussed above, the response function RO O Õ to a statistically isotropic signal is time-independent, so that the integral over time in eq. (4.40) simply results in the usual property that the SNR grows with the square root of the observation time. Finally, we factor out the uncertainty in the Hubble rate by dividing it by its rescaled value h and by considering the h 2 combination, as it is standard. This leads to

From this expression we define the "channel-channel" sensitivity

as well as the optimally weighted sum over the three channels

The total sensitivity to the ¸≠multiple is shown in figure 9 for multipoles up to ¸= 10. From this quantity, we can immediately obtain

(4.44)

We note that the curves shown in figure 9 are rescaled by Y 00 = 1/ Ô 4fi, in such a way that the curve shown for ¸= 0 coincides with the SciRD (Science Requirement Document) sensitivity curve for a homogeneous signal [203] obtained from summing over the A, E, T channels.

Sensitivity to kinematic anisotropies

Doppler anisotropies induced by the motion of the detector with respect to the SGWB rest frame count among the guaranteed features of the SGWB. In fact, already the early work [10], which sets the basis for the analysis of SGWB anisotropies with ground-based GW interferometers, estimated the prospects for ground-based detectors to measure the kinematic dipole of the SGWB. In this subsection we briefly consider the same question in the context of LISA.

The size and properties of kinematic anisotropies depend on the frequency profile of the rest-frame SGWB energy density GW (f ). This fact can be important for enhancing the amplitude of kinematic anistropies in certain early-universe scenarios where the SGWB has rich features, as the ones discussed in section 2.

We consider two cosmological frames: the first, denoted with S Õ , is comoving with the SGWB rest frame; the second, denoted with S, moves with constant velocity with respect to the rest frame S Õ . We assume that the SGWB density parameter in the rest frame, Õ GW (f ), is perfectly isotropic and depends only on frequency f . A boost transformation relates the SGWB density parameter in the rest frame S Õ to the one in the moving one S. We indicate JCAP11(2022)009 with v = -v (where -= v in units with c = 1) the velocity of the frame S with respect to the rest frame S Õ .

In the technical appendix C we derive the resulting expression of an anisotropic SGWB energy density GW (f, n) as a function of the rest-frame density Õ GW (f ). Assuming that the parameter -is small, we can Taylor expand up to second order in -and write

The functions of frequency M , Q, D, control respectively the contributions of kinematic e ects to the monopole, dipole, and quadrupole of GW energy density in the detector frame. They read

In analogy with CMB literature, we introduce the SGWB spectral tilts

)

The expressions (4.46), (4.47), (4.48) quantitatively demonstrate that enhanced spectral tilts can amplify kinematic anisotropies in certain scenarios. We plot in figure 10 the SNR for LISA to detect the kinematic dipole and quadrupole induced by a scale-invariant profile of Õ GW (f ) = constant in the SGWB rest frame. Notice the di erent vertical scale in the two plots, due to the fact that LISA sensitivity to the quadrupole is a factor ≥ 10 3 -better than that to the dipole, as discussed in the previous sections.

JCAP11(2022)009

For the sake of generality, we consider a power-law SGWB spectrum peaking at a fiducial multiple L only, parameterized by the logarithmic amplitude log 10 A c at a pivot frequency f c = 2.5 • 10 ≠3 Hz, that is chosen close to the frequency where LISA has the best sensitivity, and by a spectral index -,

(5.2)

For each multipole ¸and channel combination OO Õ , we assume a Gaussian likelihood over the averaged data L ¸given by

where N c is the number of data segments in the analysis; the sum runs over frequencies (or frequency bins) f k , D OO Õ ,¸d enotes the averaged signal over the data segments in the channel combination OO Õ , and D th OO Õ ,¸i s the theoretical ansatz for the data,

where ROO Õ ,¸i s the frequency response of the detector and Ñ OO Õ is the noise as defined in the previous section expressed in Omega units. The variance can be expressed in terms of the theoretical ansatz as

In practice, instead of summing over channels in the likelihood, we consider a single data vector and compare it with the e ective noise combination defined by eq. (4.43) and shown in figure 9, and drop the OO Õ channel indices in what follows.

Assuming a fixed noise model, the Fisher information matrix for the likelihood defined in eq. ( 5.3) is simply

where ◊, fl are a combination of the signal model parameters log 10 A c and -, and the corresponding partial derivatives are

(5.6)

The estimated LISA sensitivity to a single-monopole power-law SGWB defined in eq. ( 5.2) has been represented on figure 12 for the monopole (¸= 0), dipole (¸= 1) and quadrupole (¸= 2), for a series of fiducial values of the SGWB amplitude log 10 A c and spectral index -. In all cases, the standard deviation for each parameter is considered marginalised over the other one (i.e. taken from the diagonal elements of the covariance matrix C ◊fl , defined as the inverse of the Fisher information matrix).

As one can see in figure 12, for ¸= 0, 2, su ciently high log-amplitudes are recovered independently of the sign of the spectral index (but enhanced by stronger indices), due to the pivot frequency being chosen to approximately coincide with the peak in sensitivity at both multipoles. In contrast, for ¸= 1 positive spectral indices enhance the recovery of the amplitude. This is on the one hand because the corresponding sensitivity peaks at slightly JCAP11(2022)009 larger frequency with respect to f c ; and on the other hand because of the milder slope of the sensitivity with respect to ¸= 0, 2 towards high frequencies, so that the power law is closer to the high-frequency noise spectrum for lower |-| in ¸= 1 than in ¸= 0, 2 (see figure 9). For all multipoles, the spectral index is obviously recovered more e ectively for higher log-amplitudes. In the optimal case of a signal amplitude of order GW (f = f c )h 2 = 10 ≠9 , a null spectral index could be reconstructed with an uncertainty of order 10 ≠3 , 10 ≠2 , 10 ≠3 for the ¸= 0, 1, 2 multipoles respectively. In the more pessimistic case of GW (f = f c )h 2 = 10 ≠13 , for ¸= 0, 2 the spectral index could be reconstructed with an uncertainty of order 0.1 or greater for largely positive or negative values of it, but this uncertainty approaches order one for SGWB spectra with a spectral index between ≠1 and 1. In such a low-amplitude JCAP11(2022)009 scenario, LISA will thus be more sensitive to models with a strongly varying SGWB spectrum. Notice how, for the same reasons described in the previous paragraph, the dipole ¸= 1 is more sensitive towards positive spectral indices, whereas for ¸= 0, 2 the accuracy is almost symmetric with respect to the sign.

Map-making method

In this section we briefly describe the maximum likelihood map-making method for stochastic backgrounds proposed in [54] and provide estimates for the noise power spectrum N ¸obtained by simulating and mapping the noise directly in the sky domain. Recently, another method to map the gravitational-wave sky with LISA has been developed and it is based on a Bayesian algorithm to map the power of the SGWB using a spherical harmonic approach [204].

The maximum likelihood map-making with GW detectors relies on the specific scan strategy of the interferometer array, which describes how the sky signal is sampled as a function of time. The reconstruction of the GW sky and the angular resolution at which it may be achieved then depend on the amount of modes sampled throughout the whole duration of the observation.

To simplify the mapping procedure we assume that the anisotropic SGWB signal intensity I has a simple power-law spectral shape which may be factored out, such that

where E(f ) = (f /f 0 ) " , and f 0 is a specific reference frequency.

As for the scan strategy, we assume the spacecrafts follow three heliocentric, quasicircular orbits remaining at a constant arm-length distance from each other, and that the noise in the detector is well understood. Specifically, as in section 4, the noise is modelled by two contributions: acceleration noise and interferometer noise. For more details, see equations (30) and (31) in [54] and the description of the noise parameters in the o cial mock data release [205].

The data vector is defined (similarly to equation (4.32)) as d = Rh + n, where the first term specifies the pure signal component, made up of the contraction between the linear detector response R (see equation (A.4) in the appendix) and the SGWB strain h, and n is the noise component. We keep the formalism general here for the sake of conciseness; note that for multiple LISA TDI channels d is a vector in TDI space. To estimate the intensity I(f 0 , n) directly, the data are considered in frequency space: d(f ) with f belonging to the appropriate frequency interval observed by LISA. We assume the noise is zero-mean and Gaussian with covariance N = n ¢ n and the signal component is also Gaussian such that the total, signal plus noise, covariance of the data is C = A Ĩ + N. Here A is the operator that describes the response of the detector to the strain intensity. This can be integrated in time and projected onto pixel or spherical harmonic space. In the case of multiple TDI channels it represents the full correlated response matrix. Also note that in the case of multiple correlated TDI channels, e.g. X, Y , and Z as considered in [54], the noise covariance is the full correlated covariance matrix with the auto-correlated noise model for the diagonal and cross-correlated model for the o -diagonal terms. Ĩ is the observed realisation of GWB intensity on the sky. The likelihood L of the data is then

JCAP11(2022)009

as described in [206] for the case of CMB mapping, it is possible to find the iterative solution which maximises L,

where F is the Fisher information matrix and D © d † ¢ d. In practical applications, given the constraints on scan strategies and response functions for gravitational wave observations, the Fisher matrix will need to be regularised in order the correct iterative solution to be found using eq. ( 6.3). Here, the intensity of the sky-map is indexed generically by -such that I -are the set of "parameters", on which the signal component of C depends, that have to be estimated. For any particular application, the indices -could stand for either map domain pixels "p" or spherical harmonic domain multipoles "¸, m". The choice here is to work in the map domain, as it is not overly expensive in this case, and leads to a clearer understanding and regularisation of the Fisher matrix, as explained in [54,186]; in the latter, tests of the method and a regularisation technique are presented. In principle, the Ĩcould be estimated for a frequency band as narrow as the resolution permits, however this would result in a poorly regularised problem. To improve this the estimation must be done using wider frequency bands. Here we simply adopt the broad-band limit by assuming a spectral shape as in eq. ( 6.1). In practice this means the traces in eqs. (6.3) and (6.4) include a sum over the full frequency response, such that the final estimate is given with respect to a single reference frequency Ĩ-(f 0 ). Note that the term ˆC ˆI-Õ = A -represents the directional quadratic response of the detector, and is equal to A = R ¢ R up to appropriate normalisation factors. A -will, in general be time dependent, presenting a sky modulation with period of one year. Hence, to apply this mapping algorithm e ectively, the data must be segmented into short observation time-frames, throughout which the sky response is assumed to be constant and for which the noise can be estimated accurately. The duration of each segment • also sets the lower bound on the observable frequency window, hence there is a trade-o between frequency and (potential) sky resolution: in principle a shorter time window • will allow access to higher pixel resolution, but it will also fix the frequency resolution to 1/• . Assuming statistical independence between time-frames, the two traces in eqs. (6.3) and (6.4) are then obtained by averaging 6 over all the frames available.

For the purpose of this paper, we use the method described above directly in pixel space to calculate the noise power spectrum N ¸, assuming the LISA noise curves as in [54], over an observation time of one year. This is simply done by running the iterative map-maker over a noise-dominated data set, such that the Fisher matrix -setting now -© p -reduces to

The statistically isotropised, angular power spectrum of the noise N ¸can then be estimated by expanding and inverting the Fisher matrix [207]:

For actual data this average would be weighted by noise estimates but here we assume the noise is constant and uncorrelated between time-frames. GW units at reference frequency f 0 = 0.01 Hz. Note the N ¸here is in good agreement with the curve shown in [208], figure 4, taking into account that the observation time considered here is one year, whereas in [208] it is four.

JCAP11(2022)009

where the linear operator Y ¸m,p © Y ¸m( p) and similar for its adjoint. This is similar in spirit to what is presented in [208], however note that here the estimate is obtained by simulating and integrating the full scan strategy, without assuming the noise is isotropic to begin with. The N ¸estimates are shown in figure 13, in units of 2 GW at reference frequency f 0 = 0.01 Hz for a quicker comparison with signal models. This matches the convention chosen in [208], and the ¸= 0 mode of the two estimated noise power spectra matches well. Note that in figure 13 an observation time of four years is assumed.

The discrepancy at odd ¸s between the noise power spectrum obtained by the inverse Fisher matrix and the analytic computation presented in [208] is in part due to the substantially di erent sampling of the m modes. While these are marginalised over at a fixed time in equation (4.28) to obtain the instantaneous ¸-mode response, in the map-making procedure these are kept into account, and through the scan strategy contribute to breaking the degeneracies in the odd ¸-modes of (instantaneous) RO O . The higher-¸section of the curve is highly dominated by the conditioned inversion; in fact, the noise covariance matrix found in this case is highly singular, and ≥90% of its eigenvalues has been discarded to produce the curve in figure 13. The conditioning has a stronger impact on the higher-¸end of the angular spectrum, as the response of the detector is weaker at higher angular scales. It is therefore di cult to compare numerical estimates of the e ective noise at di erent angular scales to the result of analytical estimates. This is an active field of research, and it is clear that a robust regularisation scheme will be required when attempting to reconstruct the higher modes of the angular power spectrum with this configuration of the LISA instrument. Inevitably, the presence of complicating factors such as non-stationarity, noise uncertainties, and scanning systematics in real data will exacerbate this issue and mapping techniques will require significant developments in order to reconstruct anisotropies as optimally as possible.

JCAP11(2022)009 7 Conclusions

The anisotropies of the SGWB represent a powerful tool to characterize and distinguish the di erent sources of GWs. We have seen in this paper how di erent GW sources are characterized by di erent angular spectra. Such anisotropies have mainly two contributions: one directly related to the production mechanism of each particular GW source, and one being an e ect of the propagation of GWs on our perturbed Universe, which is common for all the GWs sources. We have made an overview of the main cosmological and astrophysical sources characterized by anisotropies, that are expected to be present in the LISA frequency band. We have presented the angular spectrum for di erent cosmological backgrounds (i.e., inflation, phase transition, PBH and cosmic strings) and an astrophysical one (Solar Mass Black Hole Binaries). We have then built a SNR estimator to quantify the sensitivity of LISA to di erent multipoles. To do this, we have computed the responses of LISA in harmonic space as functions of frequency for the AET TDI channels. We have also derived the analytic form of the responses in the low frequency limit. It is important to stress that, when anisotropic signals are considered, both the auto-correlation responses (i.e., AA, EE) and cross-correlation ones (i.e., AE, AT) are di erent from zero. We have shown how LISA will have a better sensitivity to detecting a quadrupole (i.e., ¸= 2) than it will for the dipole (i.e., ¸= 1). We have quantified the SGWB energy density required to observe the kinematic dipole and quadrupole induced by the motion of the LISA detector with respect to the SGWB rest frame. We found that an -GW ≥ 2 ◊ 10 ≠11 is required to observe a dipolar signal, while the sensitivity to the quadrupole is a factor ≥ 10 3 -larger than that to the dipole. We have also performed a forecast of the detectability of the lowest multipoles of the SGWB angular power spectrum through a Fisher matrix analysis. We have shown that for ¸= 0, 2 su ciently high amplitudes are recovered independently of the sign of the spectral index (but enhanced by stronger indices). On the other hand for ¸= 1, positive spectral indices enhance the recovery of the amplitude; conversely the spectral index is recovered more e ectively for higher log-amplitudes, for all multipoles. Finally, taking into account the LISA motion and the sky scan strategy, we have applied a maximum likelihood map-making technique to extract the noise angular power spectrum N ¸as a function of the multipole ¸. The LISA sensitivity and angular resolution will allow to detect the anisotropies of the SGWB, opening the possibility to use them in the process of characterization of the SGWB, and also to study their correlation to other cosmological tracers such as the Cosmic Microwave Background [102,110,112,120] and galaxies as tracers of the Large-Scale Structure [209][210][211][212][213]. This represents an exciting possibility to use LISA to explore our universe in a completely new perspective.

A Properties of the anisotropic response function

We insert the expression (4.1) into (4.6) and perform the line of sight integration, obtaining where we have defined hA

Lengthy but straightforward algebra then leads to the TDI combinations (defined in eqs. (4.9) and (4.10)):

In this expression, f ú is the frequency defined in eq. (4.13), and we have introduced the function

as well as the function W which is di erent for the two TDI combinations:

for TDI 1.0 e ≠2ik/fú ≠ 1 , for TDI 1.5 (A.6) The correlation between the TDI measurements in eq. (A.3) is expressed by eq. (4.12). As stated in the main text, the anisotropic LISA response function in eq. ( 4 To prove the first property, we consider a rigid rotation of the instrument, for which the position of the three satellites changes according to xi ae Rx i .

We perform an analogous rotation on the integration variable in eq. (4.15), and, accounting for the fact that scalar products of two vectors are invariant under a rotation we arrive to

The behavior of the polarization operators under a rotation can be found in eq. (A.17) of ref. [214]. Using that result, we can see by direct computation that, for any two unit vectors û, v,

As a consequence, the rotation matrix is eliminated from the last two lines of eq. (A.7), and one is left with the rotation of the spherical harmonic, from which eq. (4.16) is obtained. Inserting the expression in eq. (A.2) for M in eq. (A.5), we see that

. An identical property is shared by the GW polarization operators, and therefore by the functions G A . As a consequence,

We start from eq. (4.15) for R¸m ji . We send k ae ≠ k in the integrand, and we use the property that we have just proven. We arrive to an expression that is identical to the r.h.s. of eq. (4.15), with the only di erence that the argument of the spherical harmonic is = k. From the transformation of the spherical harmonics under parity we then obtain the property in eq. (4.18).

Let us now prove the property in eq. (4.21). We place the LISA satellites in the xy plane, with the center of LISA at the origin, and we simultaneously send the positions of the satellites xi ae ≠x i , and change sign to the integration variable k in eq. (4.15). These two operations do not change the scalar products k • l entering in the integrand of eq. (4.15). Therefore, they do not modify the first factor nor the second line of the integrand of eq. (4.15), JCAP11(2022)009 but only a ect the spherical harmonics. Next, we rotate the LISA triangle and the integration variable by 180 ¶ around the z-axis. These two operations only a ect the spherical harmonic in the integrand of (4.15). Under both sets of operations, the spherical harmonic changes to , as well as eq. (A.9). We end up with the conjugate of the r.h.s. of eq. (4.15) times the factor (≠1) ¸+m . From the last property that we have proven, we know that the response function is non vanishing only if ¸+ m is even, namely only if this additional factor is one. This proves the property in eq. (4.20).

JCAP11(2022)009

As we did for the expectation value, we can then substitute the functions " • with the Dirac delta-function, since the time • is much greater than the inverse frequencies. We then perform the integrals over f and f Õ , the sums over O ÕÕ and O ÕÕÕ , and we relabel k ae f and t av ae t in the resulting expression (B.10) Using the fact the noise is an even function of f , this expression can be finally written as eq. (4.37) of the main text.

Starting for the expressions eqs. (4.36) and (4.37), for, respectively, the expectation value and the variance of the estimator (4.34), it is convenient to define where, making use of eqs. (4.36) and (4.5),

We then see that the SNR is maximized by Q OO Õ (t, f ) = c ◊ " ú OO Õ (f, t), where c is an arbitrary constant that we can set to one. This leads to eq. (4.38) of the main text.

We conclude this appendix by providing the LISA noise functions used in our computations, referring the interested reader to ref. [9] for a detailed discussion of these quantities. We derive the expressions for the anisotropies of the SGWB induced by a boost transformation. We use the same methods as in refs. [215][216][217][218]. We consider two frames: the first, denoted with S Õ , is comoving with the SGWB rest frame; the second, denoted with S, moves with constant velocity v with respect to the rest frame S Õ . A boost transformation relates the SGWB density parameter in the rest frame S Õ to the one in the moving one S. We denote with f Õ the frequency of the GW in the SGWB rest frame. and with nÕ the unit vector denoting its direction. The frequency f in the frame in motion is associated with f Õ by a Lorentz transformation reading

where v = -v is the relative velocity of the two frames, and -= v in units with c = 1.

In order to compute how the GW energy density changes under a Doppler boost, we work in terms of the GW distribution function, denoted with Õ (f Õ ). We assume for simplicity it only depends on the frequency f Õ in the SGWB rest frame (i.e. the SGWB is perfectly isotropic in the frame S Õ ). We express the number of gravitons for unit of phase space in the rest-frame S Õ as:

where dV Õ corresponds to the infinitesimal volume containing gravitons with propagation vector nÕ in the element of measure df Õ d 2 nÕ . It is not di cult to prove that the combination f Õ2 df Õ d 2 nÕ dV Õ is invariant under boosts. In fact, we have the relations f Õ = D ≠1 f , d 2 nÕ = D 2 d 2 n, dV Õ = D dV (see refs. [217,218]). On the other hand, the number of gravitons (C.3) is independent of the frame, and dN Õ = dN . Hence [215] Õ (f Õ ) = (f ) .

(C.4)

The GW distribution function can be used to define the energy density of GW in the rest frame as energy per unit volume and unit solid angle:

This definition allows us to express the GW density parameter Õ GW (Ê Õ , nÕ ) in the rest frame

Using eq. (C.4), we have the equality

Hence, we find that the GW density parameter in the moving frame S is related with the corresponding quantity in the frame S Õ at rest through the general formula