Improved reconstruction of a stochastic gravitational wave background with LISA
Raphael Flauger, Nikolaos Karnesis, Germano Nardini, Mauro Pieroni, Angelo Ricciardone, Jesús Torrado
TTTK-20-29
Improved reconstruction of a stochastic gravitationalwave background with LISA
Raphael Flauger a , Nikolaos Karnesis b , Germano Nardini c , Mauro Pieroni d , AngeloRicciardone e,f , Jes´us Torrado g (For the LISA Cosmology Working Group) a UC San Diego, Department of Physics, 9500 Gilman Rd, La Jolla, CA, 92093, USA b APC, Universit´e Paris Diderot, CNRS/IN2P3, CEA/Irfu, Observatoire de Paris, Sorbonne ParisCit´e, 10 rue Alice Domont et L´eonie Duquet, 75013 Paris, France c Department of Mathematics and Physics, University of Stavanger, NO-4036 Stavanger, Norway d Blackett Laboratory, Imperial College London, SW7 2AZ, UK e Dipartimento di Fisica e Astronomia “G. Galilei”, Universit`a degli Studi di Padova, via Marzolo8, I-35131 Padova, Italy f INFN, Sezione di Padova, via Marzolo 8, I-35131 Padova, Italy g Institute for Theoretical Particle Physics and Cosmology (TTK), RWTH Aachen University, D-52056 Aachen, Germany
Abstract:
We present a data analysis methodology for a model-independent recon-struction of the spectral shape of a stochastic gravitational wave background with LISA.We improve a previously proposed reconstruction algorithm that relied on a single Time-Delay-Interferometry (TDI) channel by including a complete set of TDI channels. As inthe earlier work, we assume an idealized equilateral configuration. We test the improvedalgorithm with a number of case studies, including reconstruction in the presence of twodifferent astrophysical foreground signals. We find that including additional channels helpsin different ways: it reduces the uncertainties on the reconstruction; it makes the globallikelihood maximization less prone to falling into local extrema; and it efficiently breaksdegeneracies between the signal and the instrumental noise. Project coordinator and corresponding author: [email protected] a r X i v : . [ a s t r o - ph . C O ] S e p ontents A.1 Polarization tensors 26A.2 TDI variables 27A.3 Response functions 28A.4 Noise spectra 30
The
Laser Interferometer Space Antenna (LISA) experiment [1], a European Space Agency(ESA) mission with National Aeronautics and Space Administration (NASA) partnership,will be the milestone for Gravitational Wave (GW) detection in space. It is planned tolaunch in the mid-2030s and it will consist of a constellation of three satellites forming anearly equilateral triangle with 2.5 million km length arms. By monitoring the relativedisplacements among the three satellites, LISA will perform three correlated interferome-try measurements (typically dubbed XYZ channels), which can be transformed into threeuncorrelated data streams (typically dubbed A, E, and T). LISA will probe GWs in the– 1 –illi-Hertz regime which is not accessible with present and future ground-based detec-tors [2–7], and will open a new window for GW physics.The sky is extremely rich in the LISA sensitivity band, with order of ten thousandsources that will be individually resolved during the (at least) four-year mission duration [1].Depending on the astrophysical formation scenarios [8–12], LISA is going to detect fewhundreds of massive ( ∼ –10 M (cid:12) ) black hole binary mergers with Signal-to-Noise Ratio(SNR) up to thousand. These likely constitute the strongest sources that LISA will resolveindividually. Besides them, LISA will resolve between a few to a few hundred extrememass ratio inspirals per year [13] and up to tens of thousands of galactic white dwarfbinaries [14–18].The existence of so many resolvable sources is accompanied by the presence of a largenumber of events that will not be resolved individually, leading to the generation of aStochastic GW Background (SGWB). Given the LISA detection performances, there aretwo guaranteed astrophysical background components in the LISA band: the componentsourced by the unresolved compact Galactic binaries, overcoming the instrumental noiseat the frequencies 2 × − − × − Hz [19, 20], and the component due to neutron starand stellar-origin black hole mergers [21], relevant in the range 2 × − − × − Hz.Beyond these, it is plausible to expect an extra component originated by the extreme massratio inspirals, possibly above the LISA sensitivity at about 10 − − − Hz [22].On top of these astrophysical sources, LISA will potentially be sensitive to SGWBs re-lated to physics of the early universe or physics of particles beyond the standard model [23].(For recent reviews of non-astrophysical sources of SGWB, see e.g. refs. [24, 25]). Com-monly considered examples include superradiance effects [26, 27], cosmological first orderphase transitions [28, 29], networks of topological defects (e.g. cosmic strings) [30] andeven inflationary models [31] with and without subsequent SGWB due to primordial blackholes [32]. In practice, the actual SGWB signal may be a superposition of any of thesesources — and possibly of previously unknown ones — with the SGWB components ofastrophysical origin.In order to extract as much information about the processes sourcing the SGWB aspossible, a detailed signal characterization is required. This relies on an accurate recon-struction of the overall SGWB signal and the capability of precisely breaking it into itscomponents. The Galactic component is expected to be anisotropic and to have a yearlymodulation which helps in (partially) disentangling it from the isotropic components [33]. The astrophysical extra-Galactic component is instead expected to be almost isotropic andstationary, so that separating it from other isotropic and stationary signals is challenging.However, it has recently been shown that it is possible to employ techniques to perform anaccurate separation of this component from other SGWBs [34]. It is worth stressing thatcaution is necessary when moving towards the treatment of realistic LISA SGWB data dueto the complexity of the instrumental noise, the way LISA measures the signals, and thepossibility that an unexpected, unmodeled SGWB contaminates the data. For example, The angle between the Galactic disk and the plane defined by the three LISA’s satellites varies whileLISA orbits around the Sun. This is at origin of the modulation. – 2 –e should expect instrumental noise non-stationarities such as noise transients [35, 36] andspectral lines (as seen already in LISA Pathfinder data [37]), and slow variations of thenoise power spectrum due to thermal gradients inside the spacecrafts [38, 39]. In addition,non-stationarities and/or non-Gaussianity can arise from the non-perfect subtraction ofbright GW transient signals [40]. There are various strategies to mitigate, or take intoaccount these effects, such as adopting higher-tail likelihood functions, or fitting a moresophisticated noise power spectrum model. In this work, we shall assume a simplified sce-nario, where the instrumental noise is Gaussian and stationary, and that we are workingwith perfect residuals, which means that bright sources have been subtracted perfectlyfrom the data-stream.In the same spirit as ref. [41], and following up on ref. [42], we employ a model-independent approach to reconstruct the frequency shape of an unknown SGWB. Ourprocedure is based on and extends the
SGWBinner code [42], which divides the frequencyrange of LISA into sub-intervals ( bins ) in which we the signal is approximated by a powerlaw. By starting with an arbitrary number of initial bins and subsequently combining themaccording to information criteria to avoid over-fitting, our procedure is robustly able tocapture the frequency dependence of the spectrum, and leads to an accurate reconstructionof the injected signal [42]. The analysis carried out in ref. [42] relied on a single data stream(in practice ref. [42] used the X channel). In this work we consider the full set of correlatedchannels X, Y, and Z. As usual, we perform the diagonalization of these data streams anddescribe the signal and noise within the data in the so-called AET basis [43, 44]. As weshow in this work (where some simplified assumption are made), this allows us to breakdegeneracies between the signal and noise spectra which helps in disentangling the twocomponents. We test the model independent reconstruction, using the AET channels, on three wellmotivated SGWB signals: a simple power law that can mimic an astrophysical GW signalcoming from numerous unresolved binary mergers [21, 45] and/or a cosmological signalcoming from some inflationary mechanisms [31] or cosmic strings [30]; a broken power-lawsignal approximating the SGWB from a first order phase transition in the early Uni-verse [28, 29]; and finally a bump signal, which is expected mostly from post-inflationarysources of GWs, like preheating [31]. In the present paper our primary focus is the de-velopment of a model-independent method to reconstruct the overall homogeneous andisotropic SGWB signal, without aspiring to address the issue of the component separation.Nevertheless, in order to prepare a method that remains robust even when backgroundsand foregrounds coexist, we perform some tests on the method’s reconstruction capabilitiesalso in the presence of some illustrative foregrounds. Further improvements which wouldmake the analysis more realistic and could degrade our reconstructions are mentioned. Notice that in log-log scale this is the equivalent to performing a series of Taylor expansions truncatedat linear order. In this analysis we work assuming that the detector arm lengths are equal. This allows to state that inthe AET basis the noise is orthogonal if identical power spectral density in every noise component (e.g. thenoise reduction system and phasemeters) on each spacecraft are considered. The case where the noise isnot orthogonal is described in ref. [44]. – 3 –his paper is structured as follows. In section 2 we present our formalism for modelingthe LISA data stream in terms of a signal and a noise component, and we introduce ourmock signals. In section 3 we explain our data analysis techniques and outline the mainalgorithm of the
SGWBinner code introduced in ref. [42]. In section 4 we present the resultsobtained by applying our techniques to a set of mock signals and in section 5 we test therobustness of the results in the presence of foregrounds. Finally in section 6 we draw ourconclusions. Some technical details are provided in appendix A.
The main observable of interest for the detection and characterization of an isotropic SGWBis its power spectrum P h ( k ). In this section we provide a brief summary on how LISAmeasures the SGWB power spectrum. We begin with the assumption that we are working with ‘perfect’ residuals, which meansthat all transient signals and glitches in the noise have been subtracted from the timestream. After this procedure, the time stream d ( t ) contains noise, which we denote by n ( t ),and a residual stochastic signal s ( t ). We assume that both the noise and the residual signalare stationary. Because of periodic antenna re-pointing and other operational interruptions,the data are expected to be broken up into segments of length T . For the sake of comparisonwith the previous study [42], we take T = 11 . In practice the data stream d ( t ) = s ( t )+ n ( t ) is sampled at a finite rate, but we model itas a real-valued function on the interval [ − T / , T /
2] to keep the notation more appealing.We assume that signal and noise are uncorrelated. We can then treat them separatelyand similarly. First, we define the Fourier transform˜ s ( f ) = (cid:90) T/ − T/ dt e πift s ( t ) . (2.1)For stationary s ( t ), the ensemble average of Fourier modes must obey (cid:104) ˜ s ( f )˜ s ∗ ( f (cid:48) ) (cid:105) ≡ δ (cid:0) f − f (cid:48) (cid:1) S ( f ) , (2.2)for real and positive S ( f ). Furthermore, since s ( t ) is real, the Fourier transform obeys thereality condition ˜ s ∗ ( f ) = ˜ s ( − f ), so that S ( f ) = S ( − f ). For this reason, it is conventionalto work with the “one-side” power spectrum S ( f ) defined only for positive (physical)frequencies, and to include a factor 1 / The value of T is still under discussion, and in practice not all segments may have the same duration.In the stationary regime, the duration T of the segments is fictitious. We however prefer to keep this formatto prepare the formalism and the code for future developments where the stationary assumption between asegment and the subsequent one is broken. – 4 –ill provide three TDI channels, called { X,Y,Z } or { A,E,T } , depending on the basisadopted (see sec. 2.3 for the relationship between the two bases). Equation (2.2) thengeneralizes to (cid:104) ˜ s i ( f )˜ s ∗ j ( f (cid:48) ) (cid:105) ≡ δ (cid:0) f − f (cid:48) (cid:1) S ij ( f ) , (2.3)where the latin indices run over the basis { X,Y,Z } or { A,E,T } , and the spectral densities S ij ( f ) are characterized by a Hermitian matrix-valued function that is related to the powerspectrum P h ( f ) through the respective response functions R ij : S ij ( f ) = R ij ( f ) P h ( f ) . (2.4)The response functions for the different channel cross-spectra are shown in section 2.2, andtheir derivation is provided in appendix A. Notice that eq. (2.4) quotes the relationshipin the limit of negligible relative motion of the LISA’s spacecrafts. For our purpose, thematrix R ij ( f ) is real and symmetric (see e.g. ref. [46] for some peculiar scenarios wherethis property does not occur).Likewise, for stationary (and real) noise n ( t ) one concludes (cid:104) ˜ n i ( f )˜ n ∗ j ( f (cid:48) ) (cid:105) = 12 δ (cid:0) f − f (cid:48) (cid:1) N ij ( f ) , (2.5)where N ij ( f ) is a Hermitian matrix function of single-sided noise power spectra . For the SGWB of primordial origin, it is common practice to predict the signal interms of Ω GW ( f ), the energy density per logarithmic frequency interval scaled by thecritical density. The noise spectrum in Ω GW ( f ) is related to the noise spectrum in eq. (2.5)by N Ω ij ( f ) = 4 π f H /h ) N ij ( f ) , (2.6)with H denoting today’s Hubble parameter whose dimensionless parameter normalizationis h . The power spectrum of the signal can be recast in units of Ω GW ( f ) in the same way. So far our discussion of signal and noise power spectral densities generally applies to anygravitational wave interferometer. We now specialize our discussion to LISA, beginningwith the noise model. The TDI variables are designed to eliminate the dominant sourceof noise caused by fluctuations in the central frequency of the laser as well as the noisecaused by displacements of the optical benches. In this simplified model, the residualnoise components that enter into each TDI channel can be grouped into two effectivequantities, dubbed “Interferometry Metrology System” (IMS) noise, which, for example,includes shot noise, and “acceleration” noise associated with the random displacements ofthe proof masses caused, for example, by local environmental disturbances. Our current In the literature the functions N ij ( f ) are often denoted by P n . Their dimension is Hz − since both˜ n ( f ) and δ have dimension Hz − . – 5 –nderstanding of the LISA noise is based on the LISA Pathfinder experiment [37] andlaboratory tests. The IMS and acceleration noise power spectra are given by [47] P IMS ( f, P ) = P pm Hz (cid:34) (cid:18) f (cid:19) (cid:35) (cid:18) πfc (cid:19) ,P acc ( f, A ) = A fm s Hz (cid:34) (cid:18) . f (cid:19) (cid:35) (cid:34) (cid:18) f (cid:19) (cid:35) (cid:18) πf (cid:19) (cid:18) πfc (cid:19) . (2.7) Figure 1 : Left panel: the IMS and acceleration noise power spectra expressed in eq. (2.7).Right panel: The IMS and acceleration contributions to N TT weighted with their pre-factors as in eq. (2.20). In both panels we fix P = 15 and A = 3.ESA’s mission specifications require the amplitudes to be P = 15 and A = 3 with ± P and A .The precise determination of the noise properties is one of the main technical challengesof the LISA mission, and it seems premature to attempt to keep track of the full complexityhere. For the IMS contribution, we make the simplifying assumption that the noise spectrafor all links are identical, stationary, and uncorrelated. Similarly, for the acceleration noise,we assume that the fluctuations of the masses are isotropic and stationary, that the powerspectra for all test masses are equal, and that the fluctuations of the different massesare uncorrelated. Furthermore, we assume that the three satellites form an equilateraltriangle, L = L = L = L = 2 . × m. Under these assumptions, the total powerspectral density for the noise auto-correlation is N aa ( f, A, P ) = 16 sin (cid:18) πf Lc (cid:19) (cid:26)(cid:20) (cid:18) πf Lc (cid:19)(cid:21) P acc ( f, A ) + P IMS ( f, P ) (cid:27) , (2.8)and the noise cross-spectra are N ab ( f, A, P ) = − (cid:18) πf Lc (cid:19) cos (cid:18) πf Lc (cid:19) [4 P acc ( f, A ) + P IMS ( f, P )] , (2.9)where a, b ∈ { X,Y,Z } and a (cid:54) = b . Notice, in particular, that for our assumptions, the noisecovariance matrix is real. For completeness, we include a derivation in appendix A.4.– 6 –s we explain below, we marginalize over the amplitude of the IMS and accelerationnoise power spectra in our analysis. However, we take the functional form of the noisemodel used to generate the mock data to be the same as the model used in fitting thedata. So differences between the functional form of the instrumental noise and the noisemodel would introduce a bias and will have to be closely monitored. This is particularlytrue for LISA, for which calibration of the noise must happen at the same time as thereconstruction of the SGWB and other signals. So far we have defined the signal and noise spectra in section 2.1 in an arbitrary basis, andhave introduced the LISA noise spectra in the XYZ basis in section 2.2. We now introduceanother commonly used basis of TDI channels, the AET basis, which diagonalizes thesignal and noise covariance matrices [43, 44]. As we saw in section 2.1, for a stationaryand isotropic SGWB signal and stationary noise (that is uncorrelated with the signal), theauto- and cross-spectra of the different channels read (cid:104) ˜ d i ˜ d ∗ j (cid:105) (cid:48) = R ij P h ( f ) + N ij ( f ) , (2.10)where the prime indicates that we have stripped the frequency δ -function and the factor1 /
2, and we have used eq. (2.2) and eq. (2.5). Under the assumptions made in section 2.2the three interferometers have identical properties, and the noise spectra and responsefunctions obey N XX = N YY = N ZZ , N XY = N YZ = N XZ , (2.11) R XX = R YY = R ZZ , R XY = R YZ = R XZ . (2.12)The resulting symmetry properties of (cid:104) ˜ d i ˜ d ∗ j (cid:105) (cid:48) imply that it can be diagonalized by˜ d A = 1 √ (cid:16) ˜ d Z − ˜ d X (cid:17) , ˜ d E = 1 √ (cid:16) ˜ d X − d Y + ˜ d Z (cid:17) , ˜ d T = 1 √ (cid:16) ˜ d X + ˜ d Y + ˜ d Z (cid:17) , (2.13)or any other transformation related by a cyclic permutation in X, Y, Z, where ˜ d A , ˜ d E , ˜ d T and˜ d X , ˜ d Y , ˜ d Z are the data in the AET and XYZ streams, respectively, and we have suppressedthe argument to streamline our notation.The spectra in the AET basis become (cid:104) ˜ d A ˜ d ∗ A (cid:105) (cid:48) = (cid:104) ˜ d E ˜ d ∗ E (cid:105) (cid:48) = ( R XX − R XY ) P h ( f ) + N XX − N XY , (2.14) (cid:104) ˜ d T ˜ d ∗ T (cid:105) (cid:48) = ( R XX + 2 R XY ) P h ( f ) + N XX + 2 N XY , (2.15) (cid:104) ˜ d A ˜ d ∗ E (cid:105) (cid:48) = (cid:104) ˜ d A ˜ d ∗ T (cid:105) (cid:48) = (cid:104) ˜ d E ˜ d ∗ T (cid:105) (cid:48) = 0 . (2.16)– 7 – igure 2 : Left panel: LISA noise spectra in the XYZ and AET basis. Right panel: squareroot of the LISA strain noise for the different XYZ/AET combinations.We see that the noise spectra and response functions in the AET basis are given by N AA = N EE = N XX − N XY , N TT = N XX + 2 N XY , (2.17) R AA = R EE = R XX − R XY , R TT = R XX + 2 R XY . (2.18)Substituting the noise auto- and cross-spectra, eq. (2.8) and eq. (2.9) into eq. (2.17) weobtain N AA ( f, A, P ) = N EE ( f, A, P ) == 8 sin (cid:18) πf Lc (cid:19) (cid:26) (cid:20) (cid:18) πf Lc (cid:19) + cos (cid:18) πf Lc (cid:19)(cid:21) P acc ( f, A ) ++ (cid:20) (cid:18) πf Lc (cid:19)(cid:21) P IMS ( f, P ) (cid:27) , (2.19)and N TT ( f, A, P ) = 16 sin (cid:18) πf Lc (cid:19) (cid:40) (cid:20) − cos (cid:18) πf Lc (cid:19)(cid:21) P acc ( f, A ) ++ (cid:20) − cos (cid:18) πf Lc (cid:19)(cid:21) P IMS ( f, P ) (cid:27) . (2.20)These power spectrum densities are displayed in the left panel of fig. 2.Concerning the measurement of the noise amplitude parameters, from the left panelof fig. 1 one might naively expect that data at f (cid:46) × − contain precise information onthe A parameter, and in particular, to avoid signal contamination, A would be measuredmost efficiently via the T channel. However, the low-frequency expansion of eq. (2.20)yields N TT ( f, A, P ) (cid:39) (cid:18) ff ∗ (cid:19) sin (cid:18) ff ∗ (cid:19) (cid:34)(cid:18) ff ∗ (cid:19) P acc ( f, A ) + P IMS ( f, P ) (cid:35) , (2.21)– 8 –ith f ∗ ≡ (2 πL/c ) − (cid:39) .
019 Hz. So, at low frequencies, the contribution coming from P acc is suppressed by a factor ( f /f ∗ ) with respect to the contribution coming from P IMS , whichmeans that the shape of N TT is dominated by P IMS even at low frequencies , as can be seenin the right panel of fig. 1. Thus, despite the noise-only character of the TT channel, itslow frequency band only provides minimal information on A.
Figure 3 : Left panel: geometrical contribution to the LISA response function. Rightpanel: full LISA response function.As we discuss in more detail in appendix A.3, the quadratic response functions for ourchoice of TDI variables can be written as R ij ( f ) = 16 sin (cid:18) πf Lc (cid:19) (cid:18) πf Lc (cid:19) ˜ R ij ( f ) . (2.22)The first factor reflects the choice of TDI variables, the factor 2 πf L/c appears becauseLISA measures frequency changes rather than time delays (see appendix A.1 for details),and ˜ R ij is a factor that depends on the geometry of the detector. These geometrical factorsand the corresponding response functions are shown in fig. 3.In practice, we work with the numerical expressions shown in fig. 3, but for somepurposes simple analytic approximations are helpful. For the XX channel (as well as theYY and ZZ channels), the geometrical factor is well approximated by [48]˜ R XX ( f ) = 310 11 + 0 . (cid:16) πfLc (cid:17) . (2.23)Similarly, for the AA channel (as well as the EE channel) and for the TT channel, we find˜ R AA ( f ) = 920 11 + 0 . (cid:16) πfLc (cid:17) , (2.24)˜ R TT ( f ) = 920 (cid:16) πfLc (cid:17) . × + 0 . (cid:16) πfLc (cid:17) . (2.25)– 9 – igure 4 : Left panel: ratios between the strain sensitivities for the different channels andthe effective strain defined in eq. (2.31). Right panel: LISA sensitivity (see eq. (2.27)) forthe different channels and for their inverse variance weighted combination.From the power spectrum density and the response functions we can define the strainsensitivities S n ij ( f, A, P ) = N ij ( f, A, P ) R ij ( f ) = N ij ( f, A, P )16 sin (cid:16) πfLc (cid:17) (cid:16) πfLc (cid:17) ˜ R ij ( f ) . (2.26)Plots of these quantities are shown in the right panel of fig. 2.Putting everything together, we arrive at the noise model used for the generation ofthe noise realizations and extraction of noise parameters used in our analyses, in units ofthe energy density parameterΩ n,ij h ( f, A, P ) ≡ π f H /h ) S ijn ( f, A, P ) , (2.27)where H /h (cid:39) . × − Since both the signal and noise covariance in the AET basis are diagonal, it is also straight-forward to form an optimal (unbiased and with minimal variance) power spectrum estima-tor from the three channels [49]. Reserving the lower case latin indices for formulae thatapply for both the XYZ and the AET basis, and introducing greek indices for the AETbasis, α ∈ { A , E , T } , we can define the time stream variables ˆ d α = ˜ d α / √R αα , with unitresponse to the SGWB signal (cid:104) ˆ d α ˆ d ∗ α (cid:105) (cid:48) = P h ( f ) + S n α , (2.28)where S n α = N αα /R αα is the effective noise power spectral density or sensitivity. Interms of these variables the unbiased power spectrum estimator for a single channel isˆ P α = | ˆ d α | − S n α , and we can define an unbiased estimator by combining the channelsˆ P = (cid:88) α w α ˆ P α , with (cid:88) α w α = 1 . (2.29)– 10 –s usual, the optimal linear combination is obtained by inverse variance weighting w α = ( P h + S n α ) − (cid:80) α ( P h + S n α ) − . (2.30)By calculating the variance of this estimator, we see that in the limit of small signal-to-noiseratio per mode, the effective noise power spectral density for this combination of channelsis given by S n = 1 (cid:114)(cid:80) α S − n α . (2.31)A comparison between this quantity and the standard noise spectra is presented in fig. 4.The effective sensitivity for the optimal combination provides a simple way to comparesensitivities for a single channel as used in ref. [42] and three channels. This is shownin fig. 4. Beyond this comparison, we will not work with this combination and will continueto use the XYZ and AET bases, which, as we will see, has practical advantages. However,this combination would be a natural choice to visually present the results of a measurementof a SGWB with LISA. In this section we discuss the data analysis techniques employed in this work. As usual,we compute the posterior probability for the model parameters using the Bayes theorem: p ( (cid:126)θ, (cid:126)n | D ) = π S ( (cid:126)θ ) π N ( (cid:126)n ) L ( D | (cid:126)θ, (cid:126)n )Z N,S ( D ) , (3.1)where L ( D | (cid:126)θ, (cid:126)n ) is the likelihood of the data D , π N and π S are the prior for the noiseparameters (cid:126)n and for the signal parameters (cid:126)θ , respectively, and Z N,S ( D ) is the modelevidence. Since we are not interested in the model evidence (or the normalization ofthe posterior probability in general), we work with the unnormalized posterior P ( (cid:126)θ, (cid:126)n | D ).Omitting D for simplicity, we have P ( (cid:126)θ, (cid:126)n ) = π S ( (cid:126)θ ) π N ( (cid:126)n ) L ( D | (cid:126)θ, (cid:126)n ) ∝ p ( (cid:126)θ, (cid:126)n | D ) . (3.2)In this work (consistent with the analysis in ref. [42]) we assume a Gaussian prior for thenoise parameters: (cid:126)n = ( A, P ) ∼ N ( (cid:126)µ, Σ ) , with (cid:126)µ = ( ¯ A, ¯ P ) = (3 , , and Σ = diag(0 . (cid:126)µ ) . (3.3)The prior on the signal parameter is model-dependent, but in all subsequent analyses in thispaper will be taken as uniform (or log-uniform, for parameters representing amplitudes).The Maximum A Posteriori (MAP) parameters (denoted with (cid:126)θ b and (cid:126)n b ) are definedas ∂ j ln P ( (cid:126)θ, (cid:126)n ) (cid:12)(cid:12)(cid:12) (cid:126)θ b ,(cid:126)n b ≡ , (3.4)– 11 –here the j index defines one equation per parameter, both signal and noise ones. As it iscustomary, the Fisher information matrix is defined as I ij = C − ij ≡ −(cid:104) ∂ i ∂ j ln P ( (cid:126)θ, (cid:126)n ) (cid:12)(cid:12)(cid:12) (cid:126)θ b ,(cid:126)n b (cid:105) . (3.5)Notice that with this definition the Fisher matrix also contains information about thepriors. For our simulations, we will assume that the data is taken during a full mission durationof 4 years with 75% efficiency, giving an effective observation time of 3 years. As discussedin section 2.1, we will assume the TDI data to be divided into roughly 95 chunks (ordata segments) of 11 . × − Hz to 5 × − Hz, thisamounts to around 5 × data points per frequency, with a frequency resolution of ≈ − Hz, for a total of roughly 5 × data points.For each chunk l and channels i, j , we generate the data D lij ( f ) = d li ( f ) d l ∗ j ( f ) directlyin the frequency domain, following the procedure described in ref. [42]. Generation in thefrequency domain allows us to ignore window effects and overlapping segments, simplifyingthe process.We compress the simulated data by averaging over the chunks and binning the data.Specifically, we define the average¯ D ij ( f ) ≡ N c N c (cid:88) l =1 D lij ( f ) , (3.6)where N c = 95 is the number of data segments. We then bin the data in the frequencydomain into a new data set (cid:16) f ( k ) ij , D ( k ) ij (cid:17) . In practice, we use inverse variance weighting f ( k ) ij ≡ (cid:88) f ∈ { f ( k ) } w ( k ) ij ( f ) f , (3.7) D ( k ) ij ≡ (cid:88) f ∈ { f ( k ) } w ( k ) ij ( f ) ¯ D ij ( f ) , (3.8)where (cid:8) f ( k ) (cid:9) are the (discrete) frequencies in bin k , and the weights are w ( k ) ij ( f ) = σ − ij ( f ) (cid:80) f ∈ { f ( k ) } σ − ij ( f ) , (3.9)where σ ij ( f ) is an estimate of the variance of the segment-averaged data ¯ D ij ( f ). For theanalyses presented in this work we assume the variance to be well approximated by LISA’ssensitivity introduced in eq. (2.27). Since the variance depends on the combination of It can be shown that this procedure would be optimal in the case of noise-dominated data with realnoise parameters equal to their face values. In the case of a SGWB detection with signal-to-noise ratiocomparable to or larger than unity per mode, we would have to update the variance of the data to accountfor it. – 12 –hannels, the inverse variance weighting implies that the discrete frequencies also dependon the channel combination.In practice, the second step is performed by requiring a maximum of 1000 linearlyspaced points per frequency decade. For different bins the coarse-grained points are ob-tained by summing over different numbers of data points. To keep track of the weightscarried by the different coarse-grained data points, we denote as n ( k ) ij the number of pointswithin bin k for the cross-spectrum of channels i and j .As a starting point to build a likelihood for the data, it is tempting to assume a simpleGaussian likelihoodln L G ( D | (cid:126)θ, (cid:126)n ) = − N c (cid:88) i,j (cid:88) k n ( k ) ij (cid:34) D thij ( f ( k ) ij , (cid:126)θ, (cid:126)n ) − D ( k ) ij D thij ( f ( k ) ij , (cid:126)θ, (cid:126)n ) (cid:35) , (3.10)where the indices i, j run over the different channel combinations, the index k runs overthe coarse-grained data points, and we have defined D thij ( f, (cid:126)θ, (cid:126)n ) ≡ R ij Ω GW h ( f, (cid:126)θ ) +Ω n,ij h ( f, (cid:126)n ) as our model for the data. Notice that since in the XYZ basis the chan-nels are correlated, D ij is not diagonal in channel space. On the other hand in the AETbasis this likelihood reduces to the sum of the diagonal elements in channel space.However, it is known that the Gaussian likelihood in eq. (3.10) would give a result whichis systematically biased [50–53], due to the mild non-Gaussianity of the full likelihood ofthe non-averaged data. In order to correct for this bias we also introduce a log-normallikelihood [50, 51]ln L LN ( D | (cid:126)θ, (cid:126)n ) = − N c (cid:88) i,j (cid:88) k n ( k ) ij ln (cid:34) D thij ( f ( k ) ij , (cid:126)θ, (cid:126)n ) D ( k ) ij (cid:35) , (3.11)and define our final likelihood as [52]ln L = 13 ln L G + 23 ln L LN . (3.12)It can be shown [52] that this likelihood accounts for the skewness in the full non-Gaussianlikelihood, giving a more accurate result for the model parameters.As we just saw, the likelihood depends on a theoretical spectrum. This could be a modelfor one or several sources of a SGWB. Here we instead employ a non-parametric approachwhich aims to reconstruct the frequency dependence of the SGWB without assuming anyprevious knowledge of its spectral shape. We define a piece-wise model for the signal ona set of frequency intervals (bins). Within each of these bins the signal is assumed to bewell approximated by a power law: h Ω GW (cid:16) f, (cid:126)θ i (cid:17) = (cid:88) i α i (cid:32) f (cid:112) f min ,i f max ,i (cid:33) n t,i Θ ( f − f min ,i ) Θ ( f max ,i − f ) , (3.13)where Θ is the Heaviside step function, the index i runs over the bins, f min , i and f max , i arerespectively the minimum and maximum frequencies in bin i , f ∗ ,i ≡ (cid:112) f min ,i f max ,i is thepivot frequency for bin i and (cid:126)θ i = { α i , n t,i } are the signal parameters for bin i (respectivelydenoting amplitude at the pivot, and tilt). As we discuss below, the number of bins willbe adjusted dynamically. – 13 – .2 Algorithm for binned reconstruction of a SGWB signal Here we present an update of the algorithm described in ref. [42], and implemented in the
SGWBinner code, to reconstruct the frequency shape of the SGWB, using the binned power-law signal model and the noise model presented in section 2.2. Given some pre-determinedbinning, this problem reduces to finding the solution of eq. (3.4) over all parameters si-multaneously (the two signal parameters per bin, and the two noise parameters which arecommon to all the bins) together with some estimation of the uncertainty of the recon-struction.This non-trivial minimization problem may be tackled with Monte Carlo (MC) meth-ods. However, if we wish to leave the number of bins as a free parameter, a MC samplewould have to be obtained for each binning configuration, and then the optimal numberof bins would be chosen by e.g. comparing Bayesian evidences of different bin numbers.This process is computationally expensive. For this reason, similarly to ref. [42], we adopta different strategy:1. We build a prior for the noise parameters using the TT-channel (more on this below).2. We bin the signal in frequency in a given initial number of bins, which will be themaximum allowed for this run.3. We minimise the posterior for channel AA independently for each bin, imposing theTT-defined noise prior, for the two independent signal and two independent noiseparameters.4. For all pairs of neighboring bins, we iteratively check whether merging the two binsis statistically favoured, using the Akaike Information Criterion [54] as described inref. [42]. This procedure reduces the chance of overfitting, improving the quality ofthe reconstruction. If one or more merges are deemed favourable, we re-define thebins and go back to the minimisation step.5. When the number of bins has converged, we estimate the error on the reconstructionand provide a visualization of the 1 σ and 2 σ regions for the reconstructed signal,using the procedure described in ref.[42].The advantage of this methodology is that it greatly simplifies the problem from the nu-merical point of view, making it easily solvable in a reasonable amount of computationtime. At this point we should stress that step 1 is crucial for two main reasons: • A wrong estimate of the noise parameters may bias the signal reconstruction. • Since the noise parameters are actually common to all bins, this procedure may onlybe trusted if the results obtained in the per-bin analysis are consistent among allbins.In ref. [42], in order to define suiting priors for the noise parameters, we took advantageof the morphology of the spectra of the data. Based on the reasonable assumption that– 14 –t very low and very high frequencies the sensitivity of the detector is significantly weakerthan at the central O (10 − ) Hz region, we used these outer frequency regions to definestrong priors on the instrument noise for the interval of the frequency range in which weattempt to reconstruct the signal.That procedure provided good results when working with a single TDI channel. Inthe present work, however, we work with the noise-orthogonal A, E, and T TDI combi-nations [49]. As explained in section 2.2, at low frequencies the TT channel combinationwould greatly suppress any GW signal, and would therefore allow us to use it for a defini-tion of the instrument noise prior density, without the need to limit the frequency rangeon which we reconstruct the background signal. In practice, we find the posterior in theTT channel combination of a power-law signal and the noise parameters, and build a prioron the P noise parameter (the best-determined one in TT) by marginalizing over the rest.This is one of the main improvements of the present analysis with respect to the one ofref. [42].As an optional final step in the algorithm, once the optimal number of bins has beenset, we may then run an MC sampler on the total posterior for all bins and all channelcombinations. We have implemented this step in the SGWBinner code by interfacing withthe
Cobaya
MC framework [55]. One would expect this procedure to provide a moreaccurate estimation of the uncertainty on the reconstruction, but in practice we foundlittle difference with the results of the original algorithm, as discussed in section 4.5.
In this section we apply the extension of the
SGWBinner algorithm described in the previoussection to a number of mock data sets generated as described in section 3.1 in whichwe have injected large-amplitude SGWB signals. The aim of this section is to asses thecapability of our pipeline to perform a model-independent reconstruction of SGWB signalswith distinct frequency dependence and make a comparison with the reconstruction usingone LISA single-channel that was done in ref. [42]. In all the examples in this section, wehave started from log-spaced bins (as opposed to equal-SNR bins).
Our first benchmark signal is a simple power-law described by h Ω GW ( f ) = 10 α ∗ (cid:18) f .
001 Hz (cid:19) n t (4.1)where we have chosen an amplitude α = − .
352 at the pivot frequency f ∗ = 0 .
001 Hz anda tilt n t = 2 /
3. A power-law spectrum with a similar slope can be used to approximatesome inflationary models [31] as well as the astrophysical SGWB from stellar-origin blackhole binaries and neutron star binaries [21] (see section 5 for more details).The result of the reconstruction algorithm is shown in fig. 5. Due to the smooth natureof the signal, we started with a small number of bins, 10, which, as expected being the signala smooth power law across the whole LISA band, converged to a single bin and produced a– 15 – AP Figure 5 : Reconstruction of the power-law signal described in section 4.1 in the AETchannels using 10 initial bins, which converged to a single one. On the bottom right, weshow the contour plot of the posterior of the parameters of the single bin, including theMaximum A Posteriori (MAP).successful reconstruction. We also show the 1 σ and 2 σ posterior contours of the amplitudeand slope of the reconstructed signal in the single bin. When translated from the pivot scaleof the power law to the pivot of the single bin, the amplitude and the slope fall comfortablywithin the 1- σ C.L. contour. The SNR of the reconstructed signal is approximately 483.The reconstructed signal (yellow line) matches the input signal (dashed-dot blue) very wellwithin the error bars quantified by the thickness of the light-blue curves, as can be seen inthe zoomed region around the pivot in fig. 5.For the power-law case we test the procedure also in the X channel for a comparison.The result of the reconstruction with the same input parameters is given in fig. 6. In thiscase we also recover the power-law parameters at 1- σ , and the SNR of the reconstructedsignal is approximately 474. The reconstructed signal here is also very well within theerror bars of the reconstruction. The error bar of a reconstruction performed using thethree TDI 1.5 channels would be ∼ √ Our second benchmark signal is a broken power law with the following functional form: h Ω GW ( f ) = 10 α (cid:18) ff ∗ (cid:19) n t (cid:34) (cid:18) ff ∗ (cid:19) n t − n t (cid:35) (4.2)evaluated at α = − f ∗ = 10 − Hz, with slopes n t = 3 and n t = −
3. This kindof spectral shape might arise from stochastic GW produced during a first order phasetransition around the TeV energy scale (see e.g. refs.[28, 29] for a dedicated analysis for– 16 – igure 6 : Reconstruction of the power-law signal in fig. 5 using X-channel data only.Notice that there is no reconstruction in the outer bins, which must be used to estimatethe noise.LISA). The result of the reconstruction is shown in fig. 7. We have started from 20 initialbins and converged to 8. The reconstructed signal in this case has a large SNR (cid:39) . × .Despite the peculiar shape of the signal, the reconstruction is very accurate, also in theexternal bins, as can be inferred from the small size of the C.L. contours, which are shownin fig. 7b. As a last benchmark model we have chosen a bump signal characterized by an amplitude A ∗ and a width ∆ according to h Ω GW ( f ) = A ∗ exp (cid:26) − [log ( f /f ∗ )] ∆ (cid:27) . (4.3)This kind of signal is expected from cosmological sources such as non-perturbative effectsduring post-inflationary preheating, strong first order phase transitions during the thermalera of the universe, or merging of PBH’s during the early universe. These phenomenatypically generate single- or multi-peaked spectral signals that can be described as a bump.As input parameters for the signal we have chosen { A . = 10 − . , ∆ = 0 . } .The result of the reconstruction procedure is shown in fig. 8. We have produced itchoosing 60 initial bins, which has converged to 7. The reconstruction is more accurate inthe central part of the frequency band, where the signal peaks, and worsens in the outerbins. The SNR of the reconstructed bump signal is approximately 339. One might naively expect that the uncertainties of the parameter reconstruction scaleroughly with the square root of the number of independent channels we consider. This guess– 17 – Frequency [Hz]10 h G W Data (used by the binner)LISA SciRDLISA PLS 4.0y, 0.75% eff, SNR=10Input signalReconstructed sensitivityReconstructed signalBin extremesSignal 1 regionSignal 2 region (a) log of the amplitude at f * )1.961.982.002.022.042.062.08 n t ( s p e c t r a l t il t ) Bin 2 ( f * = 1.87×10 Hz)
1 Fisher2 Fisher1 region2 region
MAP log of the amplitude at f * )0.970.980.991.001.011.021.031.04 n t ( s p e c t r a l t il t ) Bin 4 ( f * = 4.93×10 Hz)
1 Fisher2 Fisher1 region2 region
MAP
MAP log of the amplitude at f * )1.561.581.601.621.64 n t ( s p e c t r a l t il t ) Bin 3 ( f * = 3.03×10 Hz)
1 Fisher2 Fisher1 region2 region
MAP
MAP log of the amplitude at f * )1.041.021.000.980.960.94 n t ( s p e c t r a l t il t ) Bin 7 ( f * = 2.12×10 Hz)
1 Fisher2 Fisher1 region2 region
MAP (b)
Figure 7 : (a) Reconstruction of the broken power-law signal described in section 4.2 inthe AET channels using 20 initial bins, which converged to 8. (b)
Contour plots of theposterior of the parameters of the inner bins, showing the maximum a posteriori (MAP).underestimates the gain in constraining power obtained from a multiple-channel analysis.A key advantage of our three-channel analysis is indeed its capability to break degeneracies.Figure 9 clearly illustrates this aspect. It considers the (peculiar but didactic) scenario inwhich in a given channel, say X, the sum of the injected signal and noise matches thenoise curve with some acceleration and interferometry-metrology-system parameters largerthan the nominal ones. Had we exclusively analysed the XX data, the separation betweenthe signal and the noise would have been completely dependent on our prior on the noiseparameters. In addition, despite the huge SNR (calculated with respect to the nominalsensitivity curve), the errors on the parameters would have been large. On the contrary,thanks to the difference between response functions in the different channel combinations of– 18 – igure 8 : Reconstruction of the bump signal described in section 4.3 in the AET channelsusing 60 initial bins, which converged to 7. In the zoomed regions we can see how thereconstruction accommodates the injected signal within the 1- σ band. Figure 9 : Reconstruction of a signal which is degenerate with noise in XX. Breaking thedegeneracy between signal and noise in this case is only possible in a 3-channel analysis(see main text).the AET basis, no signal can be fully degenerate with the noise curve in every channel. Thethree-channel analysis exploits this feature, and, as can be seen in fig. 9, the degeneracy isbroken and the reconstruction is reasonably satisfactory.
As explained in section 3.2, in this work we implement a full-posterior MC as a final step ofthe algorithm, once the number of bins has been fixed. We expect this step to produce more– 19 – Frequency [Hz]10 h G W Data (used by the binner)LISA SciRDLISA PLS 4.0y, 0.75% eff, SNR=10Input signalReconstructed sensitivityReconstructed signalBin extremesSignal 1 regionSignal 2 region
Figure 10 : Alternative reconstruction of the broken power law in fig. 7, sampling fromthe full posterior after having fixed the number of bins. The main effect of considering thefull posterior is the 1- and 2- σ bands in the outer bins being slightly thinner, due to thesmall increase in constraining power given by fitting a single noise model to all bins.accurate uncertainty bands for the reconstruction, since, contrary to the main algorithm,it fits to all bins a common noise model.In practice, we found the effect of the per-bin fit of the noise model to affect thereconstruction uncertainties only mildly, at least for moderately high signal-to-noise SGWBshapes (the only ones for which a spectrum reconstruction is expected to produce significantresults versus a simple model fit). This is possibly due to the null effect of the noise inthe central, high signal-to-noise bins, and to the expected lack of correlation between theleft- and rightmost bins, where the noise is the dominant contribution, since the noise ineach of these outer bins is affected by only one of the noise parameters: A for the leftmost,low-frequency one, and P for the rightmost, high-frequency one. As en example, in fig. 10we show the equivalent result to fig. 7. In this case, as expected, the noise paramters arerecovered at 1- σ . In this section we study the capability of LISA to disentangle different components con-tributing to the SGWB. In particular, we show that the techniques for the model-independentreconstruction described in the previous sections can be extended to the case in which thesignal is a superposition of at least two components, one of which has no reliable model.For this purpose, we consider two mock cases in which the SGWB observed by LISA is thesum of an isotropic astrophysical component (in the following foreground ) and an unex-pected cosmological component. While the latter has no signal model, the spectral shapeof the former is assumed to be well known. We consider two types of foregrounds, each The current models of the foregrounds for LISA still suffer of large uncertainties. It is however reasonableto expect that the foreground models will be quite precise at the end of the LISA mission when, among – 20 –f them expected to capture the qualitative features of some foregrounds really foreseen inLISA:
Extra-galactic foreground:
The incoherent superposition of all the extra-galactic com-pact object mergers that LISA will be unable to resolve, introduce a foreground inthe LISA data. Stellar-origin black hole and neutron star binaries should be themain contributors to this extra-galactic astrophysical signal [45]. For a LIGO-likesensitivity, the foreground of these binaries has a power spectrum (in units of en-ergy density Ω GW ) behaving like a power law with power index 2/3 and amplitudeΩ GW = 8 . +12 . − . × − at f = 25 Hz [21]. Extrapolating this information from theLIGO to the LISA frequency band leads to the foreground model [21, 45] h Ω GW = 10 α FG (cid:18) ff ∗ (cid:19) / . (5.1)From the mean of the above constraint at LIGO frequencies, we obtain α ∗ F G (cid:39) − . f ∗ = 10 − Hz, while α − F G (cid:39) − .
72 from the lower limit. We use these values tobuild the Gaussian prior π F G ( α F G ) ∼ N ( α ∗ F G , σ ) with σ (cid:39) . Galactic (averaged) foreground:
The unresolved sub-threshold mergers of galactic bi-naries constitute a foregrounds with a yearly modulation. A proper treatment of sucha signal should take into account the variation in each chunk, and consequently, asshown in ref. [33], the periodicity information can be exploited to partially subtractthe foreground. For our illustrative purposes, here we do not follow this approachbut consider the signal averaged over the year. This eventually looks like an isotropicand stationary signal with model that in first approximation reads [48] h Ω GW = 10 α FG f / e − a f + a f sin( a f ) { a ( f k − f )] } . (5.2)For concreteness, we take the values quoted in Table 1 of ref. [48] for the parameters a , . . . , a , f k , choosing four years of observation time. We assume these parametersto be known so that no fit on them will be performed. The fit is instead carriedout over the amplitude parametrised by α F G . As a toy example, we build the prioron α F G following some reasonable but somehow arbitrary criteria. We assume theprior to be Gaussian and, following ref. [48], we take it to peak at α ∗ F G = − .
95. Toset the expected deviation interval, we use an uncertainty similar to the one of thesensitivity curve. This roughly corresponds to impose α F G ∈ [ − . , − .
79] at 68%C.L. . The Gaussian prior then becomes π F G ( α F G ) ∼ N ( α ∗ F G , σ ) with σ (cid:39) . other things, the population of the sources will be bound by the measurements of the resolvable events. Several subtleties will have to be addressed to make this model more accurate. It is for instance notobvious that the different resolution power of LISA and LIGO play no relevant role, or that the binaryeccentricity evolution is a minor effect. Notice that the 20% margins on the two noise parameters induce a 45% variation on the level of thesensitivity curve. The values of σ that follows, makes the prior slightly broader than the estimate prediction.This choice is motivated by the large uncertainties in the prior definition, but should not play a big role inthis proof-of-concept example. – 21 –s for the cosmological component, we inject a broken power-law signal, as describedin eq. (4.2), with amplitude α = − f ∗ = 10 − Hz, and slopes n t = 6 and n t = −
4. Asstated, the model for the foreground is assumed to be known, so we are not interested inreconstructing its shape but only its amplitude. For this reason, the foreground parametercan in practice be treated as an extra noise parameter. Following the procedure describedin section 4.5, after we have determined the final number of bins, we perform a final MCsampling in which the foreground amplitude is fit simultaneously for all bins, as the noiseparameters are. In the examples below we have used
PolyChord [56, 57] as the MC sampler,via its interface with
Cobaya [55]. The results were analyzed and the contour plots wereproduced using
GetDist [58].In fig. 11 we show the reconstruction of the broken power-law SGWB signal describedabove in the presence of the extra-galactic foreground component described above. As canbe seen in fig. 11a, wherever the injected signal is sufficiently large with respect to theLISA sensitivity, the signal reconstruction is practically unaffected by the presence of theforeground. Contour plots for the parameters of the binned model are shown in fig. 11b.Except for the two outer bins, where the amplitude and the tilt of the signal are degenerate,the measurements of the log amplitude and of the tilt are consistently accurate. A triangleplot for the two LISA noise parameters together with the log of the foreground amplitudeis shown in fig. 11c. Only mild degeneracy are found and all the parameters are consistentwith the injected values at 2 σ . Notice that the measurement of the foreground amplitudedoes not improve over its prior, since in this example the cosmological signal is much largerthan the foreground where the sensitivity is highest. How this may change if we considerdifferent shapes and amplitudes for the cosmological component is beyond the scope of thiswork, and is left to future studies on the topic.The reconstruction of the same cosmological signal, now in the presence of the galacticforeground described above, is shown in fig. 12. In this case, the cosmological and theforeground SGWB components are the dominant contribution in different segments ofthe LISA frequency band. This allows both for a very accurate reconstruction of thecosmological component, and a precise measurement of the amplitude of the foregroundcomponent (well beyond its prior), as well as of the noise parameters. This is clear fromthe contour plots and from the triangle plot in figs. 12b and 12c. In this paper we have extended the work initiated in ref. [42] from a single channel to the useof the three TDI 1.5 channels for the LISA space-based interferometer. We have presenteda formalism and an algorithm to extract a SGWB signal from the LISA data, assuminga simplified model for the noise power spectral density and two illustrative foregrounds.We have worked in the AET channel combination, leading to vanishing off-diagonal termsin channel space when the LISA’s arms have equal length, and the noise power spectraldensities of all test masses are the same. Under these idealised conditions, we are ableto simplify the numerical computations by converting to the noise orthogonal AET TDI– 22 – Frequency [Hz]10 h G W Binned reconstruction (8 bins)
Data (used by the binner)LISA SciRDLISA PLS 4.0y, 0.75% eff, SNR=10Input signalReconstructed sensitivityReconstructed foregroundReconstructed signalBin extremesForeground 1 regionForeground 2 regionSignal 1 regionSignal 2 region (a) (bin 1) n s ( b i n ) (bin 2) n s ( b i n ) (bin 3) n s ( b i n ) (bin 4) n s ( b i n ) (bin 5) n s ( b i n ) (bin 6) n s ( b i n ) (bin 7) n s ( b i n ) (bin 8) n s ( b i n ) (b) A F G P P FG (c) Figure 11 : (a) Simultaneous reconstruction of a broken power-law signal (see section 4.2),an extra-galactic foreground component modelled as the power law in eq. (5.1), and thenoise spectrum. The initial number of bins was 23. Contour plots for the bin parametersand a triangle plot of the marginalized posterior for the two noise and one foregroundparameters are shown respectively in (b) and (c) .– 23 – Frequency [Hz]10 h G W Binned reconstruction (8 bins)
Data (used by the binner)LISA SciRDLISA PLS 4.0y, 0.75% eff, SNR=10Input signalReconstructed sensitivityReconstructed foregroundReconstructed signalBin extremesForeground 1 regionForeground 2 regionSignal 1 regionSignal 2 region (a) (bin 1) n s ( b i n ) (bin 2) n s ( b i n ) (bin 3) n s ( b i n ) (bin 4) n s ( b i n ) (bin 5) n s ( b i n ) (bin 6) n s ( b i n ) (bin 7) n s ( b i n ) (bin 8) n s ( b i n ) (b) A F G ) P P FG ) (c) Figure 12 : (a) Simultaneous reconstruction of a broken power-law signal (see section 4.2),an galactic foreground component with a model given in eq. (5.2), and the noise spectrum.The initial number of bins was 23. Contour plots for the bin parameters and a triangleplot of the marginalized posterior for the two noise and one foreground parameters areshown respectively in (b) and (c) . – 24 –ariables. In addition, we find that under these assumptions, the results are identical tousing the standard XYZ TDI variables.While in the code we have worked with the numerical response functions, we haveprovided the analytic form for the XX channel (as well as the YY and ZZ channels), andfor the AA (as well as the EE) and for the TT channels. We found that in the TT channelboth the signal response and the acceleration noise power spectrum are suppressed upto O (10 − ) Hz. As a consequence, TT can be used to impose a strong constraint on theinterferometric component of the noise, thus making possible to use the full frequency rangefor the reconstruction, which was not possible in our previous single-channel approach. Forthe same reason, only two channels contribute to the constraints on the SGWB signal, whichare thus improved by √ SGWBinner algorithm accordingly, we haveapplied it on three physically-motivated scenarios (i.e. a power-law, a broken power-lawand a bump GW signal), where we find it to be able to reconstruct the injected signalsaccurately. We check the per-bin optimisation approach of the
SGWBinner against an MCexploration of the full posterior and we found a good matching of the results. Finally, inorder to test the robustness of our reconstruction algorithm in presence of backgroundsand (astrophysical) foregrounds, we have performed some tests on signals given by thesum of an isotropic astrophysical component and a (unexpected) cosmological component.As astrophysical foregrounds we have considered an extra-galactic-like component anda galactic-like component, which are both foreseen in the LISA band. Besides being arobustness check, such a test is also relevant for testing disentangling ability for differentsignal contributions. We have seen that, when the injected signal is sufficiently large withrespect to the LISA sensitivity, the signal reconstruction is practically unaffected by thepresence of the foregrounds.In future work, we plan to apply our methodology to more realistic scenarios, both froman instrumental and an astrophysical point of view. The performance of our techniqueswill be tested against data sets that contain spurious and confusion signals as well, suchas the compact binaries in the vicinity of our galaxy [18]. Moreover, the case of varyingconstellation arm lengths, and unequal power spectra densities of the test-masses noiseswill be tested. This would further increase the complexity of the analysis by adding extradegrees of freedom and degrade the constraining capacity in our reconstruction. Neverthe-less we expect the simplified configuration studied here to be a reasonable approximationfor exploratory studies of SGWB.
Acknowledgments
We thank the LISA Publication and Presentation committee, in particular Nelson Chris-tensen for useful comments. R.F. was supported in part by the Department of Energy under– 25 –rant No. DE-SC0009919, and a grant from the Simons Foundation/SFARI 560536. Thework of M.P. was supported by Science and Technology Facilities Council consolidatedgrant ST/P000762/1. M.P. was supported in part by the National Science Foundationunder Grant No.NSF PHY-1748958. R.F. and M.P. would like to thank the Kavli Institutefor Theoretical Physics at UC Santa Barbara for the kind hospitality during part of thiswork. A.R. acknowledges funding from Italian Ministry of Education, University and Re-search (MIUR) through the “Dipartimenti di eccellenza” project Science of the Universe.N.K. was supported by the CNES DIA-PF post-doctoral fellowship program.
A Signal and detector technicalities
A.1 Polarization tensors
We work in natural units and in the Lorentz transverse-traceless gauge. We consider thecoordinate system { ˆ e x , ˆ e y , ˆ e z } arbitrarily oriented and at rest with respect to the isotropicSGWB we analyse. We neglect the motion of LISA in this frame; taking it into account isconceptually straightforward [46] but adds some technical complications that are unneces-sary for the purpose of the present paper. In this reference frame, the k wave-vector of anincoming plane GW sets the orthonormal basis [46, 59]ˆ u (ˆ k ) = ˆ k × ˆ e z | ˆ k × ˆ e z | , ˆ v (ˆ k ) = ˆ k × ˆ u , (A.1)(where ˆ k is unit vector in the direction of k and we will denote its magnitude by k = | k | )that we introduce to define the “plus” (+) and “cross” ( × ) polarization tensors: e (+) ab (ˆ k ) = ˆ u a ˆ u b − ˆ v a ˆ v b √ , e ( × ) ab (ˆ k ) = ˆ u a ˆ v b + ˆ v a ˆ u b √ . (A.2)Since ˆ u ( − ˆ k ) = − ˆ u (ˆ k ) and ˆ v ( − ˆ k ) = ˆ v (ˆ k ), the tensors e (+) ab and e ( × ) ab fulfill the conditions e + / × ab (ˆ k ) = e + / × ∗ ab (ˆ k ) , e + ab (ˆ k ) = e + ab ( − ˆ k ) , e × ab ( − ˆ k ) = − e × ab (ˆ k ) ,e + ab (ˆ k ) e + ab (ˆ k ) = 1 , e × ab (ˆ k ) e × ab (ˆ k ) = 1 , e + ab (ˆ k ) e × ab (ˆ k ) = 0 . (A.3)Analogously, we introduce the “right-handed” (R) and “left-handed” (L) polarization ten-sors: e Rab (ˆ k ) = ˆ u a + i ˆ v a √ u b + i ˆ v b √ , e Lab (ˆ k ) = ˆ u a − i ˆ v a √ u b − i ˆ v b √ . (A.4)They are related to the plus and cross polarization basis by means of the relationships e Rab (ˆ k ) = e + ab + i e × ab √ , e Lab (ˆ k ) = e + ab − i e × ab √ . (A.5)The superposition of all GWs reaching the position x at the time t can be expressedin terms of incoming plane waves and reads h ab ( x , t ) = (cid:90) + ∞−∞ d f (cid:90) Ω dΩ ˆ k e πif ( t − ˆ k · x ) (cid:88) A ˜ h A ( f, ˆ k ) e Aab (ˆ k ) , (A.6)– 26 –here f is the frequency of each plane wave, dΩ ˆ k is the infinitesimal solid angle from whichthe incoming wave with wave-vector k arrives, and finally ˜ h A ( f, ˆ k ) ≡ k ˜ h A ( k ). The index A can be used to label either the plus and cross polarizations or the the left- and right-handed polarizations, so that the sum can run indifferently over one or the other basis(i.e. A = + , × or A = L, R ). For an isotropic background the intergral over dΩ ˆ k only leadsto an overall factor, so that we can write h ab ( x , t ) = (cid:90) d k e − πi k · x (cid:88) A (cid:104) e πikt ˜ h A ( k ) e Aab (ˆ k ) + e − πikt ˜ h ∗ A ( − k ) e A ∗ ab ( − ˆ k ) (cid:105) . (A.7) A.2 TDI variables
Consider two test masses inside two of the three LISA spacecrafts. These test masses,labelled 1 and 2, are located at x and x and separated by the vector L ˆ l , with ˆ l =( x − x ) / | x − x | and L = 2 . × m. When a GW crosses the detector, a photon goingfrom a test mass to the other is received with the time shift [60, 61]∆ T ( t ) = ˆ l a ˆ l b (cid:90) L d s h ab ( t ( s ) , x ( s )) . (A.8)For the photon path, the leading order approximations t ( s ) = t − L + s and x ( s ) = x + s ˆ l suffice. It follows∆ T ( t ) = L (cid:90) d k e − πi k · x (cid:88) A (cid:104) e πikt M ( k , ˆ l ij ) ˜ h A ( k ) G A (ˆ k, ˆ l )+ e − πikt M ∗ ( − k , ˆ l ij ) ˜ h ∗ A ( − k ) G A ∗ ( − ˆ k, ˆ l ) (cid:105) , (A.9)with G Ai (ˆ k, ˆ l ij ) ≡ ˆ l aij ˆ l bij e Aab (ˆ k ) , M ( k , ˆ l ij ) ≡ e πikL (1 − ˆ k · ˆ l ij ) sinc (cid:104) kL (1 − ˆ k · ˆ l ij ) (cid:105) , (A.10)and sinc( x ) ≡ sin( πx ) / ( πx ).In practice, LISA does not observe this time shift, but instead observes the corre-sponding Doppler effect. Let us assume that the test mass 2 emits photons with frequency ν . The test mass 1 receives such photons with the fractional Doppler frequency shift∆ F ( t ) ≡ ∆ ν ( t ) /ν = − d∆ T ( t ) / d t .In LISA neither the frequency ν nor the positions of the test masses are not perfectlyknown. For this reason, the observation of a GW requires to perform time delay interferom-etry (TDI), that is, to compare the fractional frequency shifts of the same light pulse splitinto different paths. Moreover, to take into account the uncertainties on the positions ofthe test masses, the relative frequency shifts are measured after conveying the laser pulsesalong paths that are quite complex. Without reaching the level of complexity of the pathsadopted in LISA, which are still under discussion, hereafter we focus on paths that containall the conceptual elements of the realistic ones.A key element in common to every TDI measurement is the fractional Doppler shift ofa photon following paths such that photon frequency noise and motion of optical benches– 27 –ancels. For concreteness let us consider the case of the closed path 1 → →
1, with aphoton starting from x at time t − L and returning back to x at time t . In this case thefinal fractional Doppler shift is∆ F ( t ) ≡ ∆ F ( t − L ) + ∆ F ( t )= (cid:90) d k e − πi k · x ikf (cid:63) (cid:88) A (cid:104) e πik ( t − L ) T ( k , ˆ l ) ˜ h A ( k ) G A (ˆ k, ˆ l ) − e − πik ( t − L ) T ( k , ˆ l ) ˜ h ∗ A ( − k ) G A ∗ ( − ˆ k, ˆ l ) (cid:105) , (A.11)where f (cid:63) = (2 πL ) − is the detector characteristic frequency, and T is the so called transferfunction defined as T ( k , ˆ l ij ) ≡ e πikL (1 − ˆ k · ˆ l ij ) sinc (cid:104) kL (1 + ˆ k · ˆ l ij ) (cid:105) + e − πikL (1+ˆ k · ˆ l ij ) sinc (cid:104) kL (1 − ˆ k · ˆ l ij ) (cid:105) . (A.12)We can now compare the fractional Doppler shifts of two laser beams of identical frequency,with the first one following the path 1 → → → →
1. Their relative Doppler shift measured at 1 at the time t , here denoted with ∆ F ( t ),results ∆ F ( t ) = ∆ F ( t ) − ∆ F ( t ) . (A.13)Closer inspection reveals that this variable is only a good TDI variable for an equilat-eral configuration. For configurations that depart from the equilateral configuration morecomplicated TDI variables are required. For a configurations that depart from an equi-lateral configuration, but for which the time variation of the armlengths is slow, a path1 → → → → → → → → F . ( t ) = ∆ F ( t − L ) + ∆ F ( t ) (TDI 1.5 variable) . (A.14)We will work with these TDI 1.5 variables and introduce the shorthand notation X =∆ F . . The TDI variables Y and Z are obtained by cyclic permutation. However, ourtechniques can easily be adapted to more complex TDI variables that also account for thetime evolution of the armlengths. A.3 Response functions
By substituting eq. (A.11) into eq. (A.13), one can see that the signal contribution to∆ F ( t ) is∆ F ( t ) = (cid:90) d k e − πi k · x (2 πikL ) (cid:88) A (cid:104) e πik ( t − L ) R A ( k , ˆ l , ˆ l )˜ h A ( k )+ − e − πik ( t − L ) R A ∗ ( − k , ˆ l , ˆ l )˜ h ∗ A ( − k ) (cid:105) , (A.15)with R Ai ( k , ˆ l ij , ˆ l ik ) ≡ G A (ˆ k, ˆ l ij ) T ( k , ˆ l ij ) − G A (ˆ k, ˆ l ik ) T ( k , ˆ l ik ) . (A.16)– 28 –he functions R Ai are the so called LISA response function, and describe how the detectorresponds to a plane wave of wave-vector k when the LISA test masses i and j are orientedalong the direction ˆ l ij . It then immediately follows that∆ F . ( t ) = (cid:90) d k e − πi k · x (2 πikL ) (cid:88) A (cid:104) e πik ( t − L ) W ( kL ) R A ( k , ˆ l , ˆ l )˜ h A ( k )+ − e − πik ( t − L ) W ∗ ( kL ) R A ∗ ( − k , ˆ l , ˆ l )˜ h ∗ A ( − k ) (cid:105) , with W ( kL ) ≡ e − πikL −
1. Notice that | W ( kL ) | = 2 [1 − cos(4 πkL )] = 4 sin (2 πkL ). Thetwo-point correlation functions of ∆ F i ( jk ) ( t ) can then be expressed as (cid:10) ∆ F i ( jk ) ( t )∆ F l ( mn ) ( t ) (cid:11) = (cid:90) d k e − πi k · ( x i − x l ) πkL ) | W ( kL ) | ×× (cid:88) A,A (cid:48) R Ai ( k , ˆ l ij , ˆ l ik ) R A (cid:48) l ∗ ( k , ˆ l lm , ˆ l ln ) P AA (cid:48) h ( k )4 πk , (A.17)where we have used (cid:104) ˜ h A ( k )˜ h ∗ A ( − k ) (cid:105) = δ ( k + k ) P A A h ( k )4 πk , (cid:104) ˜ h A ( k )˜ h A ( k ) (cid:105) = 0 , (A.18)and we have taken P A A h ( k ) to be real. Moreover, since we restrict ourselves to a back-ground with vanishing Stokes parameters Q and U [62], P A A h ( k ) is diagonal in the L/R basis: P LLh ( k ) = I − V , P
RRh ( k ) = I + V . (A.19)Therefore, V must be vanishing too for the assumed, non-chiral background ( i.e. P LLh ( k ) = P RRh ( k )).Finally eq. (A.17) can be written as (cid:10) ∆ F i ( jk ) ( t )∆ F l ( mn ) ( t ) (cid:11) = (cid:88) A (cid:90) d k πkL ) | W ( kL ) | δ AA (cid:48) ˜ R AA (cid:48) il ( jk )( mn ) ( k ) P Ah ( k ) , (A.20)where ˜ R AA (cid:48) il ( jk )( mn ) ( k ) ≡ π (cid:90) d ˆ k e − πi k · ( x i − x l ) R Ai ( k , ˆ l ij , ˆ l ik ) R A (cid:48) ∗ l ( k , ˆ l lm , ˆ l ln ) . (A.21)Since for LISA we have ˜ R LLil ( jk )( mn ) = ˜ R RRil ( jk )( mn ) ≡ ˜ R il ( jk )( mn ) , for V = 0 we obtain (cid:10) ∆ F i ( jk ) ( t )∆ F l ( mn ) ( t ) (cid:11) = (cid:90) d k πkL ) | W ( kL ) | ˜ R il ( jk )( mn ) ( k ) P h ( k ) , (A.22)with P h ( k ) ≡ P Lh ( k ) = P Rh ( k ). At this point it is convenient to define the integrand (up tothe factor P h ( k )) as R ij ( k ) ≡ πkL ) | W ( kL ) | ˜ R il ( jk )( mn ) ( k ) . (A.23)– 29 – .4 Noise spectra For completeness, we also include a review of the calculation of the noise spectra used inthe main text. We assume that the laser frequency noise has been removed through time-delay interferometry. The dominant residual sources of noise are then associated with theinterferometer measurement system (IMS), in particular shot noise for the light travelingbetween spacecrafts, and acceleration noise, fluctuations in the proof mass velocities, causedby residual electrostatic and magnetic forces on the proof masses.For concreteness, we assume that the laser beam is reflected off the proof mass uponarrival for transmission between optical benches on different spacecrafts, and before de-parture for transmission between optical benches on the same spacecraft. In this case, theacceleration noise contribution to ˜ d X is [63, 64]˜ d accX = (cid:16) e πifL (cid:17) (cid:16) e πifL − (cid:17) ˆ l v L + 2 e πifL (cid:16) − e πifL (cid:17) ˆ l v L − (cid:16) e πifL (cid:17) (cid:16) e πifL − (cid:17) ˆ l v R − e πifL (cid:16) − e πifL (cid:17) ˆ l v R , (A.24)where L ij is the distance between spacecrafts i and j , which is assumed to be time-independent, ˆ l ij is a unit vector pointing from spacecraft j to i as before, and v iL and v iR are the fluctuations in the velocities of the proof mass in the left and right opticalbench on spacecraft i , respectively. The expressions for the variables ˜ d accY and ˜ d accZ can beobtained by cyclic permutation.If we assume that the fluctuations of the different probe masses are uncorrelated (evenfor masses on the same spacecraft) and are all characterized by the same power spectrum (cid:104) ˆ n ˜ v L ( f )ˆ n ˜ v ∗ L ( f (cid:48) ) (cid:105) (cid:48) = (cid:104) ˆ n ˜ v R ( f )ˆ n ˜ v ∗ R ( f (cid:48) ) (cid:105) (cid:48) = · · · = P acc ( f ) , (A.25)the acceleration noise contribution to the power spectra are (cid:104) ˜ d X ( f ) ˜ d ∗ X ( f (cid:48) ) (cid:105) (cid:48) = 8 (cid:0) (2 πf L ) + 2 sin (2 πf L )+ sin (2 πf ( L − L )) + sin (2 πf ( L + L )) (cid:1) P acc ( f ) , (A.26) (cid:104) ˜ d X ( f ) ˜ d ∗ Y ( f (cid:48) ) (cid:105) (cid:48) = − e πif ( L − L ) sin (2 πf L ) sin (2 πf L ) cos (2 πf L ) P acc ( f ) . (A.27)We can proceed in the same way for IMS noise, and will assume it is dominated by shotnoise. In this case, we can neglect the contribution from transmission between the opticalbenches on a single spacecraft, and will collect all sources of IMS noise for transmissionfrom spacecraft j to i into a single term ˜ n ij . The IMS noise contribution to ˜ d X then is˜ d IMSX = (cid:16) e πifL − (cid:17) ˜ n + e πifL (1 − e πifL )˜ n − (cid:16) e πifL − (cid:17) ˜ n − e πifL (1 − e πifL )˜ n . (A.28)The corresponding expressions for the variables ˜ d accY and ˜ d accZ can again be found by cyclicpermutation.If we assume that the IMS contributions for the different Doppler variables uncorrelatedso that the noise covariance matrix is diagonal, and that the auto-spectra have identical– 30 –tatistical properties so that (cid:104) ˜ n ( f )˜ n ( f (cid:48) ) (cid:105) (cid:48) = (cid:104) ˜ n ( f )˜ n ( f (cid:48) ) (cid:105) (cid:48) = · · · = P IMS ( f ) , (A.29)the IMS noise contribution to the power spectra are (cid:104) ˜ d X ( f ) ˜ d ∗ X ( f (cid:48) ) (cid:105) (cid:48) = 8 (cid:0) sin (2 πf L ) + sin (2 πf L ) (cid:1) P IMS ( f ) , (A.30) (cid:104) ˜ d X ( f ) ˜ d ∗ Y ( f (cid:48) ) (cid:105) (cid:48) = − e πif ( L − L ) sin (2 πf L ) sin (2 πf L ) cos (2 πf L ) P IMS ( f ) . (A.31)So far we have assumed that L ij are time-independent but not necessarily equal. Thenoise spectra in the main text assume an equilateral configuration, and can be obtained bysetting L = L = L = L . Collecting both acceleration and IMS noise, we find (cid:104) ˜ d X ( f ) ˜ d ∗ X ( f (cid:48) ) (cid:105) (cid:48) = 16 sin (2 πf L ) { (3 + cos (4 πf L )) P acc ( f ) + P IMS ( f ) } , (A.32) (cid:104) ˜ d X ( f ) ˜ d ∗ Y ( f (cid:48) ) (cid:105) (cid:48) = − (2 πf L ) cos (2 πf L ) [4 P acc ( f ) + P IMS ( f )] . (A.33)We see that for the equilateral configuration the noise covariance matrix for ˜ d X , ˜ d Y , and˜ d Z becomes real. We also see that we can approximate the configuration as equilateral if2 πf ∆ L/c (cid:28) References [1]
LISA collaboration,
Laser Interferometer Space Antenna , .[2] LIGO Scientific collaboration,
Advanced LIGO , Class. Quant. Grav. (2015) 074001[ ].[3] F. A. et al, Advanced Virgo: a second-generation interferometric gravitational wave detector , Class. Quant. Grav. (2015) 024001.[4] S. Hild et al., Sensitivity Studies for Third-Generation Gravitational Wave Observatories , Class. Quant. Grav. (2011) 094013 [ ].[5] KAGRA collaboration,
Detector configuration of KAGRA: The Japanese cryogenicgravitational-wave detector , Class. Quant. Grav. (2012) 124007 [ ].[6] KAGRA collaboration,
Interferometer design of the KAGRA gravitational wave detector , Phys. Rev. D (2013) 043007 [ ].[7] LIGO Scientific collaboration,
Exploring the Sensitivity of Next Generation GravitationalWave Detectors , Class. Quant. Grav. (2017) 044001 [ ].[8] A. Sesana, F. Haardt, P. Madau and M. Volonteri, Low - frequency gravitational radiationfrom coalescing massive black hole binaries in hierarchical cosmologies , Astrophys. J. (2004) 623 [ astro-ph/0401543 ].[9] A. Sesana, F. Haardt, P. Madau and M. Volonteri,
The gravitational wave signal frommassive black hole binaries and its contribution to the LISA data stream , Astrophys. J. (2005) 23 [ astro-ph/0409255 ].[10] E. Barausse,
The evolution of massive black holes and their spins in their galactic hosts , Mon. Not. Roy. Astron. Soc. (2012) 2533 [ ]. – 31 –
11] A. Klein et al.,
Science with the space-based interferometer eLISA: Supermassive black holebinaries , Phys. Rev. D (2016) 024003 [ ].[12] M. Bonetti, A. Sesana, F. Haardt, E. Barausse and M. Colpi, Post-Newtonian evolution ofmassive black hole triplets in galactic nuclei – IV. Implications for LISA , Mon. Not. Roy.Astron. Soc. (2019) 4044 [ ].[13] S. Babak, J. Gair, A. Sesana, E. Barausse, C. F. Sopuerta, C. P. Berry et al.,
Science withthe space-based interferometer LISA. V: Extreme mass-ratio inspirals , Phys. Rev. D (2017) 103012 [ ].[14] G. Nelemans, L. Yungelson and S. Portegies Zwart, Short- period AM CVn systems asoptical, x-ray and gravitational wave sources , Mon. Not. Roy. Astron. Soc. (2004) 181[ astro-ph/0312193 ].[15] A. J. Ruiter, K. Belczynski, M. Benacquista, S. L. Larson and G. Williams,
The LISAGravitational Wave Foreground: A Study of Double White Dwarfs , Astrophys. J. (2010)1006 [ ].[16] T. Marsh,
Double white dwarfs and LISA , Class. Quant. Grav. (2011) 094019[ ].[17] V. Korol, E. M. Rossi, P. J. Groot, G. Nelemans, S. Toonen and A. G. Brown, Prospects fordetection of detached double white dwarf binaries with Gaia, LSST and LISA , Mon. Not.Roy. Astron. Soc. (2017) 1894 [ ].[18] V. Korol, E. M. Rossi and E. Barausse,
A multimessenger study of the Milky Way’s stellardisc and bulge with LISA, Gaia, and LSST , Mon. Not. Roy. Astron. Soc. (2019) 5518[ ].[19] P. Bender and D. Hils,
Confusion noise level due to galactic and extragalactic binaries , Class.Quant. Grav. (1997) 1439.[20] C. R. Evans, I. Iben and L. Smarr, Degenerate dwarf binaries as promising, detectablesources of gravitational radiation , Astrophys. J. (1987) 129.[21]
LIGO Scientific, Virgo collaboration,
Search for the isotropic stochastic backgroundusing data from Advanced LIGOs second observing run , Phys. Rev.
D100 (2019) 061101[ ].[22] M. Bonetti and A. Sesana,
Gravitational wave background from extreme mass ratio inspirals , .[23] L. V. Dall’Armi, A. Ricciardone, N. Bartolo, D. Bertacca and S. Matarrese, The Imprint ofRelativistic Particles on the Anisotropies of the Stochastic Gravitational-Wave Background , .[24] C. Caprini and D. G. Figueroa, Cosmological Backgrounds of Gravitational Waves , Class.Quant. Grav. (2018) 163001 [ ].[25] E. Barausse et al., Prospects for Fundamental Physics with LISA , .[26] E. Barausse, R. Brito, V. Cardoso, I. Dvorkin and P. Pani, The stochastic gravitational-wavebackground in the absence of horizons , Class. Quant. Grav. (2018) 20LT01 [ ].[27] S. Alexander, E. McDonough, R. Sims and N. Yunes, Hidden-Sector Modifications toGravitational Waves From Binary Inspirals , Class. Quant. Grav. (2018) 235012[ ]. – 32 –
28] C. Caprini et al.,
Science with the space-based interferometer eLISA. II: Gravitational wavesfrom cosmological phase transitions , JCAP (2016) 001 [ ].[29] C. Caprini et al., Detecting gravitational waves from cosmological phase transitions withLISA: an update , JCAP (2020) 024 [ ].[30] P. Auclair et al., Probing the gravitational wave background from cosmic strings with LISA , .[31] N. Bartolo et al., Science with the space-based interferometer LISA. IV: Probing inflationwith gravitational waves , JCAP (2016) 026 [ ].[32] J. Garcia-Bellido, M. Peloso and C. Unal, Gravitational waves at interferometer scales andprimordial black holes in axion inflation , JCAP (2016) 031 [ ].[33] M. R. Adams and N. J. Cornish, Detecting a Stochastic Gravitational Wave Background inthe presence of a Galactic Foreground and Instrument Noise , Phys. Rev.
D89 (2014) 022001[ ].[34] M. Pieroni and E. Barausse,
Foreground cleaning and template-free stochastic backgroundextraction for LISA , .[35] T. Robson and N. J. Cornish, Detecting gravitational wave bursts with lisa in the presence ofinstrumental glitches , Phys. Rev. D (2019) 024019.[36] M. Armano, H. Audley, J. Baird, P. Binetruy, M. Born, D. Bortoluzzi et al., Beyond therequired lisa free-fall performance: New lisa pathfinder results down to µ Hz,
Phys. Rev.Lett. (2018) 061101.[37] M. Armano, H. Audley, G. Auger, J. T. Baird, M. Bassan, P. Binetruy et al.,
Sub-femto- g free fall for space-based gravitational wave observatories: Lisa pathfinder results , Phys. Rev.Lett. (2016) 231101.[38] F. Gibert, M. Nofrarias, N. Karnesis, L. Gesa, V. Mart´ın, I. Mateos et al.,
Thermo-elasticinduced phase noise in the LISA pathfinder spacecraft , Classical and Quantum Gravity (2015) 045014.[39] M. Armano, H. Audley, J. Baird, P. Binetruy, M. Born, D. Bortoluzzi et al., Temperaturestability in the sub-milliHertz band with LISA Pathfinder , Monthly Notices of the RoyalAstronomical Society (2019) 3368[ https://academic.oup.com/mnras/article-pdf/486/3/3368/28536406/stz1017.pdf ].[40] Y. B. Ginat, V. Desjacques, R. Reischke and H. B. Perets,
The Probability Distribution ofAstrophysical Gravitational-Wave Background Fluctuations , .[41] N. Karnesis, A. Petiteau and M. Lilley, A template-free approach for detecting a gravitationalwave stochastic background with LISA , .[42] C. Caprini, D. G. Figueroa, R. Flauger, G. Nardini, M. Peloso, M. Pieroni et al., Reconstructing the spectral shape of a stochastic gravitational wave background with LISA , JCAP (2019) 017 [ ].[43] C. J. Hogan and P. L. Bender, Estimating stochastic gravitational wave backgrounds withSagnac calibration , Phys. Rev. D (2001) 062002 [ astro-ph/0104266 ].[44] M. R. Adams and N. J. Cornish, Discriminating between a Stochastic Gravitational WaveBackground and Instrument Noise , Phys. Rev. D (2010) 022002 [ ]. – 33 –
45] A. Sesana,
Prospects for Multiband Gravitational-Wave Astronomy after GW150914 , Phys.Rev. Lett. (2016) 231102 [ ].[46] V. Domcke, J. Garcia-Bellido, M. Peloso, M. Pieroni, A. Ricciardone, L. Sorbo et al.,
Measuring the net circular polarization of the stochastic gravitational wave background withinterferometers , .[47] S. Babak and A. Petiteau, LISA Data Challenge Manual , Tech. Rep.LISA-LCST-SGS-MAN-002, APC Paris, 7, 2020.[48] T. Robson, N. J. Cornish and C. Liu,
The construction and use of LISA sensitivity curves , Class. Quant. Grav. (2019) 105011 [ ].[49] T. A. Prince, M. Tinto, S. L. Larson and J. W. Armstrong, Lisa optimal sensitivity , Phys.Rev. D (2002) 122002.[50] J. Bond, A. H. Jaffe and L. Knox, Radical compression of cosmic microwave backgrounddata , Astrophys. J. (2000) 19 [ astro-ph/9808264 ].[51] J. Sievers et al.,
Cosmological parameters from Cosmic Background Imager observations andcomparisons with BOOMERANG, DASI, and MAXIMA , Astrophys. J. (2003) 599[ astro-ph/0205387 ].[52]
WMAP collaboration,
First year Wilkinson Microwave Anisotropy Probe (WMAP)observations: Parameter estimation methodology , Astrophys. J. Suppl. (2003) 195[ astro-ph/0302218 ].[53] S. Hamimeche and A. Lewis,
Likelihood Analysis of CMB Temperature and PolarizationPower Spectra , Phys. Rev. D (2008) 103013 [ ].[54] H. Akaike, A new look at the statistical model identification , IEEE Transactions onAutomatic Control
19 (6) (1974) 716.[55] J. Torrado and A. Lewis,
Cobaya: Code for Bayesian Analysis of hierarchical physicalmodels , .[56] W. J. Handley, M. P. Hobson and A. N. Lasenby, polychord: nested sampling for cosmology. , Mon.Not.Roy.Astron.Soc. (2015) L61 [ ].[57] W. J. Handley, M. P. Hobson and A. N. Lasenby,
POLYCHORD: next-generation nestedsampling , Mon.Not.Roy.Astron.Soc. (2015) 4384 [ ].[58] A. Lewis,
Getdist: a python package for analysing monte carlo samples , .[59] N. Bartolo, V. Domcke, D. G. Figueroa, J. Garca-Bellido, M. Peloso, M. Pieroni et al., Probing non-Gaussian Stochastic Gravitational Wave Backgrounds with LISA , JCAP (2018) 034 [ ].[60] F. B. Estabrook and H. D. Wahlquist,
Response of Doppler spacecraft tracking togravitational radiation , Gen. Relativ. Gravit. (1975) 439.[61] J. D. Romano and N. J. Cornish, Detection methods for stochastic gravitational-wavebackgrounds: a unified treatment , Living Rev. Rel. (2017) 2 [ ].[62] M. Maggiore, Gravitational Waves. Vol. 2: Astrophysics and Cosmology . Oxford UniversityPress, 3, 2018.[63] F. B. Estabrook, M. Tinto and J. W. Armstrong,
Time delay analysis of LISA gravitationalwave data: Elimination of spacecraft motion effects , Phys. Rev. D (2000) 042002. – 34 –
64] M. Tinto and S. V. Dhurandhar,
Time-Delay Interferometry , Living Rev. Rel. (2014) 6.(2014) 6.