Measuring the primordial gravitational wave background in the presence of other stochastic signals
PPrepared for submission to JCAP
Measuring the primordialgravitational wave background inthe presence of other stochasticsignals
D. Poletti
International School for Advanced Studies (SISSA),Via Bonomea 265, 34136, Trieste, ItalyInstitute for Fundamental Physics of the Universe (IFPU),Via Beirut 2, 34014 Trieste, ItalyNational Institute for Nuclear Physics (INFN) Sezione di Trieste,Padriciano, 99, 34149 Trieste, ItalyE-mail: [email protected]
Abstract.
Standard methodologies for the extraction of the stochastic gravitational wavebackground (SGWB) from auto- or cross-correlation of interferometric signals often involvethe use of a filter function. The standard optimal filter maximizes the signal-to-noise ratio(SNR) between the total SGWB and the noise. We derive expressions for the optimal filter andSNR in the presence of a target SGWB plus other unwanted components. We also generalizethe methodology to the case of template-free reconstruction. The formalism allows to easilyperform analyses and forecasts that marginalize over foreground signals, such as the typical Ω GW ∝ f / background arising from binary coalescence. We demonstrate the methodologywith the LISA mission and discuss possible extensions and domains of application. a r X i v : . [ g r- q c ] J a n ontents Measuring the stochastic gravitational wave background (SGWB) is among the main objec-tives of existing and planned gravitational wave observatories. It is expected to be producedby a wide variety of astrophysical or cosmological phenomena, see [1–4] for reviews. An astro-physical SGWB certainly exists. In the recent years, we have witnessed a series of detectionsof gravitational waves transient signals produced by binary back holes and binary neutronstars mergers [5–11]. The signal produced by these systems, integrated over the whole historyof the universe, must constitute a background of gravitational waves. Other contributions areexpected to come from core-collapse supernovae, non-axisymmetric spinning neutron starsand magnetars (see, e.g., [12–15]). The background of gravitational waves produced throughthese mechanisms offers a window on a large number of astrophysical processes over the entirehistory of the universe.The SGWB is also considered a major probe of inflation [16–19]. By postulating a periodof accelerated expansion before the standard Big-Bang universe, inflation is capable of solvingthe horizon, flatness and monopole problems. The standard inflationary mechanism typicallyinvolves a particle beyond the Standard Model dominating the energy content of the infantuniverse. Its quantum fluctuations—blown up by the quasi-De Sitter expansion—naturallyact as seeds of the cosmological structures we observe in the low redshift universe as wellas the acoustic waves observed in the cosmic microwave background (CMB) at z (cid:39) .– 1 –uantum fluctuation of the gravitational field itself are expected to undergo a similar fateand additional GW contributions may arise from the presence of fields in addition or in placeof the simple scalar-field inflaton. The resulting SGWB should be still present today at thefrequencies of existing and planned GW direct detection experiments. Likewise these GWwould leave an imprint in the CMB, best observable in the sought-after B-mode polarization,key target of the ultimate (or near ultimate) CMB observatories [20, 21].Topological strings [22, 23] and superstrings [24, 25] attracted a lot of attention in thiscontext. They also generate a SGWB resulting from the superposition of powerful bursts ofGWs produced by cusps and kinks propagating on string loops. Again, this type of stochasticbackground can be constrained with both CMB [26] and direct detection experiments [27].Either cosmological or astrophysical, no detection of SGWB was obtained so far. Never-theless, the measurements of individual gravitational wave events [5–11] allowed to refine thepredicted amplitude of the SGWB from compact binary coalescences [28, 29] at 25 GHz—the frequency where the aLIGO and aVIRGO observatories are most sensitive to the SGWB.Moreover, at the other end of the frequency spectrum—in the 1-100 nHz range—the common-spectrum stochastic signal detected by the NANOGrav Collaboration across an array of 45pulsars [30] may be sourced by a SGWB, even though we have to wait for conclusive evi-dence of the quadrupolar correlations [31] before claiming a detection of SGWB and makinginferences on its origin.It is clear that some form of SGWB exists and, given the wide variety of possible origins,the analyses (as well as forecasts) will have to account for the possibility that multiple back-ground sources are contributing at the same time. The problem has been already explored[32–34]. In [34], in particular, it was recognised that under some conditions the SGWB ex-traction and foreground rejection reduces to a linear problem. This paper effectively pushesfurther this consideration. We take a different perspective that we believe will make the prob-lem (and solution) more natural and clear to the readers familiar to the popular approach tothe foreground-free problem by [35, 36]. In particular, we describe their procedure in Section3, after briefly summarizing the main quantities involved in the studies of SGWBs (Section2). It is based on the correlation of interferometry signals employing a frequency-dependentfilter, which determines the optimality of the correlation. The optimal filter is known to beessentially equal to the ratio between then strain of the SGWB and the strain-equivalent noisevariance of the two signals being correlated (Section 3.1). However, if we are interested inonly one of the components of the SGWB, the filter can be amended to optimally marginalizeover the unwanted components—in addition to optimally coadding the contributions of thedifferent frequencies, as done by the standard filter. We propose this new filter in Section 3.2(see Appendix A for the demonstration) and apply it on the LISA example in Section 5. InSection 4 and Appendix C we generalize the approach to allow the reconstruction of SGWBsthat are not known a priori. We summarize here the key definitions about isotropic and unpolarized SGWB, to whichwe restrict ourselves in the present paper. The background of gravitational waves can beexpressed as a plane wave expansion of the metric perturbation in the transverse-tracelessgauge as h ab ( t, (cid:126)x ) = (cid:88) P (cid:90) ∞−∞ df (cid:90) S d ˆ n h p ( f, ˆ n ) e i πf ( t − ˆ n · (cid:126)x/c ) e Pab (ˆ n ) , (2.1)– 2 –here the wave vector has been expressed in terms of its normalization and the unit vectorspecifying the direction of propagation, (cid:126)k ≡ πf ˆ n/c . We project the polarization state of thegravitational wave on the “plus” and “cross” bases, P = + , × , e + ab (ˆ n ) = ˆ m a ˆ m b − ˆ n a ˆ n b , (2.2) e × ab (ˆ n ) = ˆ m a ˆ m b + ˆ n a ˆ n b , (2.3)where ˆ m and ˆ n are unit vectors orthogonal to ˆ n and to each other. Note that the basisis normalized such that e Pab e P (cid:48) ,ab = 2 δ P P (cid:48) . In this article we assume that the observedgravitational waves belong to a SGWB that is Gaussian distributed, with zero mean andvariance (cid:104) h ∗ P ( f, ˆ n ) h P (cid:48) ( f (cid:48) , ˆ n (cid:48) ) (cid:105) = 18 π δ (ˆ n, ˆ n (cid:48) ) δ P P (cid:48) δ ( f − f (cid:48) ) S h ( f ) , (2.4)where the delta functions reflect the assumption that the background is isotropic, unpolarizedand stationary. The spectral density function S h ( f ) is a real function satisfying S h ( − f ) = S h ( f ) . Eq. (2.4) also sets our convention for its normalization, which varies across theliterature and in our case satisfies (cid:88) P P (cid:48) (cid:90) d ˆ n d ˆ n (cid:48) (cid:104) h ∗ P ( f, ˆ n ) h P (cid:48) ( f (cid:48) , ˆ n (cid:48) ) (cid:105) = δ ( f − f (cid:48) ) S h ( f ) . (2.5)The information carried by S h is often expressed in terms of Ω gw ( f ) ≡ ρ cr dρ gw ( f ) d ln f , (2.6)where ρ cr = 3 c H / πG it the critical density of the universe today and ρ gw = c πG (cid:104) ˙ h ab ( t, (cid:126)x ) ˙ h ab ( t, (cid:126)x ) (cid:105) (2.7)is the energy density in gravitational waves. Given our convention, Ω gw and S h are relatedby Ω gw ( f ) = 4 π H f S h ( f ) (2.8) The output of a gravitational wave detector I is a time stream d I ( t ) = s I ( t ) + n I ( t ) . (3.1)The signal s is a real function of time and is sourced by the gravitational waves hitting thedetector s I ( t ) = (cid:88) P =+ , × (cid:90) ∞−∞ df (cid:90) d ˆ n F PI (ˆ n, f ) h p ( f, ˆ n ) e i πf ( t − ˆ n · (cid:126)x I /c ) , (3.2)where (cid:126)x is the location of the detector and the response function F PI (ˆ n, f ) depends on prop-erties of the detector such as the geometry and orientation—which we assume constant in– 3 –ime. The noise term n I is also a real quantity, assumed to be stationary and Gaussian, withvariance given by (cid:104) ˜ n I ( f )˜ n ∗ J ( f (cid:48) ) (cid:105) = 12 δ ( f − f (cid:48) ) N IJ ( f ) . (3.3)The / prefactor is conventional but, unlike the one of S h , this normalization is consistentacross the literature. Note that in Eq. (3.3) we have considered the possibility that I and J are different detectors with correlated noise, as in the case of the XY Z channels of the LISAmission [37].The searches of stochastic gravitational wave backgrounds typically involve the correla-tion of (possibly different) detectors, I and J , both collecting signal over a time Tx IJ ≡ (cid:90) T/ − T/ dt (cid:90) T/ − T/ dt (cid:48) (cid:18) d I ( t ) d J ( t (cid:48) ) − N IJ ( | t − t (cid:48) | ) (cid:19) Q IJ ( t, t (cid:48) ) . (3.4) N IJ ( | t − t (cid:48) | ) / (cid:104) n I ( t ) n J ( t (cid:48) ) (cid:105) = (cid:82) ∞−∞ df e i πf ( t − t (cid:48) ) N IJ ( f ) / sources the noise bias. The filterfunction Q can be arbitrarily chosen in order to maximize the signal-to-noise of the crosscorrelation. In Section 3.1 and 3.2, we discuss in detail the optimization of Q . For timebeing, we only note that because of the stationarity of both signal and noise, Q must be afunction of | t − t (cid:48) | . We assume now that Q ( | t − t (cid:48) | ) is non-negligible only for time differencesmuch smaller than the observation time. This allows to push to infinity the extremes of oneof the integrals in Eq. (3.4), and thus to get x IJ = (cid:90) ∞−∞ df (cid:90) ∞−∞ df (cid:48) δ T ( f − f (cid:48) ) ˜ d ∗ I ( f ) ˜ d J ( f (cid:48) ) ˜ Q IJ ( f (cid:48) ) − T (cid:90) ∞−∞ df N IJ ( f ) ˜ Q IJ ( f ) . (3.5)The function δ T ( f ) = sin( πf T ) /πf has the following properties. It converges to the Dirac δ function as T → ∞ . In the same limit, δ T converges to T δ . For finite T , δ T (0) = T .The expected value and variance of x can be computed from those of ˜ d ∗ I ( f ) ˜ d J ( f (cid:48) ) , whichin turn have simple expressions thanks to Eqs. (2.5), (3.2) and (3.3)—assuming that all thenon-stationary, anisotropic or polarized signals have been removed or are negligible—, (cid:104) x IJ (cid:105) = T (cid:90) ∞ df S h ( f ) R IJ ( f ) ˜ Q IJ ( f ) , (3.6)where we have defined the response function R IJ ( f ) ≡ π (cid:90) S d ˆ n (cid:88) P F ∗ PI (ˆ n, f ) F PJ (ˆ n, f ) e − i πf ˆ n · ( (cid:126)x J − (cid:126)x i ) /c , (3.7)which depends on the properties of the detectors and their relative distance and orientation.The variance receives a contribution from both the signal and noise, but the noise poweris typically much higher than the signal. Therefore, we can ignore s in the time streams ofthe detectors and get (cid:104) x IJ (cid:105) − (cid:104) x IJ (cid:105) = T (cid:90) ∞ df (cid:0) N II ( f ) N JJ ( f ) + N IJ ( f ) (cid:1) | ˜ Q ( f ) | , (3.8)The following two sections discuss how to optimize the signal-to-noise ratio. SNR = (cid:104) x IJ (cid:105) (cid:113) (cid:104) x IJ (cid:105) − (cid:104) x IJ (cid:105) (3.9) The tilde denotes the Fourier transform, defined as ˜ g ( f ) ≡ (cid:82) ∞−∞ dt e − πfti g ( t ) – 4 – .1 Standard optimal filtering In this section we derive the standard expression of the optimal Q [35, 36], thus summarizingthe approach commonly adopted in the literature and paving the road for the new filter thatwe propose in the next section.We want to find the Q ( f ) that maximises the SNR. In order to do it, it is very convenientto express Eqs. (3.6) and (3.8) in term of the following inner product of complex functions a · b = T (cid:90) ∞ df a ∗ ( f ) b ( f ) N ( f ) (3.10)where we have defined for convenience N ( f ) ≡ N II ( f ) N JJ ( f ) + N IJ ( f ) , (3.11)In contrast with the entire existing literature, we include T in the definition of the innerproduct. We denote with boldface the objects (lowercase for vectors and upper case formatrices) that obey to this inner product. For convenience we define the following vectors h ≡ S h ( f ) R IJ ( f ) N ( f ) (3.12) q ≡ ˜ Q ( f ) (3.13)The expression for the SNR becomes SNR = ( q · h ) q · q , (3.14)which is obviously independent of the normalization of q and is maximized when the vector q is aligned to h . The optimal filter is thus ˜ Q ( f ) ∝ R IJ ( f ) S h ( f ) N ( f ) (3.15)and the signal-to-noise achieved is SNR = 2 T (cid:90) ∞ df |R IJ ( f ) | S h ( f ) N ( f ) . (3.16)Other expressions in the literature do not have the factor 2. This is either because theirintegral ranges from −∞ to ∞ or due to the fact they are considering an auto-correlation ofa channel I , in which case N ( f ) = 2 N II ( f ) . The filter presented in the previous section has one limitation: It can only optimize the detec-tion of a single, global gravitational wave background. In particular, it can not accommodatefor the presence of multiple components with unknown relative amplitudes, or with ampli-tudes known up to some uncertainty. We now propose a new filter capable of handling thisgeneralized SGWB.We study explicitly the case of a two-components SGWB S h ( f ) = S p ( f ) + αS a ( f ) , (3.17)– 5 –here we want to optimize the filter ˜ Q ( f ) for the measurement of a primordial SGWB S p ( f ) ,in the presence of an astrophysical contribution with known shape S a ( f ) and unknown ampli-tude α , which is the only free parameter of the model. We comment on models with non-linearfree parameters in Appendix B. We allow for some external information on α to be available.We indeed assume its expected value to be ¯ α and its variance σ . The lack of such externalinformation is represented by the limit σ → ∞ .It is natural to amend the cross-correlation estimator to remove the expected contribu-tion from astrophysical sources y IJ = x IJ − ¯ α T (cid:90) ∞ df S a ( f ) R IJ ( f ) Q IJ ( f ) . (3.18)Note that the noise bias removed in Eq. (3.5) is not any different from the removal of anastrophysical component with σ = 0 .In this more general data model, the optimal filter becomes ˜ Q ( f ) ∝ R IJ ( f ) S p ( f ) N ( f ) − R IJ ( f ) S a ( f ) N ( f ) 2 T (cid:82) ∞ df (cid:48) |R IJ ( f (cid:48) ) | S a ( f (cid:48) ) S p ( f (cid:48) ) N − ( f (cid:48) ) σ − + 2 T (cid:82) ∞ df (cid:48) |R IJ ( f (cid:48) ) | S a ( f (cid:48) ) S a ( f (cid:48) ) N − ( f (cid:48) ) . (3.19)and the corresponding SNR is SNR = 2 T (cid:90) ∞ df |R IJ ( f ) | S p ( f ) N ( f ) − (cid:2) T (cid:82) ∞ df |R IJ ( f ) | S a ( f ) S p ( f ) N − ( f ) (cid:3) σ − + 2 T (cid:82) ∞ df |R IJ ( f ) | S a ( f ) N − ( f ) . (3.20)We report the derivation in Appendix A, where we consider the even more general case of anarbitrary number of contributions to S h ( f ) .In both Eq. (3.19) and (3.20), the first term corresponds to the standard filter and SNR,reported in Eqs. (3.15) and (3.16). The second term is the correction that accounts for thepresence of the astrophysical component in the SGWB. Note that in Eq. (3.20) this term isalways negative, reflecting the intuition that the presence of astrophysical sources can onlydegrade the SNR.Looking at the numerator more closely, the integral (and thus the SNR degradation)is maximum when the primordial and astrophysical components have the same frequencydependence. When this happens the SNR becomes SNR = 2 T (cid:82) ∞ df |R IJ ( f ) | S p ( f ) N − ( f )1 + σ T (cid:82) ∞ df |R IJ ( f ) | S a ( f ) N − ( f ) for S a ( f ) ∝ S p ( f ) , (3.21)which is zero for σ → ∞ : When the primordial and astrophysical SGWB have the samefrequency dependence, nothing can be said on the primordial component if no external infor-mation on the astrophysical one is available.Focusing now on the denominator of the second term in Eq. (3.20), it represents thetotal information about the amplitude of the astrophysical SGWB. It is clearly separated intothe external and the internal contribution, with the latter always overtaking the former forsufficiently long observational times. Building on the approach illustrated in the previous section, we now extend the discussion tothe reconstruction of the primordial SGWB without assuming a template for it. More details– 6 –re provided in Appendix C, together with the general expression for an arbitrary number offoreground components.We consider the following quantity, z IJ ( f ) ≡ T R IJ ( f ) (cid:90) ∞ df (cid:48) (cid:0) d ∗ I ( f ) d J ( f (cid:48) ) + d I ( f ) d ∗ J ( f (cid:48) ) (cid:1) δ T ( f − f (cid:48) ) − N IJ ( f ) R IJ ( f ) − ¯ αS a ( f ) . (4.1)It is essentially an elaboration on Eq. (3.5), without the integration over f . The two termsinside the integral (instead of one) have the only effect of folding the integral, which rangesover positive values of f (cid:48) . z IJ ( f ) can already be considered a template-free reconstruction,as the expected value is (cid:104) z IJ ( f ) (cid:105) = S p ( f ) : All the terms outside of the integral just removeadditive and multiplicative biases. Its variance is V ( f, f (cid:48) ) = N ( f ) δ ( f − f (cid:48) )2 T |R IJ ( f ) | + σ S a ( f ) S a ( f (cid:48) ) . (4.2)and its inverse is equal to F ( f, f (cid:48) ) = 2 T |R IJ ( f ) | N ( f ) δ ( f − f (cid:48) ) − (cid:2) T |R IJ ( f ) | S a ( f ) N − ( f ) (cid:3) (cid:2) T |R IJ ( f (cid:48) ) | S a ( f (cid:48) ) N − ( f (cid:48) ) (cid:3) σ − + 2 T (cid:82) ∞ df (cid:48)(cid:48) |R IJ ( f (cid:48)(cid:48) ) | S a ( f (cid:48)(cid:48) ) N − ( f (cid:48)(cid:48) ) . (4.3)In the Gaussian approximation F coincides with the Fisher matrix, the information availablefor any SGWB model for the given experimental configuration. In any case, it is related tothe SNR on a specific model in Eq. (3.20) by SNR = (cid:90) ∞ df (cid:90) ∞ df (cid:48) F ( f, f (cid:48) ) S p ( f (cid:48) ) S p ( f ) . (4.4) z IJ is typically very noisy and it should be used as an intermediate step towards itsprojection onto a subspace defined by a set of model spectra. These modes can be motivatedby the theory (e.g., SGWB spectra from a class of physically motivated models), by a targetproperty of the reconstructed spectrum (e.g., polynomials of varying nature or some typeof smooth functions to avoid sharp features) or by the constraints that the experiment canprovide over them. The principal component analysis (PCA) is part of this last category andits usage in this context was already proposed by [34]. Our explicit expression for the Fisherinformation allows to better understand (and exploit) its structure, and to properly computeits eigenvectors and eigenvalues (see Appendix C). PCA consists in selecting the k largesteigenvalues of F , λ i , and projecting z IJ onto the corresponding eigenvectors, v i ( f ) , whichrepresent the uncorrelated modes in our measurement . The estimated spectrum is ˆ S p ( f ) = k (cid:88) i =1 v i ( f ) (cid:90) ∞ df (cid:48) v i ( f (cid:48) ) z IJ ( f (cid:48) ) , (4.5)and has a total SNR equal to SNR = k (cid:88) i =1 λ i (cid:20)(cid:90) ∞ df (cid:48) v ( f (cid:48) ) z IJ ( f (cid:48) ) (cid:21) . (4.6) λ i and v i ( f ) can be equivalently computed either from the Fisher or the covariance matrix. The covariancehas an easier expression while the Fisher matrix handles more naturally σ → ∞ . – 7 –hese expressions are particularly simple compared to the corresponding expressions for anyother basis of models on which to project z IJ (compare with Eqs. (C.5-C.6)). Moreover, F is a diagonal matrix plus a (very-)low-rank correction, which makes the computation of theeigenvectors cheap even for a large number of frequencies. However, we believe that the PCAis not very interesting in this context. The reason is that this structure of F also makes theeigenvectors very similar to delta functions, which in turn makes the PCA only marginallydifferent from the selection of a frequency interval in the trough of the noise curve. Projectingon a different basis of models is reasonably simple anyway: It is essentially a generalized leastsquared estimation (see Appendix C for more details and Section 5.4 for an example). In this section we show an example of how the formalism of the previous section can beapplied. We forecast the performance of the LISA mission attempting to detect a specificmodel of primordial SGWB in the presence of an astrophysical signal.
The Laser Interferometer Space Antenna (LISA) [38] is the most advanced project for a grav-itational wave antenna in space. It is an ESA L-class mission with a NASA partnership andit is planned for launch in the early/mid 2030s. The experiment consists in three spacecraft,2.5 million km apart, in a triangular configuration. Laser links between the spacecrafts willprovide three interferometry measurements which will monitor the relative distance betweenthe test masses in the spacecrafts, allowing to probe the space-time distortion due to incominggravitational waves. Each of the three pair of arms provide an interferometry measurement—a time series as the one in Eq. (3.1). These are the so-called
XY Z channels. These signalsare correlated but, thanks to the symmetry of the configuration, it is easy and natural toextract their eigenmodes, dubbed
AET channels, given by d A = ( d X − d Y + d Z ) / √ (5.1) d E = ( d X − d Z ) / √ (5.2) d T = ( d X + d Y + d Z ) / √ . (5.3)While A and E have identical noise properties and response function, those of the T channelare very different and make this latter channel much less sensitive than the former ones (seeFigure 1). We also convert the A/E noise PSD to energy density Ω GW with Eq. (2.8) anddisplay it in Figure 3. The calculation of the noise curves as well as the response functions follows [39] and make the same simplifying assumptions. For example, we assume that thereis no gap in the data, the noise is perfectly stationary and we ignore the time dependence ofthe response function due to the orbital motion of the spacecrafts. We refer the reader to thisarticle for more details. We consider two types of of backgrounds of astrophysical origin. The first one, is sourcedby the incoherent superposition of the emission of compact binaries that are not individu-ally resolved by the experiment—mostly stellar-origin black holes and neutron star binaries For the response functions we use the tabulated numerical values that the authors made available at https://doi.org/10.5281/zenodo.3341817 – 8 – f [Hz]10 ( f ) f [Hz]10 N ( f ) f [Hz]10 f S n ( f ) AA , EETT
Figure 1 . The main properties of the A , E and T channels of LISA, the experimental configurationconsidered in the demonstration in Section 5. Left:
Response function.
Center:
Noise power spectraldensity.
Right:
Characteristic strain sensitivity. (BBH+BNS). This signal is well approximated by a power law [4] Ω GW ( f ) h = Ω ∗ (cid:18) ff ∗ (cid:19) / , ( BBH + BN S ) (5.4)The reference value we use is Ω ∗ = 8 . × − at f ∗ = 25 Hz [40] but we consider it anunknown parameter to marginalize over. We assume that a prior on it is available, it iscentered around the reference value and has a standard deviation given by σ BBH+BNS
The second type of astrophysical background that we consider is sourced by unresolvedgalactic binaries (UGB). We model its contribution with [41, 42] S ( f ) = Af − / e − f α + βf sin κf [1 + tanh( γ ( f k − f ))] Hz − , ( U GB ) (5.5)which we convert to Ω GW with Eq. (2.8). We fix the free parameters to the following values,which refer to a 4 year-long LISA mission, α = 0 . , β = − , κ = 521 , γ = 1680 and f k = 0 . [42]. Our reference value for the normalization factor A is × − but, alsohere, we treat it as a free parameter with some prior constraint centered at the reference valueand a standard deviation given by σ UGB .Finally the primordial SGWB that we consider is the AX1 model of [43]—which isproduced by a spectator axion-SU(2) model. The origin and properties of such a backgroundare, however, irrelevant for our analysis. We choose this model only for illustration purposesand are motivated mostly by its amplitude, which makes such a background within reach forthe instrumental configuration of our choice.We display the three components of the measured background in Figure 3. Their dif-ferent frequency dependence is the main feature that our methodology exploits, while theiramplitudes constitute the free parameters of the model.
We now apply the formalism of Section 3.2 to the LISA A channel and study the new filterand SNR. The results for the E channel would be of course identical, while the T channel isnoise dominated and we do not report results about it.We consider values of σ UGB ranging from − and − and values of σ BBH+BNS rangingfrom − and − . As we will see, the lower values correspond to the case where the externalinformation completely constrains the astrophysical component, the higher value is equivalentto the case where no external information is available. We stress that the range of values for– 9 – BBH + BNS S N R UGB = 0.0001
UGB = 0.001
UGB = 0.01 10 BBH + BNS S N R UGB = 0.0001
UGB = 0.001
UGB = 0.01
Figure 2 . Left:
Signal-to-noise ratio from either the A or E channel, for varying values of σ UGB and σ BBH+BNS , computed with Eq. (3.20)—actually Eq. (A.11), because we are dealing with multipleastrophysical sources. For reference, the black dashed line reports the SNR obtained with the standardformula Eq. (3.16) ignoring the astrophysical components.
Right:
Signal-to-noise ratio from the A and E channels combined. The dashed line reports what one would get with a naive co-addition ofthe SNR of the individual channels (i.e. multiplying them by √ ). σ UGB and σ BBH+BNS is chosen only to best illustrate the phenomenological properties of ourformalism and does not stem from an instrumental or theoretical forecast. Also the inclusionof this level of UGB is not fully realistic, as this signal can be removed exploiting its timedependence owing to the motion of the spacecrafts [44].
The filter
We start by focusing on the effect that the astrophysical components have onthe filter Eq. (3.19)—and Eq. (A.10). The three panels of Figure 3 refer to three differentamounts of external information on the extra-galactic binaries, while the three colors representdifferent priors on the galactic binaries. The dashed line reports the standard filter, whichonly performs inverse noise co-addition: it is always positive and significantly larger than zeroonly between 0.4 mHz and 10 mHz. The blue line in the top panel represents the case inwhich the amplitude of both the astrophysical backgrounds is very constrained by externalinformation. This case reduces to the standard filter, as mentioned in Section 3.2. When σ UGB increases to 0.001 (orange line) the filter give less weight to the frequencies around 1mHz in order to reduce the response to the UGB signal—which peaks at these frequencies,see the figure on the left. When σ UGB reaches 0.01 (green line) the amplitude of the signalfrom galactic binaries is so uncertain that the filter is required to have zero response to theUGB signal shape. This can only be achieved with a negative region in the filter: The trougharound 1 mHz that the green lines have in all the panels.A completely analogous discussion can be done about the BBH and BNS signal, whoseamplitude with respect to the other components increases with frequency—even though it isshadowed by the noise above 10 mHz. This is the reason why the filter with σ BBH+BNS = 0 . (middle panel) gives slightly less weight compared to the case with σ BBH+BNS = 0 . (toppanel) and the same region becomes negative for σ BBH+BNS = 0 . .As a final remark, we note that the peak of the filter increases when the amount ofexternal information is reduced. This is a consequences of the fact that when displaying these– 10 – f [Hz]10 h G W ( f ) Noise A/EUGBBBH + BNSAX1 01 Q ( f ) BBH + BNS = 0.001
UGB =0.01
UGB =0.001
UGB =0.0001 Q ( f ) BBH + BNS = 0.01 f [Hz]025 Q ( f ) BBH + BNS = 0.1
Figure 3 . Left:
Components of the SGWB considered in application of the methodology. The noiseequivalent Ω GW for the A (or E ) channel is also reported. Right:
Shape of the filter functions fordifferent values of σ UGB (color) and σ BBH+BNS (subplots). For comparison, the standard filter is thedashed back line. filters we impose that they have the same response to the target primordial signal . Therefore,the decreased (or negative) weight in the frequencies dominated by the astrophysical emissionis compensated with an increased weight around the minimum of the astrophysical emission. The SNR
We now study the SNR achieved for the same configuration and external in-formation about the astrophysical components. The result is reported in Figure 2, wherewe also draw the standard SNR (black dashed line) computed with Eq. (3.16) ignoring thepresence of astrophysical sources. First note that this value coincides with the SNR obtainedfor σ BBH+BNS = 0 . and σ UGB = 0 . , when the amplitude of the astrophysical signals iscompletely constrained by the external prior and therefore there is no information loss whilemarginalizing over them (blue line, left end). Looking at the other extreme, if the ampli-tude of the astrophysical components is completely unconstrained ( σ BBH+BNS = 0 . and σ UGB = 0 . , green line, right end), the rejection of their signal severely degrades the SNRfrom 12 to 2.4, which corresponds to a 96% information loss. This should be largely ascribedto the BBH+BNS, as even with σ BBH+BNS = 0 . (blue line) the information loss is still93%. The reason is simple: The frequency dependence of the primordial signal is much moreorthogonal to the one of the UGB than the one of the BBH+BNS. The only purpose of this section is to stress a supposedly obvious fact that is however easyto forget: Using the same prior information across multiple measurements correlates theconstraints they produces, and they can not be combined as if they were independent.One of the most attractive properties of working with the AET channels compared tothe XYZ channels is the fact that the cross-channel power is zero and they are statisticallyequivalent to independent experiments, at least in our idealised treatment. Moreover, the Aand E channels have the same noise level and, therefore, combining them effectively doublesthe statistical information provided by the individual channel and, in the standard case Eq.(3.16),
SNR A + E = SNR A + SNR A = 2SNR A . We remind that the normalization of the filter is arbitrary and does not affect its optimality – 11 – f [Hz]10 h G W ( f ) AX1 x 1 x 100 Single simulationMean simulations68% simulations Single simulation data
Figure 4 . Template-free reconstruction. The thick blue line is the target signal. The orangerepresents the estimators that assume the same σ UGB and σ BBH+BNS used for simulating the data,while the green assumes one hundred times larger values. The solid lines are the median of simulated reconstructions and the shaded areas covers from the 16 th to the 84 th percentile. As anexample, we also show the data (dots) and the reconstruction (dashed lines) for a single simulation. The formalism in Section 3.2, however, can accommodate external information, whichmay play a significant role in
SNR A and SNR E computed with Eq. (3.20). When this happensthe two SNRs are correlated and can not be added in quadrature. Instead of accounting forthe correlation in their coaddition, it is easier to properly compute the joint SNR. It can beshown easily that the optimal combination of the two channels boils down to the averageof the two auto-correlations, which has the same expected value and half the variance theindividual channels and is therefore equivalent to the A (or E ) channel with a N dividedby √ . SNR A + E can be computed with the Eq (3.20) applied to A (or E ) but multiplyingevery integral by 2. Multiplying SNR A by two altogether, doubles the external information σ − —which, on the contrary, stays constant in SNR A + E .In the right panel of Figure 2 we compare SNR A + E with √ A . The two coincidewhenever the the σ s are either very low or very high (the endpoints of the blue and greenlines). When a σ is very low, the correction term in Eq. (3.20) is negligible in both cases,while when the σ is very high the correction term is not affected by its exact value.Summarizing, the naive co-addition of the SNRs from independent measurements issignificantly incorrect when the external information plays a significant but not overwhelmingrole. Note that this example of combination of the A and E channel completely analogous tothe combination of separate frequency bins. Naively co-adding the SNR from arbitrarily smallbins produces a total SNR arbitrary close to the standard SNR with no contaminants—whichis of course unphysical (and wrong). Without assuming prior knowledge on the primordial SGWB, we perform its reconstructionfrom simulated data of the A channel following the prescription of Section 4 (and AppendixC). We simulate z AA at linearly spaced frequencies, emulating what we would get if we wereto compute it from the FFT of d A . The signal in z AA is the same primordial SGWB ofthe previous sections. On the top of it we add Gaussian noise generated according to the– 12 –ovariance matrix Eq. (4.2) assuming σ UGB = 0 . and σ BBH+BNS = 0 . . We create andanalyze independently realizations. The reconstruction is done using 8 top-hat functionswith equal logarithmic width, covering the range from 0.1 mHz to 0.1 Hz. More sophisticatedsets of functions are possible, of course, but this simple choice already allows to show someof the main features of our approach. The variance assumed by the template-free estimatoris either the same of the simulated data, or one that assumes one hundred times larger σ UGB and σ BBH+BNS . This latter case essentially coincides with the absence of external informationand, therefore, its inverse-variance down-weights (or filters) more aggressively the signals withthe frequency dependence of the astrophysical components.The results of the analysis are shown in Figure 4, which shows all the quantities ( z AA included) in terms of energy density spectrum Ω GW . Focusing on the single simulation rep-resented by the dots, it clearly shows an over-subtraction of the galactic signal and an under-subtraction of the extragalactic astrophysical signal. Our estimator differs from a simplebinning because it not only does an inverse-noise-weighted average inside the bins, but alsodown-weights the spectral shapes that match the ones of the astrophysical signals—the larger σ , the stronger the down-weighting. The reconstructions (dashed lines) are indeed closer tothe target signal than the original data. In particular, the estimator that assumes larger σ (green) filters out the spectral shapes of astrophysical origin and, therefore, it seems tomitigate better the under- and over-subtraction of the astrophysical signals. However, whenmany simulations are averaged (solid lines), filtering translates into a bias because the samemodes are suppressed from the target signal. Instead, when the estimator uses the correctvalues of σ the mean reconstruction is unbiased, even though in the single simulation theastrophysical residuals do not seem to be removed as effectively—which is also reflected inthe larger scatter of the simulations (orange boxes, compared to the green boxes).As a last remark, we want to stress that inferring uncertainties from the scatter ofindividual frequency bins across simulations (the shaded areas in Figure 4) can be highlymisleading: They represent only the diagonal of the covariance matrix, while the bins arehighly correlated. Once the full covariance is taken into account (see Appendix C) one realisesthat the estimator that uses the correct values of σ is not only unbiased even on the singlesimulation but also it is truly minimum variance. The scatter in the simulations is dominatedby modes that have the same frequency dependence of the astrophysical signals because thesemodes are highly uncertain. When they are completely filtered out not only the scatter butalso the signal in those modes is removed and therefore their relative uncertainty is infinity. The existence of some form of SGWB is well motivated both from the theory and the directGW observations of the past years. The standard detection techniques have focused ondistinguishing a template SGWB signal from the noise in the auto- or cross-correlation ofinterferometric signals. From the data analysis point of view, there is freedom on the filterfunction involved in the correlation and the optimal choice is equal to the template signaldivided by the noise power variance—both diagonal in the frequency-domain.We extend this approach to the case in which other sources are present in additionto the target SGWB. We derive the filter that optimally balances between the inverse-noiseweighting done by the standard filter and the marginalization over the unwanted components,taking into account possible external priors. Both the filter and the corresponding SNR arealmost as simple as the the standard ones—and reduce to them either if the amount of– 13 –xternal information is very large, or if the spectral shape of the unwanted components isvery different from the target SGWB. In particular the optimal filter and SNR have closed-forms and all the terms are easy to compute and interpret. It is extremely simple, for example,to forecast the sensitivity of an experimental configuration to a primordial SGWB producedby an exotic inflationary model while, at the same time, marginalizing over a background Ω GW ∝ f / produced by unresolved binary black holes and neutron stars [see 43]. We applythe formalism to the LISA mission and show that neglecting the marginalization over the theastrophysical signals—as done by the standard estimator—can grossly overestimate the SNR.We remind that in the text we have used the terms “primordial” or “astrophysical SGWB”only for illustration purposes. The approach is perfectly applicable to contexts where thetarget signal has astrophysical origin or where both the target and unwanted component(s)are primordial.Our methodology is derived and applied in the context of isotropic, unpolarized back-ground signals, but it can be easily generalized to polarized and anisotropic signals. The onlyimportant caveat is that it is based on the assumption (standard in the literature) that thevariance in the auto- or cross-correlation is dominated by the noise at all frequencies. Thisassumption is key in obtaining nice and simple closed-form for the optimal filter and SNR,but may be too stringent in very high SNR settings.Compared to techniques that blindly attempt to reconstruct the SGWB, the new filterwe have illustrated does require a template of the target signal. However, the filter can berelevant also to procedures that perform a more flexible SGWB reconstruction. Take for ex-ample [33], the authors perform a template-free reconstruction of a primordial SGWB whilemarginalizing over a set of templates with uncertain amplitude (representing the same astro-physical components we considered in Section 5). Their MCMC samples both the parametersof their non-linear primordial SGWB model and the amplitude of the astrophysical templates.If the interest lies in the former and the latter are only marginalized over, this marginalisa-tion can be done analytically during the MCMC using the filter in Eq. (3.19) (or analogousexpressions) with a signal defined by the current value of the target signal parameters in theMCMC chain. This allows to achieve the same result with reduced computational cost, as theparameters related to the astrophysical background are not explored by the Markov chain.Moreover, we have shown how our methodology is applicable also to a template-freereconstruction of the SGWB. Provided a basis for the admissible models, the method projectthe data onto the space spanned by these models, balancing noise-weighting and foregroundremoval. In our example we consider simple top-hat functions, but any choice is admitted—such as polynomials, harmonic functions, or a basis of the possible SGWB produced by a classof inflationary models. Albeit much richer than a single template, this type of procedure allowsonly linear models, but this class of parametrizations may turn out to be sufficiently flexiblefor many applications and fit the (simulated) data with an accuracy comparable to a moreinvolved non-linear fit. Acknowledgments
The author is thankful to Paolo Campeti, Eiichiro Komatsu and Carlo Baccigalupi for use-ful comments and discussions. The author acknowledges support from the ASI-COSMOSnetwork ( ) and the INDARK INFN Initiative ( web.infn.it/CSN4/IS/Linea5/InDark ). – 14 –
Derivation of the new filter
We assume the presence of a primordial SGWB and a set of astrophysical backgrounds, dueto different populations of unresolved sources, S h ( f ) = S p ( f ) + (cid:88) i α i S a i ( f ) . (A.1)We optimize the filter for the primordial SGWB but make no use of its properties in thederivation. Likewise, the astrophysical origin of the unwanted components plays no role. Theonly relevant assumption for what follows is the fact that all the contributions have knownfrequency dependence and possibly unknown amplitude. For example, all the models thatpredict an f − / dependence can be accommodated (or approximated) in this formalism.In the notation introduced in Section 3.1, the signal can be written as h = p + Aα , (A.2)where p is defined similarly to Eq. (3.12). The same is done for the astrophysical sources,which are further collected in the columns of A . Their amplitudes are free parameters andare collected in a vector α , with expected value ¯ α and covariance Σ .We are interested in optimizing the filter for the detection of p . The cross-correlationcan be corrected for the contribution of the astrophysical SGWB, y IJ = x IJ − q · A ¯ α , (A.3)so that its expected value depends only on the primordial SGWB, (cid:104) y IJ (cid:105) = q · p . (A.4)The variance is (cid:104) y IJ (cid:105) − (cid:104) y IJ (cid:105) = q t q + q t A Σ A t q , (A.5)where a t b ≡ a · b . The new expression for the SNR is therefore SNR = ( q t p ) q t q + q t A Σ A t q . (A.6)Also in this case, the SNR is independent of the normalization of q . We choose it such that q t p = γ, (A.7)where γ is an arbitrary constant. We now minimize the denominator of Eq. (3.20) under thecondition in Eq. (A.7), which can be readily done using Lagrange multipliers. The gradientof the Lagrangian is ∇ q L ( q , λ ) = 2( + A Σ A t ) q + λ p . (A.8)Since the matrix in parenthesis is positive-defined, there is only the following stationary pointand this point is a minimum q ∝ ( + A Σ A t ) − p . (A.9) If useful, a free parameter can be introduced also in front of the primordial contribution. Nothing wouldchange in the expression of the optimal filter and its derivation. The only difference would be the multiplicationof the SNR by the expected value of such additional parameter. – 15 –his expression is practically inconvenient because it involves the inversion of a large anddense matrix. Therefore, we use the Woodbury identity to re-express Eq. (A.9) as q ∝ p − A (Σ − + A t A ) − A t p , (A.10)which is the final form of our filter. The corresponding signal-to-noise ratio is SNR = p t p − p t A (Σ − + A t A ) − A t p . (A.11)For the case of one single astrophysical component, we report the explicit expression of thisfilter and SNR in Eqs. (3.19) and (3.20). B Non-linear SGWB models
Many models of both cosmological and astrophysical SGWB contain non-linear parameters.Therefore, we generalize our model as follows S h ( f ) = S p ( f ) + (cid:88) i α i S a i ( f, β i ) . (B.1)where β i is a non-linear parameter in the model of the i -th astrophysical component. ByTaylor expanding the data model to first order S h ( f ) (cid:39) S p ( f ) + (cid:88) i α i S a i ( f, ¯ β i ) + (cid:88) i ¯ α i ∂S a i ∂β i ( f, ¯ β i ) (cid:0) β i − ¯ β i (cid:1) . (B.2)it becomes linear in the free parameters and therefore compatible with the assumptions madein the Appendix A: A acquires extra columns given by all the ¯ α i ∂S ai ∂β i ( f, ¯ β i ) , and α acquiresextra rows given by all the β i − ¯ β i coefficients. Of course, Σ acquires the same number ofextra rows and columns, that can be used to accommodate external information about thethe non-linear parameters as well as their correlation with the linear ones. Note that, if theyare infinity (i.e. no external information is available on the non-linear parameters), the ¯ α i factor multiplying the i -th new column of A is irrelevant: Only the frequency dependence ofthe new column matters, not its normalization.This natural extension of the formalism to non-linear parameters should, however, beused with care when analyzing data (or simulations). First, one should consider if the modelis overly complex for the experimental configuration being analyzed. If the uncertainty onthe linear parameters α i is already very large, there is no point in refining the model of the i -th foreground. Second, even if the constraints on those parameters are good, those on thenon-linear parameters might be loose. As a result, α i S a i ( f, ¯ β i ) + ¯ α i ∂S ai ∂β i ( f, ¯ β i ) (cid:0) β i − ¯ β i (cid:1) forthe best fit α i and β i may be quite far from what one would get form the fully non-linear fitof α i S a i ( f, β i ) —meaning that the linear Taylor expansion Eq. (B.2) is insufficient to describethe behaviour of the i -th astrophysical component in the range of plausible noise realizations.This occurrence may not be easy to know in advance, it could even depend on the specificnoise realization. This may or may not be a problem for the constraints on the primordialcomponent of the SGWB, depending on how degenerate S p ( f ) is with the difference betweenthe linearized and fully non-linear fit. The degeneracy is typically unknown from the onset It can also be a vector of parameters, the generalization is straightforward. – 16 –nd its evaluation would require to compute the non-linear fit, which would make the linearapproximation of limited interest. Third, also the high signal-to-noise has a caveat: It stillrequires to guess ¯ β i . If ∂S ai ∂β i varies noticeably between ¯ β i and the true value of β i , thelinearized and non-linear best-fit could be significantly different and bias the constraints onthe primordial SGWB.On the other hand, the applicability range of the linearized data model is much wider inthe realm of forecasting. The biases related to the noise realization and the imperfect choiceof the reference ¯ β i are indeed under control. Moreover, if the quantity of interest is only theSNR, the value obtained with the linearized data model coincides with the Fisher estimate ofthe constraints from the non-linear fit. C A closer look to the template-free reconstruction
We start from the estimator in Eq. (4.1), amended to account for an arbitrary number ofastrophysical components z IJ ( f ) ≡ T R IJ ( f ) (cid:90) ∞ df (cid:48) (cid:0) d ∗ I ( f ) d J ( f (cid:48) ) + d I ( f ) d ∗ J ( f (cid:48) ) (cid:1) δ T ( f − f (cid:48) ) − N IJ ( f ) R IJ ( f ) − (cid:88) i ¯ α i S a i ( f ) . (C.1)Its covariance is equal to V ( f, f (cid:48) ) = N ( f ) δ ( f − f (cid:48) )2 T |R IJ ( f ) | + (cid:88) i σ i S a i ( f ) S a i ( f (cid:48) ) . (C.2)Since (cid:104) z IJ ( f ) (cid:105) = S p ( f ) , V is the frequency-frequency covariance for single-sided primordialsignals. The first term is the instrumental noise, while the second is the variance due to thefact that the normalizaiton of the astrophysical component a i is known up to a standarddeviation σ i .From now on we assume that z IJ ( f ) is Gaussian-distributed. This is an approximationthat may be oversimplistic at the full sampling rate of the experiments, but it becomesrealistic if z IJ ( f ) represents an average over several independent periods or the mean of manyneighbouring frequencies. In this approximation the Fisher information matrix is the inverseof the variance. To obtain readable expressions we define N ≡ N ( f ) δ ( f − f (cid:48) ) / T |R IJ ( f ) | and Σ ≡ diag( σ i ) and collect all the astrophysical signals S a i in the columns of A . Notethat the operator A is not exactly the same of Appendix A: The basis of the frequencydomain is different. We will also use a different (and more standard) inner product: a · b ≡ (cid:82) ∞ df a ( f ) b ( f ) . In this notation V ( f, f (cid:48) ) = N + A Σ A t (C.3)and its inverse can be obtained using the Woodbury identity F ( f, f (cid:48) ) = N − − N − A (cid:0) Σ − + A t N − A (cid:1) − A t N − . (C.4) It might be useful to remember that (cid:104) n ∗ I ( f ) n J ( k ) n I ( f (cid:48) ) n ∗ J ( k (cid:48) ) (cid:105) − (cid:104) n ∗ I ( f ) n J ( k ) (cid:105)(cid:104) n I ( f (cid:48) ) n ∗ J ( k (cid:48) ) (cid:105) = 14 [ δ ( f − f (cid:48) ) δ ( k − k (cid:48) ) N II ( f ) N JJ ( k )++ δ ( f + k (cid:48) ) δ ( f (cid:48) + k ) N IJ ( f ) N IJ ( f (cid:48) )] – 17 –his expression reduces to the Eq. (4.3) in presence of only a single astrophysical source.Note that, in spite of being very large, this matrix is easy to handle because it is a diagonalmatrix plus a low-rank correction. When no external information is available ( Σ − = 0 ), theterm in parenthesis has the only role of making the A columns N − -orthonormal, so that thespectrum of F on the space spanned by the foregrounds is zero (i.e., there is no informationon any linear combination of the columns of A ).Using F as the inverse variance, we can also compute (cid:82) ∞ df (cid:82) ∞ df (cid:48) F ( f, f (cid:48) ) z IJ ( f (cid:48) ) z IJ ( f ) to try to detect the presence of a signal in z IJ that can not be explained as a statisticalfluctuation of the noise or the astrophysical signal. Unfortunately, this can hardly be usefulbecause this quantity is dominated by the immense number of low-SNR modes, which wouldlikely shadow even fairly visible signals. We can select a family of modes { m i ( f ) } i =1 ,...,k that we wish to use as a basis for the k -dimensional sub-space of admissible primordialbackgrounds. These modes can be used to probe the SNR on them one-by-one with Eq. (4.4).It is probably more interesting, however, to project z IJ on the entire subspace generated bythe models of our choice and to compute the total SNR. Collecting the basis vector in thecolumns of the matrix M , the minimum-variance fit to the data is ˆ S p ( f ) = M (cid:0) M t F M (cid:1) − M t F z , (C.5)and the SNR achieved is SNR = z t F M (cid:0) M t F M (cid:1) − M t F z . (C.6)These are the familiar expressions of the generalized least squared estimator, which arise whenestimating the parameters of a linear model under a constraint of minimum variance. Notethat, in general, the reconstructed amplitudes of the modes { m i ( f ) } i =1 ,...,k —estimated usingEq. (C.5) without the leading M —are correlated: Their covariance is (cid:0) M t F M (cid:1) − .We conclude this appendix with some computational remarks that can be relevant inthis context when the number of sampled frequencies n is large. PCA computes the k eigen-vectors { v i ( f ) } i,...,k —collected in V —and the corresponding eigenvalues { λ i } i,...,k . Whenthe eigenvectors are employed for a template-free reconstruction, we get particularly simpleexpressions ˆ S p ( f ) = V V t z (C.7)and ˆ S p ( f ) = z t V diag( λ , . . . , λ k ) V t z , (C.8)which coincide with the expressions in Eq. (4.5) and (4.6). As we said in Section 4, we do notbelieve that PCA is interesting for the type of application discussed in this paper, but still wewant to highlight that it is possible to do it when n is very large, even though an explicit calcu-lation of the matrix and its decomposition would respectively require O ( n ) storage locationsand up to O ( n ) operations. The most practical solution is to down-sample the frequenciesuntil the brute-force calculation can be afforded, but the PCA can be performed without anydown-sampling in order O ( k × n ) operations because the matrix (Fisher or covariance) is adiagonal matrix plus a low-rank correction. Such a simple structure allows to store the ma-trix with only few O ( n ) locations—the spectrum of the noise and the astrophysical templates,whose number is of order unity—and to apply the matrix on a vector in O ( n ) operations:It consists of few vector-vector and scalar-vector operations. We are not interested in thefull set of eigenvectors and eigenvalues of the matrix, just in the k best constrained by our– 18 –xperiment. The eigenvectors with the highest eigenvalues can be computed in O ( k ) appli-cations of the matrix—plus other O ( k ) vector-vector products—using the Lanczos method.The overall O ( k × n ) scaling of this PCA poses no threat for any sensible value of n and k .Note, however, that this scaling is accurate only for k (cid:28) n . When k is comparable with n ,other O ( k ) calculations internal to the Lanczos method start to be significant. The totalscaling is still better than O ( n ) but if n is small the prefactors can be significant. Even moreimportant, the explicit solution to the eigenvalue problem internally uses matrix-matrix op-erations, compared to the vector-vector operations inside the application of our sparse Fishermatrix. The former run much faster than the latter on modern CPUs (and GPUs), so thetime-to-completion of an iterative O ( k × n ) algorithm may result comparable (or even higher!)than a O ( n ) algorithm for small values of n .Lastly, we want to comment on the inversion of the covariance matrix. The analyticalexpression in Eq. (C.4) is reasonably simple to implement and handle (most important, itis diagonal with a low-rank correction, as the covariance matrix). Still one may prefer towork only with the covariance matrix instead of the Fisher matrix, in order to avoid thecomputation of (cid:0) Σ − + A t N − A (cid:1) − . Eqs. (C.5) and (C.6) can be computed using only theability of applying V on a vector by casting every application of F on a vector (or matrix) intoan inverse problem with V as system matrix. If solved with the preconditioned conjugategradient (PCG), the solution is achieved at most in a number of iterations equal to thenumber of astrophysical components, thus with few O ( n ) operations. The reason is that theCG converges to the solution in a number of iterations lower than the number of distincteigenvalues of the system matrix. Once preconditioned with the inverse of the first term inEq. (C.2), the variance becomes the identity plus a low-rank correction and, therefore, thenumber of distinct eigenvalue is equal to this rank plus one. References [1] Michele Maggiore. Stochastic backgrounds of gravitational waves.
ICTP Lect. Notes Ser. ,3:397–414, 2001.[2] Nelson Christensen. Stochastic Gravitational Wave Backgrounds.
Rept. Prog. Phys. ,82(1):016903, 2019.[3] Chiara Caprini and Daniel G. Figueroa. Cosmological Backgrounds of Gravitational Waves.
Class. Quant. Grav. , 35(16):163001, 2018.[4] Tania Regimbau. The astrophysical gravitational wave stochastic background.
Res. Astron.Astrophys. , 11:369–390, 2011.[5] B. P. Abbott et al. GW151226: Observation of Gravitational Waves from a 22-Solar-MassBinary Black Hole Coalescence.
Phys. Rev. Lett. , 116(24):241103, 2016.[6] B.P. Abbott et al. Observation of Gravitational Waves from a Binary Black Hole Merger.
Phys. Rev. Lett. , 116(6):061102, 2016.[7] B.P. Abbott et al. Binary Black Hole Mergers in the first Advanced LIGO Observing Run.
Phys. Rev. X , 6(4):041015, 2016. [Erratum: Phys.Rev.X 8, 039903 (2018)].[8] Benjamin P. Abbott et al. GW170104: Observation of a 50-Solar-Mass Binary Black HoleCoalescence at Redshift 0.2.
Phys. Rev. Lett. , 118(22):221101, 2017. [Erratum: Phys.Rev.Lett.121, 129901 (2018)].[9] B.. P.. Abbott et al. GW170608: Observation of a 19-solar-mass Binary Black HoleCoalescence.
Astrophys. J. , 851(2):L35, 2017. – 19 –
10] B.P. Abbott et al. GW170814: A Three-Detector Observation of Gravitational Waves from aBinary Black Hole Coalescence.
Phys. Rev. Lett. , 119(14):141101, 2017.[11] B.P. Abbott et al. GW170817: Observation of Gravitational Waves from a Binary NeutronStar Inspiral.
Phys. Rev. Lett. , 119(16):161101, 2017.[12] Alessandra Buonanno, Gunter Sigl, Georg G. Raffelt, Hans-Thomas Janka, and Ewald Muller.Stochastic gravitational wave background from cosmological supernovae.
Phys. Rev. D ,72:084001, 2005.[13] Konstantin N. Yakunin et al. Gravitational Waves from Core Collapse Supernovae.
Class.Quant. Grav. , 27:194005, 2010.[14] Valeria Ferrari, Sabino Matarrese, and Raffaella Schneider. Stochastic background ofgravitational waves generated by a cosmological population of young, rapidly rotating neutronstars.
Mon. Not. Roy. Astron. Soc. , 303:258, 1999.[15] Quan Cheng, Yun-Wei Yu, and Xiao-Ping Zheng. Stochastic gravitational wave backgroundfrom magnetic deformation of newly born magnetars.
Mon. Not. Roy. Astron. Soc. ,454(3):2299–2304, 2015.[16] Alexei A. Starobinsky. A New Type of Isotropic Cosmological Models Without Singularity.
Adv. Ser. Astrophys. Cosmol. , 3:130–133, 1987.[17] Andrei D. Linde. A New Inflationary Universe Scenario: A Possible Solution of the Horizon,Flatness, Homogeneity, Isotropy and Primordial Monopole Problems.
Adv. Ser. Astrophys.Cosmol. , 3:149–153, 1987.[18] Alan H. Guth. The Inflationary Universe: A Possible Solution to the Horizon and FlatnessProblems.
Adv. Ser. Astrophys. Cosmol. , 3:139–148, 1987.[19] Viatcheslav F. Mukhanov and G. V. Chibisov. Quantum Fluctuations and a NonsingularUniverse.
JETP Lett. , 33:532–535, 1981.[20] M. Hazumi et al. LiteBIRD: A Satellite for the Studies of B-Mode Polarization and Inflationfrom Cosmic Background Radiation Detection.
J. Low Temp. Phys. , 194(5-6):443–452, 2019.[21] Kevork Abazajian et al. CMB-S4 Science Case, Reference Design, and Project Plan. 7 2019.[22] T.W.B. Kibble. Topology of Cosmic Domains and Strings.
J. Phys. A , 9:1387–1398, 1976.[23] Tanmay Vachaspati, Levon Pogosian, and Daniele Steer. Cosmic Strings.
Scholarpedia ,10(2):31682, 2015.[24] Saswat Sarangi and S.H.Henry Tye. Cosmic string production towards the end of braneinflation.
Phys. Lett. B , 536:185–192, 2002.[25] Nicholas T. Jones, Horace Stoica, and S.H. Henry Tye. The Production, spectrum andevolution of cosmic strings in brane inflation.
Phys. Lett. B , 563:6–14, 2003.[26] Luca Pagano, Laura Salvati, and Alessandro Melchiorri. New constraints on primordialgravitational waves from Planck 2015.
Phys. Lett. B , 760:823–825, 2016.[27] B.P. Abbott et al. Constraints on cosmic strings using data from the first Advanced LIGOobserving run.
Phys. Rev. D , 97(10):102002, 2018.[28] B.P. Abbott et al. GW150914: Implications for the stochastic gravitational wave backgroundfrom binary black holes.
Phys. Rev. Lett. , 116(13):131102, 2016.[29] Benjamin P. Abbott et al. GW170817: Implications for the Stochastic Gravitational-WaveBackground from Compact Binary Coalescences.
Phys. Rev. Lett. , 120(9):091101, 2018.[30] Zaven Arzoumanian et al. The NANOGrav 12.5-year Data Set: Search For An IsotropicStochastic Gravitational-Wave Background. 9 2020. – 20 –
31] R.w. Hellings and G.s. Downs. UPPER LIMITS ON THE ISOTROPIC GRAVITATIONALRADIATION BACKGROUND FROM PULSAR TIMING ANALYSIS.
Astrophys. J. Lett. ,265:L39–L42, 1983.[32] Zhen Pan and Huan Yang. Probing Primordial Stochastic Gravitational Wave Background withMulti-band Astrophysical Foreground Cleaning.
Class. Quant. Grav. , 37(19):195020, 2020.[33] Raphael Flauger, Nikolaos Karnesis, Germano Nardini, Mauro Pieroni, Angelo Ricciardone,and Jesús Torrado. Improved reconstruction of a stochastic gravitational wave backgroundwith LISA. 9 2020.[34] Mauro Pieroni and Enrico Barausse. Foreground cleaning and template-free stochasticbackground extraction for LISA.
JCAP , 07:021, 2020. [Erratum: JCAP 09, E01 (2020)].[35] Bruce Allen. The Stochastic gravity wave background: Sources and detection. In
Les HouchesSchool of Physics: Astrophysical Sources of Gravitational Radiation , pages 373–417, 4 1996.[36] Bruce Allen and Joseph D. Romano. Detecting a stochastic background of gravitationalradiation: Signal processing strategies and sensitivities.
Phys. Rev. D , 59:102001, 1999.[37] Michele Vallisneri and Chad R. Galley. Non-sky-averaged sensitivity curves for space-basedgravitational-wave observatories.
Class. Quant. Grav. , 29:124015, 2012.[38] Pau Amaro-Seoane, Heather Audley, Stanislav Babak, John Baker, Enrico Barausse, PeterBender, Emanuele Berti, Pierre Binetruy, Michael Born, Daniele Bortoluzzi, Jordan Camp,Chiara Caprini, Vitor Cardoso, Monica Colpi, John Conklin, Neil Cornish, Curt Cutler,Karsten Danzmann, Rita Dolesi, Luigi Ferraioli, Valerio Ferroni, Ewan Fitzsimons, JonathanGair, Lluis Gesa Bote, Domenico Giardini, Ferran Gibert, Catia Grimani, Hubert Halloin,Gerhard Heinzel, Thomas Hertog, Martin Hewitson, Kelly Holley-Bockelmann, DanielHollington, Mauro Hueller, Henri Inchauspe, Philippe Jetzer, Nikos Karnesis, Christian Killow,Antoine Klein, Bill Klipstein, Natalia Korsakova, Shane L Larson, Jeffrey Livas, Ivan Lloro,Nary Man, Davor Mance, Joseph Martino, Ignacio Mateos, Kirk McKenzie, Sean TMcWilliams, Cole Miller, Guido Mueller, Germano Nardini, Gijs Nelemans, Miquel Nofrarias,Antoine Petiteau, Paolo Pivato, Eric Plagnol, Ed Porter, Jens Reiche, David Robertson, NornaRobertson, Elena Rossi, Giuliana Russano, Bernard Schutz, Alberto Sesana, David Shoemaker,Jacob Slutsky, Carlos F. Sopuerta, Tim Sumner, Nicola Tamanini, Ira Thorpe, Michael Troebs,Michele Vallisneri, Alberto Vecchio, Daniele Vetrugno, Stefano Vitale, Marta Volonteri,Gudrun Wanner, Harry Ward, Peter Wass, William Weber, John Ziemer, and Peter Zweifel.Laser Interferometer Space Antenna. arXiv e-prints , page arXiv:1702.00786, February 2017.[39] Tristan L. Smith and Robert Caldwell. LISA for Cosmologists: Calculating the Signal-to-NoiseRatio for Stochastic and Deterministic Sources.
Phys. Rev. D , 100(10):104055, 2019.[40] B.P. Abbott et al. Search for the isotropic stochastic background using data from AdvancedLIGO’s second observing run.
Phys. Rev. D , 100(6):061101, 2019.[41] Neil Cornish and Travis Robson. Galactic binary science with the new LISA design.
J. Phys.Conf. Ser. , 840(1):012024, 2017.[42] Travis Robson, Neil J. Cornish, and Chang Liu. The construction and use of LISA sensitivitycurves.
Class. Quant. Grav. , 36(10):105011, 2019.[43] Paolo Campeti, Eiichiro Komatsu, Davide Poletti, and Carlo Baccigalupi. Measuring thespectrum of primordial gravitational waves with CMB, PTA and Laser Interferometers. 7 2020.[44] Matthew R. Adams and Neil J. Cornish. Detecting a Stochastic Gravitational WaveBackground in the presence of a Galactic Foreground and Instrument Noise.
Phys. Rev. D ,89(2):022001, 2014.,89(2):022001, 2014.