Probing the galactic and extragalactic gravitational wave backgrounds with space-based interferometers

We employ the formalism developed in [1] and [2] to study the prospect of detecting an anisotropic Stochastic Gravitational Wave Background (SGWB) with the Laser Interferometer Space Antenna (LISA) alone, and combined with the proposed space-based interferometer Taiji. Previous analyses have been performed in the frequency domain only. Here, we study the detectability of the individual coefficients of the expansion of the SGWB in spherical harmonics, by taking into account the specific motion of the satellites. This requires the use of time-dependent response functions, which we include in our analysis to obtain an optimal estimate of the anisotropic signal. We focus on two applications. Firstly, the reconstruction of the anisotropic galactic signal without assuming any prior knowledge of its spatial distribution. We find that both LISA and LISA with Taiji cannot put tight constraints on the harmonic coefficients for realistic models of the galactic SGWB. We then focus on the discrimination between a galactic signal of known morphology but unknown overall amplitude and an isotropic extragalactic SGWB component of astrophysical origin. In this case, we find that the two surveys can confirm, at a confidence level ≳ 3σ, the existence of both the galactic and extragalactic background if both have amplitudes as predicted in standard models. We also find that, in the LISA-only case, the analysis in the frequency domain (under the assumption of a time average of data taken homogeneously across the year) provides a nearly identical determination of the two amplitudes as compared to the optimal analysis.

detecting an anisotropic Stochastic Gravitational Wave Background (SGWB) with the Laser Interferometer Space Antenna (LISA) alone, and combined with the proposed space-based interferometer Taiji.Previous analyses have been performed in the frequency domain only.
Here, we study the detectability of the individual coefficients of the expansion of the SGWB in spherical harmonics, by taking into account the specific motion of the satellites.This requires the use of time-dependent response functions, which we include in our analysis to obtain an optimal estimate of the anisotropic signal.We focus on two applications.Firstly, the reconstruction of the anisotropic galactic signal without assuming any prior knowledge of its spatial distribution.We find that both LISA and LISA with Taiji cannot put tight constraints on the harmonic coefficients for realistic models of the galactic SGWB.We then focus on the discrimination between a galactic signal of known morphology but unknown overall amplitude and an isotropic extragalactic SGWB component of astrophysical origin.In this case, we find that the two surveys can confirm, at a confidence level ≳ 3σ, the existence of both the galactic and extragalactic background if both have amplitudes as predicted in standard models.We also find that, in the LISA-only case, the analysis in the frequency domain (under the assumption of a time average of data taken homogeneously across the year) provides a nearly identical determination of the two amplitudes as compared to the optimal analysis.

Introduction
The first direct observation of Gravitational Waves (GWs), generated by the merging of compact objects in 2015 was achieved by the LIGO-Virgo collaboration [3].The observation gave rise to one of the most rapidly expanding fields in contemporary physics.The technology implemented by the collaborations using interferometers for gravitational wave detection (currently operational ones include LIGO [4], Virgo [5], and KAGRA [6]) is continuously improving.The increase in sensitivity and integration time is leading to improving upper limits on the amplitude of several models of the SGWB [7,8], with a detection of an astrophysical SGWB being a realistic possibility over the next few years.
Among the experiments that will collect GW data in the next decade, the space-based LISA mission [9] is undoubtedly one of the most promising.Consisting of three satellites orbiting the Sun in elliptical orbits while maintaining a (nearly) constant distance of 2.5 million km from each other, the instrument will be capable of observing gravitational waves in a frequency range that is inaccessible to ground-based interferometers.LISA is currently in the definition phase, where the mission system requirements are being reviewed.This phase will end with possible adoption in 2024, which would see the experiment fly in ∼2035 and begin data collection a few years later [9,10].
On the other hand, Taiji [11,12], an analogous experiment proposed by the Chinese Academy of Sciences and currently in the planning stages, may be operative by the time that LISA will be flying.The proposed Taiji experiment has slightly longer arms than LISA (3 million km) and a sensitivity that peaks in the same frequency range.
At these frequencies, we expect many physical phenomena that can produce a GW signal to be active.In this work, we focus on the class of sources that cannot be resolved as individual events, but which can produce a globally detectable effect given their large number.The signal coming from the superposition of all the individually undetectable sources of GW radiation is called the Stochastic Gravitational Wave Background (SGWB).
While only upper bounds exist at frequencies probed by ground-based interferometers [7,13], interest in SGWBs has increased significantly over the past year after the recent claims by Pulsar Timing Array (PTA) collaborations [14][15][16] of a detection of a GW signal at much lower frequencies than those probed by LISA.These measurements may be evidence of an SGWB at frequencies ∼ 10 −9 − 10 −8 Hz.Analogously, we expect a SGWB to exist at LISA and Taiji frequencies.These would be sourced from several different phenomena.On the one hand, there is an astrophysical component of the SGWB, mainly constituted by the GW signal coming from mergers of compact binary systems, such as stellar mass black holes and neutron stars.On top of this signal, there are GWs generated by other astrophysical sources, such as massive black holes (that are expected from several models explaining the formation of Supermassive Black Holes (SMBH)), other galactic binaries, and Extreme Mass Ratio Inspirals (EMRIs).[17].On the other hand, in the LISA frequency band, a cosmological GW background may also be present [18].In this case, the mechanisms capable of sourcing a GW signal can have different origins.There is the possibility to observe a background that originates from inflation, preheating, phase transitions, topological defects, primordial black holes (see [19] for a general overview), or other (known or as yet unknown) early universe phenomena.
The variety of phenomena mentioned above, together with the fact that we can only access their superposition in the background, results in many challenges from the point of view of data analysis.However, the distinct shapes in frequency and angular distribution of the Power Spectral Density (PSD) of the various signals can be exploited to better model the effect and to optimize the analysis process.In particular, thanks to the directional dependence, we may be able to distinguish between galactic and extragalactic components and search for any correlation with known tracers of structure [20][21][22][23].This paper aims to exploit the data analysis techniques described in [1] within the framework of groundbased detectors and to apply them to the LISA data analysis to produce the final forecast on the observability of some features of the SGWB.We focus on detecting the anisotropies of an SGWB [24,25], taking advantage of the periodic motion of LISA, and then of LISA in combination with Taiji in their trajectory around the Sun, along with a complete timefrequency analysis.A network of constellations of space-based interferometers has also been considered in [26]; in contrast to that work, we leverage our analysis on the specific trajectory of the satellites, with the application to the LISA -Taiji pair.
Further information on the origin of the SGWB can be obtained by studying its polarization [27][28][29][30][31][32].A coupling between a pseudo-scalar inflaton and a gauge field, can result in a fully circularly polarized SGWB [33].Still, a degree of net polarization might also appear in the astrophysical SGWB due to Poisson fluctuations in the (finite) number of unresolved sources [34] and might be potentially detectable for a network of space based interferometers [35].As explained below, in our treatment we consider an intrinsically stationary background, i.e. whose statistical proprieties do not change with time.This assumption is reasonable for cosmological and astrophysical backgrounds made of a superposition of many unresolved sources, as predicted in most theoretical models.However, analysis can be extended to non-stationary cases [36,37].
The plan of the paper is as follows.In section 2, we summarize the most useful definitions in setting up the analysis of SGWBs and introduce the theoretical parameters we would like to probe.This is done following the formalism of [1] and [2].In section 3 we go beyond what was developed in [1] to analyze the case of LISA and Taiji, showing how one can build physical observables that can be used to probe the theoretical features of a SGWB.We show how to perform the optimal analysis, taking into account the nuisance effects due to instrumental noise and the variance of the signal as we did for ground-based detectors in [1], defining an optimally filtered chi-squared in the measurement.We then produce the forecasts of two analyses that can be performed using this formalism.In section 4 we show how LISA alone and LISA-Taiji combined can constrain the anisotropic coefficients of the map of the astrophysical background of our galactic binaries [38,39], along with the total energy density.On the other hand, in section 5 we show how accurately the two networks of instruments can disentangle a galactic and an extragalactic astrophysical component of the SGWB [40], taking advantage of the different frequency dependencies of the power spectra and the different angular distribution of the sources.In the last section 6 we summarize the results of our analysis and discuss the advantages and limitations of the two experiments in both applications.

The stochastic gravitational wave background
The starting point of our analysis is the formalism that we developed in [1,41] to determine the observability of anisotropic modes of SGWBs with generalized networks of baselines.We refer the interested readers to these works for details.Generally, we consider gravitational waves at x α , the location of the detectors, as small perturbations of the metric tensor on a Minkowski background in such a way that it can be decomposed in plane waves (the usual h 00 = h 0a = h aa = ∂ a h ab = 0 gauge has been adopted, where Latin indices run only over the spatial components) where n is the unit vector on the 2-sphere, and e R ab (n), e L ab (n) are the two chiral polarization basis tensors, with the same normalization adopted in [2], i.e. e λ * ab e λ ′ ab = δ λλ ′ .The direction vector n is defined as n = cos ϕ sin θ x + sin ϕ sin θ ŷ + cos θ ẑ.
We define (2.1) in ecliptic reference coordinates, setting the ẑ axis to point along the direction perpendicular to the plane of the Earth orbit.To describe a stochastic source, we treat the complex amplitude h λ (f, n) as a random Gaussian variable with zero mean.This means that its statistics are completely specified by its variance.The dependence of the stochastic background on frequency and direction may be stated solely in terms of the expectation value of the two-point correlator for the random variable h λ (f, n) as where δ (2) D (n − n′ ) is a covariant two-dimensional Dirac delta-function on the unit two-sphere and δ D (f − f ′ ) the Dirac delta function on the frequency space.In some cases, we suppose that we have a factorized dependence of the angles and the frequencies in the power spectrum, such as where the function Ĩ00 (|f |) is introduced analogously in [2].Although commonly adopted in the literature, this factorization is not generic, and, in fact, in section 5 we consider a more complex angular and frequency dependence.It is conventional to define the fractional energy density in GW per logarithmic frequency interval and per polarisation as where ρ crit is the current energy density of the universe, which we assume to be flat, while H 0 is the Hubble constant.In the second step, we have used the fact that only the monopole contributes to the total GW energy density (after integrating over the arrival direction).
Next, we decompose the angular power spectrum in terms of spherical harmonics, writing it as where the value of δ GW 00 is set considering the convention for the monopole Y 00 = 1 √ 4π .Finally, we rewrite the frequency dependence of Ω GW (f ) as so that Ω 0 is the fractional energy density at the pivot frequency f 0 .In order to produce numerical forecasts in our work, we suppose that we have a perfect knowledge of this function.However, it is in principle possible to write it in parametric form and infer the value of the parameters in the same analysis.

SGWB effect in space-based detectors
We consider two possible configurations of space-based interferometers.We call LISA-self the network made of three satellites in a triangular-shaped configuration that constitutes the LISA experiment as in figure 1.On the other hand, we call LISA-Taiji the network made of the 3 LISA satellites along with Taiji, a proposed experiment by the Chinese Academy of Sciences [11,42] that would be built similarly to LISA, having 3 satellites in a triangularshaped configuration as in figure 1.We denote by ∆T L ij (t) the modification to the time of flight of a photon in the laser beam between the i-th and the j-th satellite of LISA due to the effect of a gravitational wave, where i, j = {1, 2, 3}.Analogously, ∆T T ij (t) will be the same time of flight distortion between two satellites in the Taiji constellation.One can show that where lL ij (t) is the unit vector going from ⃗ x L i to ⃗ x L j , the i-th and j-th satellite of LISA.An analogous expression holds for Taiji, with the suffix L replaced by the suffix T in all these quantities.In our numerical evaluations, the unperturbed distance L L between two LISA satellites is set equal to 2.5 × 10 6 km, while for Taiji the unperturbed distance L T is 3 × 10 6 km [42].The actual trajectory considered for the six satellites is explained in full detail in appendix A.
Following the formalism adopted in [2], we introduce the set of three Time-Delayed Interferometric (TDI) measurements ∆F L i (t) made by a combination of the Doppler frequency shifts between the LISA satellites (∆F T i (t) will be the same combinations for the Taiji constellation).We reference (4.11) in [2] for the precise definition of these quantities.Under the assumption of equal arm lengths (which means that both LISA and Taiji are perfect equilateral triangles; in appendix A we verify the accuracy of this assumption), we can combine the TDI data streams ∆F L/T i (t) in order to diagonalize the noise matrix.This accounts for performing the usual linear transformation in the "AET" basis After some straightforward algebra (see [2] and [43] for more details) the self-correlators between each pair of the LISA TDI data streams in the "AET" basis read with f On the other hand, the cross-correlators between the "AET" channels of LISA and the "AET" channels of Taiji are where the time dependence of ⃗ x L/T i (t) and lL/T ij (t) is understood and where the explicit expression of the functions W L/T (f ) and R λ (f n, lij , lik ) can be found in eqs.(4.14) and (A.4) of [2], respectively.As stated above, the actual measurement in the detector can be affected by some instrumental noise in such a way that where we assume that the noise ñL/T O (f, t) is drawn from a Gaussian distribution and that noises of different detectors are uncorrelated, i.e.
The functions N L/T O (f ) are given in eqs.(B14) and (B15) in [2]. 1 We can now build the final estimator from the self-correlators of LISA TDI datastreams: where QLL OO ′ (t, f ) are filter functions to be optimized to maximize the signal-to-noise ratio of the measurement.In a similar way, the estimator built by making use of the cross-correlators between LISA and Taiji datastreams is We stress a fundamental difference between the two analyses considered.In the first case (when only LISA data streams are available) a perfect knowledge of the noise curve is assumed: any error in the estimate of the noise will translate into a bias in the final estimator C LL [45].On the other hand, since it is reasonable to assume no correlation between LISA and Taiji realizations of the noise, in the second analysis an imprecise estimate of the noise within each constellation would not bias the expectation value of the estimator C LT .This is a well-known fact in analyses with ground-based detectors, where self-correlators of measurements of the same detector are typically not included in the analysis to avoid biases from the imperfect knowledge of the noise function N [46].For this reason, in the LISA-Taiji analyses we disregard any self-correlators of two measurements taken within one constellation.

Signal expectation value and chi-squared
After the same steps performed in [2], we find that the expectation value of the estimator (3.8) is where An analogous formula holds for ⟨C LL ⟩ and can be found as (4.36) in [2].Next, as done in [1], we build a chi-squared in the quantity where θ = (Ω 0 , {δ GW ℓm }) are the parameters of the theory that we want to probe and θ their fiducial values.Following the analogous steps in [1], we obtain the final expression of the optimal filtered chi-squared on the estimator for the LISA-Taiji analysis where and where We note that D OO ′ (f, t) is a real quantity.This follows from the fact that, assuming a real map for the anisotropic distribution of the galactic signal (as (B.2) is), we have In these expressions, the index O runs over the LISA channels and the index O ′ over the Taiji ones.The chi-squared evaluates to and we note that, as in the analysis of [2], the TDI factors W L,T (f ) cancel out from this expression.B14) and (B15) in [2].As discussed in the text, for LISA we have chosen a value for the parameters P = 15, A = 3, and L = 2.5 × 10 6 km, while the PSD of Taiji is computed assuming P = 8, A = 3, and L = 3 × 10 6 km.equation of (3.15) is (4.38) in [2].We perform the integration in time of (3.15) numerically, by taking advantage of the fact that the response functions R ℓm OO ′ (f, t) are periodic in time with period T s = 1 yr.Therefore we replace the upper time of integration in (3.15) with T s and we multiply the results by the number of years of observation n.Furthermore, as we are unable to perform the time integration analytically, we discretize it in a large but finite amount N t of time intervals.The numerical evaluation of the time-dependent response function is extremely expensive, and so we evaluate it once in each time interval, at the times t i = i Ts Nt with i ∈ {1, . . ., N t }.The quantity N t is limited from below by the requirement that the observed sky does not change appreciably during each time interval.More precisely, as we assume stationary signal statistics, we require that the constellations vary their angular position (in their motion around the sun) by an angle ≪ 180 • lmax , where l max is the maximum multipole that we include in our analysis (choosing smaller values of N t would smear the measurement of the higher multipoles).For l max = 30, the limiting angle is 6 • , which is spanned in about 6 days.Correspondingly, we must choose N t ≫ 60.The quantity N t is also limited from above, by the requirement that all frequencies of our interest are well probed within each time window.The smallest frequencies to which the measurement is sensitive is f min = O 10 −4 Hz .This translates into Ts Nt f min ≫ 1, which gives N t ≪ 3000.In our analysis, we choose N t = 365.With the discretization in time, (3.12) is cast in the form (3.16) In the following sections, we employ this procedure for two different analyses.Firstly, in section 4 we forecast the measurement of the energy density Ω 2 0 and of the anisotropic coefficients δ GW ℓm for a SGWB generated by the galactic sources only.As a second example, in section 5 we consider a SGWB resulting from the superposition of a galactic and an extragalactic signal, forecasting the observability of the energy density of the two components.

Reconstructing the galaxy map
As a first example, in this section we would like to show how LISA and LISA with Taiji can reconstruct the anisotropic distribution of the galactic background, assuming that this is the only component of the SGWB.To do so, we first assume that the frequency dependence of the power spectrum is [39] which is the best fit for the SGWB after the removal of the resolvable sources. 3In agreement with [38], we have a fiducial value of Ω0 = 2 × 10 −11 .As a second assumption, the model of the galaxy used to set the fiducial value of the anisotropic coefficients δGW ℓm is explained in appendix B.
Following the same procedure adopted in [1], we perform a Fisher forecast on the measurement of the GW amplitude and anisotropies, by expanding the χ 2 to quadratic order in the departure between their measured and their fiducial values.We obtain where the Fisher matrix elements are and where Ω0 and δGW ℓm are the fiducial values of the corresponding non-hatted quantities.The Fisher matrix elements turn out to be where we defined and where for the numerical evaluations we chose H 0 = 70 km s −1 Mpc −1 = 2.27 × 10 −18 Hz.From the chi-squared in (4.2) we can then define the posterior distribution of the energy density Ω 0 and the anisotropic coefficients {δ ℓm } as where π (Ω 0 , {δ ℓm }) is the prior distribution on our parameters, which we assume to be flat, and N the probability normalization factor.In figure 3 we show the forecast errors in the measurement of each coefficient δ GW ℓm after marginalizing the posterior distribution (4.6) over all the other parameters, in a similar way to what has been done in [1].Not surprisingly, a network made of LISA and Taiji would give better constraints in the measurement of the high ℓ anisotropic coefficients of the SGWB.Because of the longer baseline, the overlap functions of the instruments for each anisotropy up to ℓ = 20 have significant support at the frequencies where the LISA and Taiji will operate.On the other hand, for the lowest even multipoles (ℓ = 2, 4) LISA-alone shows a better response. 4This is due to the fact that in the opposite limit, i.e.where the baseline between the instruments is short, the response to these specific multipoles is not suppressed [41].

Disentangling the galactic and the extragalactic signals
In the previous section, we studied the detectability of the anisotropies of the SGWB without any prior knowledge of them.More precisely, the galactic morphology outlined in appendix B has been taken as a fiducial model, but its knowledge has not been assumed in the data analysis.In this section, we instead assume the specific anisotropic galaxy pattern determined by the galaxy, with the specific spatial distribution described in appendix B, both in the fiducial model and as a prior knowledge in the analysis.This means that, in the present analysis, the coefficients δ GW ℓm in (3.12) have known values, and we assume that only the overall scale Ω G of the galactic signal is unknown.We also assume the presence of an isotropic extragalactic component of unknown amplitude Ω B , so that the full signal is where ψ G (f ) is the same functional frequency dependence adopted in the previous section, and, as already mentioned the coefficients δ GW ℓm are obtained from the spatial distribution given in appendix B. Finally, the frequency dependence of the extragalactic signal is the expected one, as also considered in [40].With the assumption of perfect knowledge of the galaxy morphology, and of the frequency-dependence of the galactic and extragalactic components, the two only unknown quantities in the signal (5.1) are the two amplitudes Ω B and , while the second is obtained with ψ(f ) = ψ G (f ) as in (4.1).A total observation time of T = 10 yr and a fiducial value of Ω0 = 2 × 10 −11 at the pivot frequency f 0 = 1 mHz are assumed.Ω G , which provide the fractional energy density contributions of the two components at the pivot scale f 0 = 1 mHz5 .The chi-squared introduced in (3.12), then becomes a function of these two parameters: where the denominator D OO ′ is formally defined as in (3.13), in terms of the new function (5.3) Accordingly to [38] and [40], we set the fiducial values of the unknown parameters to ΩG = 2 × 10 −11 and ΩB = 3.78 × 10 −13 , both defined at the pivot frequency f 0 = 1 mHz.We can further simplify the notation to by introducing the dimensionless response functions (divided by the total variance) for the extragalactic, and for the galactic signal.The Fisher approximation of the chi-squared (5.4) is where the Fisher matrix elements are (5.8)

Frequency-only vs. time-frequency analysis
The filter functions Q adopted in the data analysis are weights that combine measurements obtained at different frequencies to maximize the signal-to-noise ratio (SNR) of the combination.Under the hypothesis of stationary signal and noise, and for the case of an isotropic signal, the optimal weights Q are time-independent, as the same linear combination of the measured frequencies maximizes the SNR at all times.On the contrary, the weights QLL OO ′ (f, t) and QLT OO ′ (f, t) employed in our analysis, see eqs.(3.7) and (3.8), depend also on time, accounting for the fact that the time positions of the satellites change with time, so that different linear combinations of the measured frequencies maximizes at different times the SNR associated to an anisotropic signal (that is measured differently at different times even for a stationary signal, due to the time-dependence of the satellite response functions).
A simpler analysis can be performed by ignoring this time information and performing a frequency-only analysis, in which essentially data obtained at any given frequency are averaged over the year.This will result in a sub-optimal estimator.We want to assess the improvement of the optimal time-frequency analysis done above vs. this simplified frequencyonly analysis.
The frequency-only analysis is performed by integrating along the whole time of observation the cross-correlators between the different channels (and instruments in the case of the LISA-Taiji network) first, and by then applying the filters QLL OO ′ (f ) and QLT OO ′ (f ), which are now functions of the frequency only.Concretely, this means that eqs.(3.7) and (3.8) are replaced by (5.9) We can now proceed along the lines of the analysis performed in section 5 by introducing the time-averaged detector overlap functions and by employing them in the expressions The optimal chi-squared of this frequency-only analysis (which is sub-optimal when compared to the one computed in the previous section) acquires the form where the coefficients associated with the extragalactic and the galactic term are, respectively, and where the denominators D OO ′ (f ) are related to the functions T O ′ O (f ) analogously to (3.13), and the Fisher matrix elements are in this case (5.14) We stress that the chi-squared (5.12) differs from (5.4) from the fact that here we take an average over the time of the response functions first, while in the previous analysis, we built a full measurement for each time window we considered.As we already remarked, and as we show with one explicit example in figure 4, this simplified method presented here results in a less accurate determination of the model parameters.

The posterior distribution and final forecast
In a similar way to what we did in the previous section, we introduce the posterior probability of the two parameters of our problem, Ω G and Ω B : where also in this case a flat prior π(Ω G , Ω B ) is assumed and where the chi-squared is defined in (5.4) or (5.12) depending if we want to adopt a time-frequency or a frequency-only analysis.
In figure 4 we show the contour levels of (5.15) for a total observation time of T = 10 years with either the LISA constellation alone or the LISA-Taiji networks (in this second case, we do not include correlations between two measurements performed within satellites of the same constellation).We note from the figure that, at face value, LISA alone would perform better than when combined with Taiji.However, we stress that any SGWB search performed by LISA alone provides a perfectly unbiased estimator of the signal (in (3.7))only if a perfect knowledge of the noise statistics is assumed.On the other hand, the cross-correlators between LISA and Taiji TDIs as in (3.8) are not noise-correlated.We furthermore note that a time-frequency analysis would not make significant improvements with respect to a frequency-only analysis for LISA alone (for LISA-self, the contours of the frequency-only and of the time-frequency analysis in Figure 4 are practically superimposed to each other, and cannot be distinguished by eye), while there is a significant improvement for a LISA-Taiji joint analysis.The reason is the fact that we assumed an isotropic extragalactic SGWB in (5.1), for which the response functions of LISA do not change with time.The same is not true for the LISA-Taiji analysis, due to the variation of the relative orientation of LISA and Taiji arms with respect to their baseline.This also increases the time variation of the LISA-Taiji response function to the galactic signal, with respect to the LISA-self case.In appendix C we verify that indeed, for the LISA-self case, the time-dependence of the galactic signal is not strong enough to make the time-frequency analysis significantly better than the sub-optimal frequency-only one.The red regions, with thick lines, correspond to the forecast posterior distribution of a measurement performed by the LISA constellation alone.The blue regions, with dashed lines, represent the forecast obtained using a network made of LISA and Taiji, without including self-correlations within each constellation.This second analysis provides a less precise determination of the model parameters, but it is not biased by imperfect knowledge of the instrumental noise.For comparison, a frequency-only analysis for the LISA-Taiji network leads to the contour levels displayed with dotted lines.For LISA-self, the frequency-only analysis produces contours superimposed to the solid ones.A time of observation of T = 10 yr and fiducial value of ΩG = 2 × 10 −11 , ΩB = 3.78 × 10 −13 are assumed, as expected from [38] and [40] (as discussed in the text).

Conclusions
In the present work, we applied the formalism developed in [1] and [2] to produce forecasts on the measurement of a stochastic GW background with two possible configurations of spacebased detectors.We considered two possible analyses, under the assumption of a stationary SGWB in its reference frame, for both a network made of the three LISA satellites and a second network where LISA and Taiji can operate jointly.In the first analysis, we produced a fiducial model of the distribution of the GW galactic signal and we took advantage of their anisotropic distribution in the reference frame of the Earth to forecast at which level of accuracy the sky maps of these sources can be obtained.We realized that, because of the fact that the GW emission by the unresolved galactic binaries interests just the lowest frequencies detectable by LISA and Taiji (as in (4.1)),only the very first multipoles {δ GW ℓm } of the anisotropic distribution of the sources are resolvable.This is due to the fact that the higher multipoles overlap functions of the instruments, as in (3.10) have their maximum at higher values of frequencies, where there are no unresolved galactic sources emitting.In this case, the network of LISA and Taiji outperforms LISA alone thanks to the fact that the baselines between the LISA and the Taiji satellites are greater than the distance between each pair of LISA satellites.This is clear from figure 3, where the forecast error σ δ GW ℓm represents the accuracy in determining the value of the anisotropic coefficient δ GW ℓm as introduced in (2.5).In general, the network made of LISA and Taiji displays a better sensitivity to the anisotropic coefficients up to ℓ = 20 with respect to LISA alone, because of the longer baseline between the instruments.On the other hand, LISA alone has an enhanced response to the monopole and the first even ℓ value multipoles, since the corresponding overlap functions are not suppressed in the low-baseline limit.Since the coefficients δ GW ℓm are of the order of O(10 −1 ) or smaller, it means that a high-resolution galaxy map of the galactic background with our space-based detectors is difficult to obtain.These forecast errors decrease with the square root of the total time of observation and with the squared of the power spectral density of the noise (since the data stream is almost noise-dominated).
From the results displayed in figure 3 we see that it is crucial to adopt the real model of the galaxy power spectrum as a function of the frequency, i.e. after the removal of the resolvable sources.This removal subtracts power from the SGWB at comparatively higher frequencies, thus worsening the detectability of the anisotropic dependence of the residual SGWB.While the lower panel of figure 3 is a better representation of the detectability of the SGWB of galactic origin, the results in the upper panel can be used as a proxy on the detectability of an extra-galactic component, for which individual high-frequency sources cannot be removed. 6nother novelty of this analysis is the fact that, for the first time, we have a precise forecast of the measurement of the anisotropic coefficients δ GW ℓm separately.For example, the study of the anisotropy presented in [2] does not provide a precise forecast for the measurement of the specific coefficients δ GW ℓm , but it estimates the level of their detectability under the assumption of statistical isotropy.Under this assumption, ref. [2] provided an estimate of the level of detection of the combination C GW ℓ ∝ m δ GW ℓm 2 .As shown in [2], the LISA response function for this combination is isotropic, and therefore the motion of the LISA satellites is completely irrelevant to that assessment.We, on the contrary, have studied the determination of the specific anisotropic coefficients δ GW ℓm , for which the precise time-dependent orientation of the LISA constellation must be known.The same is true when Taiji is added to the analysis.For this reason, we adopted the model for the satellite trajectories presented in appendix A, and we then computed the overlap functions of the LISA and Taiji networks as a function of time.
As a second analysis, we treated two possible sources of SGWB radiation together and we forecasted the ability of LISA and LISA+Taiji to disentangle their separate contributions.We considered a galactic background as in the previous analysis.However, in this case, we assumed perfect prior knowledge of the frequency dependence of the power spectrum and its anisotropic distribution in the data analysis, keeping as the only unknown the overall logarithmic energy density, i.e. the parameter Ω G in (5.1).On top of the galactic signal, we assumed the presence of an isotropic extragalactic component, of astrophysical nature, due to the superposition of gravitational waves sourced by a large number of compact binaries [40].We assumed the typical f 2/3 frequency dependence of this signal, keeping also in this case, a free normalization Ω B in (5.1).In figure 4 we show the posterior distribution of the parameters Ω G and Ω G , centered in the value predicted by our fiducial model.The result appears to indicate that LISA alone performs better than when correlated with Taiji.This is due to the fact that the main contributions to the data stream of both the sources of GW background radiation come from the lowest multipoles, at which LISA alone is more sensitive than when correlated with Taiji.However, we stress the fact that a measurement performed with LISA alone is perfectly unbiased only under the ideal condition of perfect knowledge of the noise of the instrument: any error in the estimate of the noise budget would translate into a bias in the final estimate.On the other hand, the cross-correlation of LISA and Taiji produces an unbiased estimator of the GW power spectrum, as the noises of the two constellations are not cross-correlated.
A final important point we would like to stress is the advantage of performing a timefrequency analysis compared to the traditional frequency domain in the case of a network made of multiple instruments.As it appears evident from the figure 4, a network made of LISA and Taiji will produce stronger constraints when time-frequency response functions are employed (the thick lines of the plot) with respect to the frequency-only analysis (the dashed lines of the plot).This is mainly due to the relative orientations of the baseline and the arms of the two instruments, which results in a very complex time dependence of the response (overlap) functions of the network.On the other hand, when LISA alone is considered, the response (overlap) functions have a dominant constant plus a subdominant time-dependent component which is very regular over time (see appendix C for more details).Assuming that data are continuously taken over the year (or, more realistically, that gaps in the data are evenly distributed across the orbit), one finds that the analysis in the frequency-only domain (with time-averaged response functions) gives only marginally worse results than the optimal time-frequency analysis.
where ω = 2π 1 year is the angular frequency of the orbit, and where To describe the position of the satellites of a given constellation (LISA or Taiji) in this coordinate system, we introduce a second rotation matrix and the tree unit vectors In terms of these quantities, the position of the three satellites is given by where the Greek index λ = 1, 2, 3 refers to one of the three satellites in the constellation, while the Latin indices (which also run from 1 to 3) are coordinate indices.The two constant phases α c and β c are specific to each constellation.The expression (A.5) can be written more explicitly in the form The implicit relation (A.5) manifestly shows that the motion of each constellation can be thought of as the composition of two rotations, both of periods equal to one year.The rotation R I [ω t + α c ] keeps the constellation at fix relative phase α c with respect to the Earth. 7Simultaneously, the three satellites rotate with R II [ω t + β c ], which is a rigid frequencies is expected to be that of unresolved galactic binaries [38,39].We build a template morphology for the gravitational wave signal of galactic binaries using a bulge and disc stellar model of the galaxy [47,48].The stellar density, as a function of galactocentric cartesian coordinates, is given by where We form a template of the gravitational-wave sky by integrating the density in logarithmic intervals along all lines-of-sight l p corresponding to each Healpix pixel p on a map of resolution N side = 32 with pixels of angular size 1.8 • and angular Nyquist scale ℓ ≈ 120.
where the origin of each line-of-sight is set to the position of the sun in the galactocentric coordinate system, [−8.1, 0, 0] kpc, and it is extended to 20 kpc in all directions around the sun location.The resulting map is shown in Figure 5.
Given the Fisher matrix estimates of (4.4) for the relative anisotropy δ GW ℓm and the spectral coefficients of the unit normalized template a M ℓm , we can produce realizations of the 'observed' sky that account for noise levels and mode coupling induced by each network of baselines.To generate the realization of the noise with the correct correlation structure, we generate a set of random, uncorrelated normal variates ξ ℓm with unit variance and rotate to a set a N lm using the hermitian square root of the inverse Fisher matrix C Time-frequency vs. frequency only analysis for LISA-self Figure 4 indicates that, for the analysis performed by the LISA constellation alone, the timefrequency domain does not bring a significant advantage with respect to the frequency-only analysis.This Appendix discusses the reason for this.Firstly, let us consider the signal from the galaxy, which, due to its angular distribution displayed in figure 5, induces a timefrequency response function (5.6) that is well fit by where G OO ′ (f ) is a function of the frequency only and A is the relative amplitude of the modulation of the galactic signal in the orbit of LISA as a function of time.Figure 6 shows the goodness of the fit for two frequencies in the A channel.We note that the modulation is the same for all frequencies (namely, A is constant, and no frequency-dependent phase is present in the argument of the cosine) as a consequence of the factorization of the angular  and frequency dependencies in the galactic signal in eq.(5.1).From our numerical fit, we obtain A ≃ 0.185, namely the amplitude of the modulated term accounts for about 20% of the overall amplitude.The frequency-only galactic response function, as introduced in the second line of eq.(5.13) for the sub-optimal analysis then evaluates to RG OO ′ (f ) ≡ namely, the term that encodes the modulation cancels out in the average (this assumes measurements that are taken homogeneously across the yearly revolution).Secondly, we note that the response functions to the isotropic extragalactic background to be used in the time-frequency and in the sub-optimal frequency-only analysis (see eq. ( 5.5) and the first line of eq.(5.13)) coincide, namely The final forecasts displayed in Figure 4 depend on the Fisher matrices, introduced in eq.(5.8) for the time-frequency analysis and in eq.(5.14) for sub-optimal frequency-only analysis, respectively.From the expressions written in this Appendix, the elements of the Fisher matrix for the time-frequency analysis acquire the form (C.4) From the above discussion, we then see that the elements of the Fisher matrix for the suboptimal frequency-only analysis are also given by eqs.(C.4), if one simply sets A = 0 in the F Ω G Ω G element.From the value of A given above, we see that this term induces an increase of only A 2 /2 ∼ 1.7% in this element in the time-frequency vs. the frequency-only analysis.We now employ these Fisher matrix elements to plot the contour levels of the two analyses, and we show that we can reproduce the near coincidence of their results observed in figure 4 for the LISA-self case.Specifically, let us consider the 3σ contours (although identical considerations can be done for all the other confidence levels), that we show in figure 7, to appreciate the difference between the two analyses.We note that the optimal analysis results in a slightly better determination of the parameter (smaller area of the 3σ region), which was not noticeable by eye with the scale shown in figure 4. From the Fisher matrix elements, the normalized posterior distribution of the two parameters Ω B and Ω G , assuming a flat prior, is  factor, leading to a discrepancy of ΩG obtained in the two case quantified as As we already mentioned, A = 0.185.We numerically find

Figure 1 :
Figure 1: Scheme of the LISA-Taiji configuration.The three LISA satellites are separated by an arm length of 2.5 × 10 6 km, while for Taiji the arm length is 3 × 10 6 km.The center of mass of each of the two constellations lays at 20 • from the position of the Earth in the ecliptic plane, and it follows the trajectory of the Earth during its rotation around the Sun.The precise trajectory of the six satellites is described in appendix A. We refer as LISA-self configuration the same in the figure where however the three Taiji satellites are not present.

Figure 2 :
Figure 2: Freuquency dependence of the noise PSD Ñ L/T O (f ) as defined in eqs.(B14) and (B15) in[2].As discussed in the text, for LISA we have chosen a value for the parameters P = 15, A = 3, and L = 2.5 × 10 6 km, while the PSD of Taiji is computed assuming P = 8, A = 3, and L = 3 × 10 6 km.

Figure 3 :
Figure 3: The expected (marginalized) forecast error of each relative anisotropic coefficient, δ GW ℓm , as observed by a network made of the three LISA satellites and the LISA-Taiji constellations.The first panel assumes the frequency dependence of the signal (cf.(2.6))

Figure 4 :
Figure4: Contour lines at 1σ, 2σ and 3σ levels of the posterior distribution of(5.15).The red regions, with thick lines, correspond to the forecast posterior distribution of a measurement performed by the LISA constellation alone.The blue regions, with dashed lines, represent the forecast obtained using a network made of LISA and Taiji, without including self-correlations within each constellation.This second analysis provides a less precise determination of the model parameters, but it is not biased by imperfect knowledge of the instrumental noise.For comparison, a frequency-only analysis for the LISA-Taiji network leads to the contour levels displayed with dotted lines.For LISA-self, the frequency-only analysis produces contours superimposed to the solid ones.A time of observation of T = 10 yr and fiducial value of ΩG = 2 × 10 −11 , ΩB = 3.78 × 10 −13 are assumed, as expected from[38] and[40] (as discussed in the text).
R b is the characteristic radius of the bulge, and R d and Z d are the characteristic radius and scale height of the disc respectively.A gives the relative weight of the stellar density in the bulge compared to the disc.The normalization of the stellar density ρ 0 is set to unity as we are only interested in the morphology of the signal.We set A = 0.25, R b = 500 pc, R d = 2500 pc, and Z d = 200 pc.

Figure 5 :
Figure 5: Mollweide projection of the template morphology for the galactic GW signal arising from unresolved binaries.The template is unit normalized and in celestial coordinates.

Figure 6 :
Figure 6: Time evolution of the response function RG AA (f, t) during one year, for two specific frequencies chosen as an example.Only the real part is shown, as the imaginary part is subdominant.The points are the numerical values used in our analysis (taken with a time interval of one day), while the solid line is the fitting function (C.1), with A = 0.185 in both cases.

Figure 7 :
Figure 7: A detail of figure 4, to visually appreciate the difference between the 3σ contours in the LISA-self analyses in the time-frequency domain (solid line) and in the sub-optimal frequency only case (dotted line).The central value of Ω B on the horizontal axis corresponds to the fiducial value ΩB .

2
On the other hand, for what concerns the LISA-self analysis, the corresponding