# Prospects for detecting anisotropies and polarization of the stochastic
gravitational wave background with ground-based detectors

G. Metasti,11footnotetext: Corresponding author. C. R. Contaldi  and M. Peloso

###### Abstract

We build an analytical framework to study the observability of anisotropies
and a net chiral polarization of the Stochastic Gravitational Wave Background
(SGWB) with a generic network of ground-based detectors. We apply this
formalism to perform a Fisher forecast of the performance of a network
consisting of the current interferometers (LIGO, Virgo and KAGRA) and planned
third-generation ones, such as the Einstein Telescope and Cosmic Explorer. Our
results yield limits on the observability of anisotropic modes, spanning
across noise- and signal-dominated regimes. We find that if the isotropic
component of the SGWB has an amplitude close to the current limit, third-
generation interferometers with an observation time of $10$ years can measure
multipoles (in a spherical harmonic expansion) up to $\ell=8$ with ${\cal
O}\left(10^{-3}-10^{-2}\right)$ accuracy relative to the isotropic component,
and an ${\cal O}\left(10^{-3}\right)$ amount of net polarization. For weaker
signals, the accuracy worsens as roughly the inverse of the SGWB amplitude.

###### Contents

  1. 1 Introduction
  2. 2 Main results
    1. 2.1 Results: forecast observations of anisotropic modes
    2. 2.2 Results: forecast observations of net polarization
    3. 2.3 Discussion
  3. 3 The stochastic gravitational wave background
  4. 4 Data stream of ground-based interferometers
  5. 5 Measurement by a network of detectors
    1. 5.1 Expectation value of the signal
    2. 5.2 Covariance matrix and optimal $\chi^{2}$
  6. 6 Detection of anisotropies
    1. 6.1 Expectation value of the signal
    2. 6.2 Covariance matrix and optimal $\chi^{2}$
    3. 6.3 Fisher forecast
  7. 7 Detection of a net polarization
    1. 7.1 Expectation value of the signal
    2. 7.2 Covariance matrix and optimal $\chi^{2}$
    3. 7.3 Fisher forecast
  8. A Signal expectation value
  9. B Covariance matrix
    1. B.1 The simplest network: two L-shaped detectors
    2. B.2 The most general network
  10. C Evaluation of the response functions
    1. C.1 Monopole response function for a chiral SGWB
  11. D Geometry of the detectors

##  1 Introduction

Following the first observation of the signal originating from the merging of
a massive, compact system [1], Gravitational Wave (GW) science has become one
of the most rapidly growing fields in fundamental physics. Since then, many
more detections have been made thanks to the constantly improving sensitivity
of the existing detector network that now includes LIGO, Virgo, and KAGRA [2].
The sensitivity of detectors is expected to keep improving in the coming
years. In the next ten to twenty years, a new generation of ground-based
interferometers, such as Einstein Telescope (ET) [3] and Cosmic Explorer (CE)
[4] will also become operational. These detectors will work alongside improved
versions of the current ones. The expectation is that these instruments will
be several orders of magnitude more sensitive to gravitational wave strain [3,
5] than existing ones. This sensitivity increase will significantly enrich our
catalogue of merger events of many different origins.

We still await the detection of a Stochastic Gravitational Wave Background
(SGWB) of any origin. Today, only upper limits exist as obtained from
analyzing LIGO–Virgo–KAGRA data. The energy density per logarithmic frequency
interval of a scale-invariant and an astrophysical-sourced SGWB are
constrained to be, respectively, $\Omega_{\rm GW}\lesssim 5\times 10^{-9}$ [2]
and $\Omega_{\rm GW}\lesssim 8.6\times 10^{-9}$ [6, 7, 8, 9].

There are several well-motivated theoretical reasons to expect that
astrophysical and cosmological phenomena can generate an SGWB. On the one
hand, at frequencies probed by ground-based detectors (1Hz-10kHz), an
astrophysical component due to the superposition of many unresolved sources,
primarily black holes and neutron star binary coalescences, is expected. On
the other hand, several mechanisms can produce a cosmological SGWB. Among them
are several non-minimal models of inflation, phase transitions, and
topological defects (see, e.g. [10, 11, 12, 13]).

Detection of an SGWB will generate significant interest in determining the
originating mechanism. Characterizing the SGWB will involve several efforts
and may be complicated by the presence of multiple underlying source
mechanisms. Firstly, to understand the nature of an SGWB and distinguish
separate components, the frequency dependence of the isotropic part of
$\Omega_{\rm GW}$ will be analyzed [12]. Secondly, a measurement of any
anisotropies in the SGWB will also be of interest for characterization and
component separation. In particular, any directional dependence will assist in
distinguishing between galactic and extra-galactic components and searching
for any correlation with known tracers of structure [14, 15, 16, 17].
Furthermore, anisotropies in the astrophysical background are expected to be
correlated with the large-scale structure due to how the GW originate and how
they propagate through a perturbed universe [18, 19, 20, 21]. Thirdly, crucial
information on the nature of the SGWB might be obtainable from studying its
polarization. Some cosmological mechanisms, for instance, a coupling between a
pseudo-scalar inflaton and a gauge field [22], can result in a fully
circularly polarized SGWB [23, 24]. A degree of polarization might also be
present in the astrophysical SGWB due to Poisson fluctuations in the (finite)
number of unresolved sources [25].

It is worth noting that an SGWB could be, in general, non-stationary in that
its statistical properties could vary in time due to the non-stationarity of
the source. In this case, the analysis of such a background is still possible
in principle [26]. For our purpose, we assume that the non-stationary signal
of the SGWB can be analyzed and removed from the data stream, considering the
fully stationary component of the background.

Previous work [27] studied the instrumental response to an anisotropic SGWB by
two L-shaped ground-based detectors. The formalism was extended in [28] and
[29] to include the case of a polarized SGWB. This paper extends the
aforementioned formalism to a more general one. We consider anisotropies and
an asymmetry in polarizations of the SGWB within a unified framework. We
generalize the results of [27] to a generic network of L-shaped and
triangular-shaped detectors. We give explicit analytic expressions for the
response of the network to the monopole, the dipole, and the quadrupole in a
spherical harmonics expansion of the SGWB, extending some results of [29]
(these expressions significantly speed up the numerical forecasts). We provide
a pedagogical derivation of the variance of the measurement, also paying
attention to the contributions for the variance of the signal, to assert the
level of accuracy that can ever be ideally achieved by a given network
configuration of arbitrary low noise [30]. Finally, we apply this formalism to
forecast the measurement of the anisotropies and the polarization of the SGWB
with a network of existing and third-generation ground-based interferometers.

The paper is structured so readers only interested in the results can find a
self-contained and complete presentation and discussion in section 2. The
remainder of the paper serves to introduce and develop the mathematical tools
necessary for this forecast. In section 3, we review the formalism for
describing SGWBs. Sections 4 and 5 introduce the formalism and notation for
the observables involved in ground-based networked detectors. The formalism is
given for a generic network made of an arbitrary number of triangular-shaped
detectors and L-shaped detectors. Section 6 applies this formalism to the
detection of anisotropic multipoles of the SGWB, whilst section 7 includes a
similar analysis for the detection of a net polarization. In both cases, we
individually highlight the signal- and noise-dominated regimes before
obtaining general Fisher forecasts for the observational uncertainty. Many of
the technical derivations underlying the main results are elaborated in
appendices A-C.

##  2 Main results

This section presents our main results on the forecast of measurements of the
anisotropies and the net polarization of the SGWB with a network of ground-
based interferometers. Our formalism assumes a generic set of L-shaped and
triangularly-shaped GW detectors, with noise specified by eqs. (4.3), (4.4),
and (4.5). In our evaluations and explicit studies, we consider a noise model
according to the advanced versions of the existing LIGO, Virgo, and Kagra
(denoted collectively as LVK in our figures) and the design sensitivity of
Einstein Telescope (ET) and Cosmic Explorer (CE). We show the noise Power
Spectral Densities for these instruments in figure 1, where the numerical
values are taken from ref. [31]. Our numerical results are obtained by
integrating in frequency in the range $f\in[1,1000]$ Hz.

Figure 1:  Predicted (square root of the) Power Spectral Density (PSD) of
Advanced LIGO, Advanced Virgo, KAGRA, Einstein Telescope, and Cosmic Explorer
as a function of the frequency. Data were taken from ref. [31]

The section is divided into three parts. Subsections 2.1 and 2.2 present our
results for, respectively, anisotropies and a net polarization of the SGWB.
These results are discussed in subsection 2.3, which also serves as a
concluding section of our paper, as the remainder of this work is devoted to
deriving the methodology used to obtain the results presented here.

###  2.1 Results: forecast observations of anisotropic modes

We consider a Gaussian SGWB characterized by

$\left\langle
h^{*}(f,\hat{n})h(f,\hat{n})\right\rangle=\frac{3H_{0}^{2}}{32\pi^{3}f^{3}}\,\Omega_{0}\left(\frac{f}{f_{0}}\right)^{\alpha}\sum_{\ell
m}\delta^{GW}_{\ell m}Y_{\ell m}(\hat{n})\,,$ (2.1)

which has a factorized dependence on frequency $f$ and on the GWs arrival
direction ${\hat{n}}$. The parameter $\Omega_{0}$ is the standard fractional
energy density in GWs at the pivot frequency $f_{0}$ (in our analysis, we fix
$f_{0}=100$ Hz). In our explicit evaluations, we fix the spectral tilt
$\alpha$ to two typical values considered in the literature, namely to
$\alpha=0$, as for a (nearly) scale-invariant cosmological inflationary
signal, and to $\alpha=2/3$, as expected for the stochastic background due to
the black hole-black hole and black hole-neutron star binary system inspirals
[26].

The last factor in eq. (2.1) encodes the non-trivial angular structure in
terms of a sum over spherical harmonics $Y_{\ell m}$ (where the $z-$axis is
taken to coincide with the rotation axis of the Earth), which generalizes the
more common isotropic studies in which only the monopole term is assumed. Due
to the normalization of the spherical harmonics, we fix
$\delta_{00}=\sqrt{4\pi}$ so that the coefficients $\delta^{GW}_{\ell m}$
encode the relative strength of the various multipoles relative to the
monopole. In our explicit analysis, we assume that only one given multipole is
present, with a magnitude much smaller than one, corresponding to a nearly
isotropic signal.

For each interferometer $i$ in the network (and for each channel in the case
of triangular-shaped detectors, which are made of multiple instruments), we
take the data stream $m_{i}(t)$ that results from the superposition of the GW
signal and instrumental noise,

$\displaystyle m_{i}(t)=s_{i}(t)+n_{i}(t)\;.$ (2.2)

We (window-)Fourier transform these functions and take every possible cross-
correlation of the data streams at different sites. We then integrate the
cross-correlations in time, using the assumption of a stationary SGWB, so that
the statistics underlying the response of the network of detectors is periodic
with period $T_{e}=1$ day. We apply optimal filters, building a set of
unbiased observables

${\cal
C}_{ij,m}\propto\int_{-\infty}^{\infty}dfm_{i}(f)m_{j}(f)\,Q_{m}(f)\;\;,\;\;i\neq
j\;,$ (2.3)

whose expectation value depends on the parameters in eq. (2.1). We consider
typical values $\alpha=0$ (from cosmology) and $\alpha=2/3$ (from
astrophysics) for the tilt of the signal. We construct a $\chi^{2}$ as a
function of the amplitude of the signal $\Omega_{0}$ and of the anisotropy
parameter $\delta_{\ell m}^{GW}$ under consideration (namely, for a specific
choice of $\ell$ and $m$). From the $\chi^{2}$, we forecast the uncertainties
under which the network will measure these two quantities under consideration.

In our explicit evaluations, we consider a fiducial isotropic signal,
$\delta^{GW}_{\ell m}=0$ (apart from the monopole term), with amplitude
${\bar{\Omega}}_{0}$. We then compute the $\chi^{2}$ as discussed above, as a
function of $\Omega_{0}$ and $\delta^{GW}_{\ell m}$. In figure 2 we show the
$1-,2-$, and $3-\sigma$ contours obtained from, an increase of, respectively,
$\delta\chi^{2}=1,\,4,\,9$, with respect to the value obtained for the
fiducial signal. 222The coefficients $\delta^{GW}_{\ell m}$ of the spherical
harmonics decomposition of eq. (2.1) are in general complex. The limits
obtained in our analyses refer to their absolute values. A measurement outside
these contours would reject the hypothesis of isotropy of the signal at the
corresponding confidence level. The different panels of the figure correspond
to a different multipole analysed and to a different choice of spectral tilt.
In all cases the fiducial amplitude is ${\bar{\Omega}}_{0}=10^{-9}$, which is
close to the current upper bound [6]. In all panels, we consider measurements
involving third-generation detectors (either correlated with the existing
second-generation ones or among themselves), as measurements from the current
detectors alone, do not constrain the anisotropy parameters below one.

Figure 2:  Forecast $1\sigma$, $2\sigma$ and $3\sigma$ contours as a function
of the amplitude $\Omega_{0}$ and of one anisotropy parameter
$\delta^{GW}_{\ell m}$ for a fiducial isotropic model with fiducial amplitude
${\bar{\Omega}}_{0}=10^{-9}$, for an observation time of $T=10$ yrs in the
LVK+ET (solid lines) and in the ET-CE network (dotted lines) network. The left
(respectively, right) panels study the dipole coefficient $\ell=m=1$ (resp.,
the quadrupole coefficient $\ell=m=2$). The upper (resp., lower) panels assume
a scale-invariant signal (resp., a blue signal with $\alpha=2/3$).

Marginalizing the results of such analyses over $\Omega_{0}$, we can then
obtain the uncertainty on the measurement of the anisotropy parameter alone.
In figure 3, we plot the forecast error in the measurements of a dipole
($\ell=m=1$) and a quadrupole ($\ell=m=2$) coefficient, as a function of the
fiducial amplitude of $\Omega_{0}$ and for different networks.

The vertical line in each panel shows the current upper limit on $\Omega_{0}$
[6]. We also show results in the excluded range of $\Omega_{0}$ because we
want to observe the change of slope in the forecast uncertainty of the
anisotropy. The left (resp., right) portion of the plot is in a regime in
which the instrumental noise (resp., the signal) dominates in the measurement
(and, at the technical level, in the variance that is present at the
denominator of the $\chi^{2}$).

One specific targeted anisotropy that we study in figure 3 is the one arising
from our particular motion. If we assume that the SGWB is isotropic in the CMB
rest-frame, then our peculiar motion with velocity $v$ in this frame produces
an anisotropic signal analogous to that in the CMB, characterized by a
kinematic dipole (suppressed by $v$ wrt the monopole), kinematic quadrupole
(suppressed by $v^{2}$) and so on [32] 333In principle, the kinematic effect
also changes the frequency dependence of the signal in the various multipoles,
invalidating the factorization of eq. (2.1)). There are two possible
contributions to this, one arising from the spectral dependence of
$\Omega_{GW}\left(f\right)$ and one from an aberration that mixes different
multipoles. The first effect is not present for a power-law
$\Omega_{GW}\left(f\right)$ [32], as we consider in eq. (3.15). The second
effect is suppressed by the smallness of the anisotropies (as it is
proportional to spatial gradients of the SGWB energy density [33]). Therefore
our analysis of the kinematic effect based on eq. (2.1) remains accurate. For
more general frequency shapes and for a large anisotropic astrophysical SGWB,
see [33, 13]. As $v/c\simeq 1.2\times 10^{-3}$ [34], the dipole is the largest
coefficient induced by this kinematic effect. In the frame of reference chosen
by our decomposition in multipoles, the kinematic dipole aligns closely with
the $\ell=1$, $m=-1$ multipole, and it has a predicted value of
$\delta^{GW}_{1,-1}\simeq 1.7\times 10^{-3}$ [27], that we indicate with a
horizontal line in figure 3.

Higher multiples are studied in figure 4, where we show the forecast errors on
various $\delta_{\ell m}^{GW}$, after marginalizing over $\Omega_{0}$. Each
line in the figure is obtained from a separate analysis. For each $\ell$
varying from $1$ to $8$, we make separate analyses for each value of $m$, and
we then show in the figure the result for the value of $m$ that is best
reconstructed for that $\ell$.

Figure 3: Forecast error in the measurement of the anisotropic coefficient
$\delta_{\ell m}^{GW}$ as a function of the fiducial amplitude
$\bar{\Omega}_{0}$ of the signal for a total observation time of $T=10$ yr and
for different networks. The left (resp. right) panels study the dipole
coefficient $\ell=m=1$ (resp., the quadrupole coefficient $\ell=m=2$). The
upper (resp., lower) panels assume a scale-invariant signal (resp., a blue
signal with $\alpha=2/3$). The horizontal and vertical lines indicate,
respectively, the expected value of the kinematic dipole and the current upper
bound on the amplitude of $\Omega_{0}$ [6].
Figure 4:  Forecast error in the measurement of the coefficients $\delta_{\ell
m}$ as a function of the fiducial amplitude $\bar{\Omega}_{0}$. Each line
corresponds to a different $\ell$ and the value of $m$ that results in the
smallest $\sigma_{\delta_{\ell m}}$ for that $\ell$). An observation time
$T=10$ yr with the LVK + ET (left panel) and the ET + CE (right panel) network
is assumed. The vertical line indicates the current upper bound on the
amplitude of $\Omega_{0}$ [6].

###  2.2 Results: forecast observations of net polarization

We now consider a Gaussian SGWB characterized by

$\left\langle
h_{\lambda}^{*}(f,\hat{n})h_{\lambda}(f,\hat{n})\right\rangle=\frac{3H_{0}^{2}}{32\pi^{3}f^{3}}\,\Omega_{0}\left(\frac{f}{f_{0}}\right)^{\alpha}\left(1+p\lambda\right)\,,$
(2.4)

where the parameter $p$ encodes the amount of net circular polarization of the
SGWB, this parameter ranges from $-1$, corresponding to a left-handed
background, to $+1$, corresponding to a right-handed background. Within this
range, the choice $p=0$ corresponds to an unpolarized background.

The study of the net polarization proceeds in the same way as for the
anisotropic search. We form cross-correlations of the data streams at
different sites and build an unbiased observable as a function of the two
varying parameters $\Omega_{0}$ and $p$. figure 5 shows the $1-$, $2-$, and
$3-\sigma$ contours in the $\left\\{\Omega_{0},\,p\right\\}$ plane for
measurements with the LVK+ET and ET+CE networks. As for the anisotropy, we do
not show contour plots from the LVK network alone as they do not yield
significant constraints on $p$.

Figure 5: Forecast $1\sigma$, $2\sigma$ and $3\sigma$ contours as a function
of the amplitude $\Omega_{0}$ and of the polarization parameter $p$ for a
fiducial amplitude ${\bar{\Omega}}_{0}=10^{-9}$, for an observation time of
$T=10$ yrs in the LVK+ET (solid lines) and in the ET-CE network (dotted lines,
barely visible in the two panels). The left (resp., right) panel assumes a
fiducial unpolarized signal with $\alpha=2/3$, (resp., a fully circularly
polarized scale-invariant signal).

In figure 6, we then show forecasts on the net polarization factor $p$ for
different values of the fiducial amplitude $\Omega_{0}$, after marginalizing
over $\Omega_{0}$ in the analysis. The left panel (resp., right panel) assumes
an astrophysical unpolarized background (resp., a cosmological fully polarized
background), and the different lines shown correspond to different network
configurations.

Figure 6: Forecast error in the measurement of the polarization factor $p$ as
a function of the fiducial amplitude $\bar{\Omega}_{0}$ for a total
observation time of $T=10$ yrs. and different networks. The left (resp.,
right) panel assumes an unpolarized fiducial signal with $\alpha=2/3$, (resp.,
a fully circularly polarized scale-invariant signal). The vertical line
indicates the current upper bound on the amplitude of $\Omega_{0}$ [6].

###  2.3 Discussion

In this work, we outline in full detail a mathematical framework to forecast
the sensitivity of a network of L-shaped and triangular-shaped ground-based
detectors to an SGWB with anisotropies or with net circular polarization.
Above, we have shown some results obtained when this formalism is applied to
different networks of present and next-generation ground-based detectors.

Our study of anisotropies in subsection 2.1 assumes that the monopole
dominates the SGWB signal, which is expected to be true for the vast majority
of astrophysical and cosmological models [32]. Moreover, we assumed that,
besides the monopole, only one additional multipole is present in the
expansion of eq. (2.1), where, parametrically, the coefficient
$\delta^{GW}_{\ell m}$ controls the ratio between the (directionally-
dependent) contribution to the fractional SGWB energy density (per logarithmic
interval in frequency) of that multipole relative to the monopole. The
amplitude of the monopole term is controlled by the standard parameter
$\Omega_{0}$, which gives the fractional GW energy density (per logarithmic
interval in frequency) $\Omega_{0}$ at the pivot scale $f_{0}=100\,{\rm Hz}$.
Our analysis provides the forecast uncertainties on determining the two
parameters $\Omega_{0}$ and $\delta^{GW}_{\ell m}$. We also presented the
uncertainty on $\delta^{GW}_{\ell m}$ by marginalizing the results of this
analysis over $\Omega_{0}$.

Not surprisingly, as visible in figures 3 and 4, the sensitivity to the
relative size of the anisotropic component $\delta^{GW}_{\ell m}$ worsens as
the overall magnitude of the signal decreases, namely as we decrease the
fiducial value of $\Omega_{0}$. More precisely, as the variance of the
measurement is noise-dominated in this range, the sensitivity scales as the
inverse of this quantity. As a reference, we show in the figures the current
limit on $\Omega_{0}$ [6]. The SGWB obtained in astrophysical models may be
close to this upper limit [35, 36], while predictions for the cosmological
SGWB are much more uncertain and model dependent [37]. Our results indicate
that if the SGWB saturates the current bound, the existing LIGO-Virgo-KAGRA
network at the designed sensitivity will barely be able to constrain even
anisotropies of order one significantly. On the other hand, the next
generation detectors, such as Einstein Telescope and Cosmic Explorer, might be
able to probe a ${\cal O}\left(10^{-3}\right)$ anisotropy at the largest
angular scales (lowest $\ell$; figure 4 shows how the sensitivity worsens for
higher multipoles).

As a comparison, anisotropies in the astrophysical SGWB can be estimated to be
of the order of $10^{-2}$ [38, 39] The anisotropies due to the propagation of
primordial GWs are expected to be below this level [18, 19, 20, 12], with a
quadrupole amplitude $\delta^{GW}_{2m}\simeq\sqrt{C_{2}}\simeq 3\times
10^{-5}$ [32]. 444Cosmological perturbations might be greater in certain
cases, see e.g. [40, 41, 42]. Finally, the kinematic dipole, induced by the
peculiar motion of the Earth, is of ${\cal O}\left(10^{-3}\right)$, as we
already discussed above in relation to figure 3.

To summarize, the anisotropies induced by astrophysical sources and the
kinematic effect might be probed by the Einstein Telescope - Cosmic Explorer
combination, while the measurement of the anisotropy from the cosmological
propagation appears to be more challenging. Significantly improving over this
will require more than just decreasing the instrumental noise. Even in the
ideal limit of vanishing instrumental noise, the sensitivity to the
subdominant anisotropic component is limited by the variance of the dominant
monopole contribution [30]. The curves shown in figures 3 and 6 show a marked
change in their slope from the left part at small ${\bar{\Omega}_{0}}$ and the
right part at large ${\bar{\Omega}_{0}}$. This change in the slope is due to
the transition between a regime in which the variance of the measurement is
dominated by the instrumental noise to a regime in which it is dominated by
the signal.555While both contributions to the variance are fully included in
the results presented in this section, in subsection 5.2 we show that the
variance can be written in terms of relatively simple analytical expressions
in the limits in which it is fully dominated by either contribution. For the
network configurations we have discussed, the measurement is in the noise-
dominated regime once the signal is constrained to be below the current upper
limit. Assuming that futuristic experiments will reach an ideal level of
negligible instrumental noise, the variance of their measurements will be due
solely to that of the signal. The sensitivity to the anisotropy will then
reach the level that can be seen from the high-${\bar{\Omega}}_{0}$ portion of
our figures. Our results assume an observation period of $10$ years. Improving
over these limits will then require increasing the observation period $T$, or
the number of detectors $N$ [30]. We recall that the expected error decreases
as the square root of the measurements. Therefore it scales as $T^{-1/2}$ and
as $N^{-1}$ since the number of cross-correlations between different detectors
scales as $N^{2}$. 666This statement is strictly true in the limit of
vanishing contribution to the variance of measurements taken at different
pairs. As we discuss in subsection 5.2, we verified that the off-diagonal
contributions of the variance between different pairs are indeed negligible
with respect to the diagonal ones between the same pair (see the structure in
eq. (5.14). We expect that this would no longer be true in the limit of large
$N$.

In subsection 2.2, we presented the results of a study of the forecast
sensitivity to a net circular polarization, which, as we discussed in the
Introduction, is predicted to be of order one in some inflationary mechanisms,
and which might also be present, although to a significantly smaller amount,
in the astrophysical signal. figure 6 shows that if the amplitude of the SGWB
is close to the current limit, an order one net polarization can be marginally
probed by the existing detectors at design sensitivity. The Einstein Telescope
- Cosmic Explorer network can improve over this by about three orders of
magnitude. As shown in the figure, the forecast sensitivities worsen as the
amplitude $\Omega_{0}$ decreases below the current limit.

Our results are in agreement with [43], where a frequency-dependent power
spectrum and a generic power law tilt have been investigated, and with [25],
where the predicted values of $\Omega_{0}$ for an astrophysical background and
the noise curves (as in figure 1 of that paper) produce an expected SNR for a
network made of ET and CE which is in agreement with our forecast error in the
measurement of a net amount of polarization. A more precise comparison with
these results, particularly for what concerns [43], may be done by allowing
the polarization parameter $p$ and the power spectrum tilt $\alpha$ to depend
on frequency, which is feasible in principle in our case. However, the
comparison between our two panels in figure 6 indicates that the result does
not vary very much in the realistic region of the parameter space, thus
supporting the agreement between our results and those of [43] on the
polarization.

The analyses carried out in this paper can be extended in several directions.
Our formalism can be immediately employed for any planned network and to
compare the performance of different configurations (as done, for instance, in
the recent ET study [44]). Another natural step will be to go beyond the
assumption of factorized angular frequency dependence of the SGWB, along with
considering simultaneously more than a single multipole $(\ell,m)$ in the
analysis of the anisotropies. It may also be interesting to produce the same
kind of forecasts for LISA [45], the space-based triangular-shaped GW
interferometer planned to be launched on $\sim$2037, improving over the
estimates of [12] that assume statistical isotropy.

##  3 The stochastic gravitational wave background

We summarise here the basic formalism for describing an anisotropic SGWB
following the notation and conventions of earlier work [27, 28]. We consider
gravitational waves at $x^{\alpha}$, the location of the detectors, as small
perturbations of the metric tensor on a Minkowski background in such a way
that

$g_{\mu\nu}(x^{\alpha})=\eta_{\mu\nu}+h_{\mu\nu}(x^{\alpha})\,,$ (3.1)

where $\eta_{\mu\nu}$ is the (mostly positive signature) metric tensor for
Minkowski spacetime and $h_{\mu\nu}$ is the gravitational wave, which is a
small perturbation of the background metric, $h\ll\eta$.

It is convenient to completely fix the gauge freedom through the standard
choice
$\partial_{\mu}h^{\mu}_{\nu}-\frac{1}{2}\partial_{\nu}\,h^{\mu}_{\mu}=h_{00}=h_{0a}=0$
(greek indices run from 0 to 3, while latin indices run from 1 to 3). With
this choice, the GW term has non-vanishing entries only in the spatial
components, where it is transverse ($\partial_{a}h_{ab}=0$) and traceless
($h_{aa}=0$), and where we decompose in plane waves

$\displaystyle h_{ab}(t,\vec{x})$ $\displaystyle=\int_{-\infty}^{\infty}df\int
d^{2}{\hat{n}}\,e^{2\pi
if\left(t-{\hat{n}}\cdot\vec{x}\right)}\sum_{s=+,\times}h_{s}(f,\hat{n})\,e^{s}_{ab}\left({\hat{n}}\right)\equiv\sum_{s}h_{ab,s}(t,\vec{x})\,,$
(3.2)

where $\hat{n}$ is the unit vector in $\mathrm{R}^{3}$, and
$e^{+}_{ab}(\hat{n})$, $e^{\times}_{ab}(\hat{n})$ are the two polarization
basis tensors, defined by

$\displaystyle
e_{ab}^{+}(\hat{n})=p_{a}p_{b}-q_{a}q_{b}\;\;,\;\;e_{ab}^{\times}(\hat{n})=p_{a}q_{b}+q_{a}p_{b}\,.$
(3.3)

For our choice of spherical coordinates, the direction vector is defined as

$\displaystyle\hat{n}=\cos\phi\sin\theta\,\hat{x}+\sin\phi\sin\theta\,\hat{y}+\cos\theta\,\hat{z}\,.$
(3.4)

The direction vector is aligned with the direction of propagation of the plane
wave ($\theta,\phi$ are polar angular coordinates, while
$\hat{x},\hat{y},\hat{z}$ are unit vectors along the three Cartesian axes),
while $\hat{p}$ and $\hat{q}$ are two directions perpendicular to $\hat{n}$
and to each other, that are conventionally chosen as in [27]

$\displaystyle\hat{p}=\sin\phi\hat{x}-\cos\phi\hat{y}\;\;,\;\;\hat{q}=\cos\phi\cos\theta\hat{x}+\sin\phi\cos\theta\hat{y}-\sin\theta\hat{z}\,.$
(3.5)

One can also consider right-handed ($\lambda=+1$) and left-handed
($\lambda=-1$) circular polarizations rewriting eq. (3.2) as

$\displaystyle h_{ab}(t,\vec{x})=\int_{-\infty}^{\infty}df\int
d^{2}{\hat{n}}\,e^{2\pi
if\left(t-{\hat{n}}\cdot\vec{x}\right)}\sum_{\lambda=\pm
1}h_{\lambda}(f,\hat{n})\,e^{\lambda}_{ab}\left({\hat{n}}\right)\;\;,\;\;e^{\lambda}_{ab}\equiv\frac{e^{+}_{ab}+i\,\lambda\,e^{\times}_{ab}}{\sqrt{2}}\,.$
(3.6)

The plane wave solution of eq. (3.6) is defined in a reference system fixed
with respect to the expanding cosmological background and in which the spatial
pattern of the perturbations of the stochastic background is assumed to be
statistically time-independent. The $\hat{z}$ vector of this system is chosen
to point along the direction of the Earth’s rotation axis. To describe a
stochastic source, we treat the complex amplitude $h_{\lambda}(f,\hat{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_{\lambda}(f,\hat{n})$ as

$\left\langle
h^{*}_{\lambda}\left(f,{\hat{n}}\right)h_{\lambda^{\prime}}\left(f^{\prime}{\hat{n}}^{\prime}\right)\right\rangle=\delta_{\lambda\lambda^{\prime}}\delta_{D}^{(2)}\left({\hat{n}}-{\hat{n}}^{\prime}\right)\delta_{D}\left(f-f^{\prime}\right)\mathcal{H}_{\lambda}\left(\left|f\right|,{\hat{n}}\right)\,,$
(3.7)

where $\delta_{D}^{(2)}\left({\hat{n}}-{\hat{n}}^{\prime}\right)$ is a
covariant two-dimensional Dirac delta-function on the unit two-sphere and
$\delta_{D}\left(f-f^{\prime}\right)$ the Dirac delta function on the
frequency space.

In principle, such a source has spectral properties which arbitrarily depend
upon amplitude and frequency. For simplicity, here we consider a factorized
dependence

$\mathcal{H}_{\lambda}\left(\left|f\right|,{\hat{n}}\right)=H_{\lambda}\left(\left|f\right|\right)P_{\lambda}\left({\hat{n}}\right)\,,$
(3.8)

which amounts to assuming that sources of the SGWB situated along a different
line of sights emit with the same ”average” spectrum in frequency. It is
conventional to define the fractional energy density in GW per logarithmic
frequency interval and per polarisation as 777In the decompositions of eqs.
(3.2) and (3.6) the frequency can have both positive and negative values. The
power $\mathcal{H}_{\lambda}$ in eq. (3.7) is a function only of the magnitude
of the frequency so that one can conveniently define the energy density
$\rho_{GW,\lambda}$ and the fractional energy density $\Omega_{GW,\lambda}$ in
terms of positive frequencies only.

$\Omega_{GW,\lambda}\left(f\right)\equiv\frac{f}{\rho_{\rm
crit}\left(f\right)}\frac{d\rho_{GW,\lambda}}{df}=\sum_{\lambda}\frac{16\pi^{3}}{3H_{0}^{2}}f^{3}H_{\lambda}\left(|f|\right)\,,$
(3.9)

where $\rho_{\rm 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

$P_{\lambda}(\hat{n})=\sum_{\ell=0}^{+\infty}\sum_{m=-\ell}^{\ell}\,\delta^{GW}_{\ell
m,\lambda}\,Y_{lm}(\hat{n})\;\;\;,\;\;\;\delta^{GW}_{00,\lambda}=\sqrt{4\,\pi}\,,$
(3.10)

where the value of $\delta^{GW}_{00,\lambda}$ is set considering the
convention for the monopole $Y_{00}=\frac{1}{\sqrt{4\pi}}$.

In this expression, we have allowed the coefficients of the decomposition
$\delta^{GW}_{\ell m,\lambda}$ to depend on the polarization. However, due to
the lack of strong theoretical motivation for the generation of large-scale
polarization-dependent anisotropies, for simplicity, we assume that

$\delta^{GW}_{\ell m,1}=\delta^{GW}_{\ell m,-1}\equiv\delta^{GW}_{\ell
m}\;\;\Rightarrow\;\;P_{+1}\left({\hat{n}}\right)=P_{-1}\left({\hat{n}}\right)\equiv
P\left({\hat{n}}\right)\;.$ (3.11)

On the other hand, we allow the monopole term to depend on the polarization in
a frequency-independent amount:

$\displaystyle
H_{\lambda}\left(|f|\right)\equiv\left(1+\lambda\,p\right)\,H\left(|f|\right)\;\;\;,\;\;\;-1\leq
p\leq 1\;,$ (3.12)

where $p=+1$ (respectively, $p=-1$) corresponds to a fully right-handed
(respectively, left-handed) circularly polarized SGWB, while $p=0$ corresponds
to an unpolarized background. With this parametrization, the quantity
$H\left(|f|\right)$ can be immediately related to the GW fractional energy
density with the standard expression

$\Omega_{GW}(f)\equiv\sum_{\lambda}\Omega_{GW,\lambda}(f)=\frac{32\pi^{3}}{3H_{0}^{2}}f^{3}H(|f|)\;,$
(3.13)

irrespective of the degree of polarization of the SGWB (namely, for any value
of $p$). We also note that one can write

$\displaystyle\Omega_{GW,\lambda}\left(f\right)=\frac{1+\lambda\,p}{2}\,\Omega_{GW}\left(f\right)\;.$
(3.14)

Finally, we assume that, in the range which affects the measurements, the GW
energy density has a power-law frequency dependence

$\Omega_{\rm
GW}\left(f\right)={\Omega}_{0}\left(\frac{f}{f_{0}}\right)^{\alpha}\;,$ (3.15)

where ${\Omega}_{0}$ is the fractional energy density at the pivot scale of
$f_{0}=100\,{\rm Hz}$, and where, as discussed in subsection 2.1 we fix
$\alpha=0,\,2/3$. From the above relations, it is clear that the same power-
law dependence applies to $H_{\lambda}\left(|f|\right)$.

##  4 Data stream of ground-based interferometers

As mentioned above, we consider a network of $N$ L-shaped ground-based
detectors (for instance, LIGO, Virgo, and KAGRA) and $M$ ground-based
detectors with triangular configuration (for instance Einstein Telescope). For
each L-shaped interferometer of the network, we measure the difference $\Delta
T$ for light-travel times across separate arm paths. The measurement is
affected by the instrument noise and any GW signal present,

$m_{i}\left(t\right)=\frac{\Delta
T_{i}}{T_{0}}=n_{i}\left(t\right)+s_{i}\left(t\right)\;\;\;,\;\;\;i=1,\dots,N\;,$
(4.1)

where $T_{0}$ is the time the laser would need to travel back and forth
without any signal and noise.

We define the set of measurements made by each triangular-shaped
interferometer of the network similarly. In this case, each triangular
detector has an equilateral triangular configuration with three Michelson
interferometers at its vertices $\alpha=X,\,Y,\,Z$, therefore

$m_{i\alpha}\left(t\right)=\frac{\Delta
T_{i\alpha}}{T_{0}}=n_{i\alpha}\left(t\right)+s_{i\alpha}\left(t\right)\;\;\;,\;\;\;i=1,\dots,M\;\;\;,\;\;\;\alpha=X,\,Y,\,Z\;.$
(4.2)

We then assume that the noise of each network instrument is not correlated
with the one affecting the others. We also assume that the noise of each
detector is Gaussian distributed with zero mean. These assumptions imply a
vanishing correlation between L-shaped and triangular-shaped detectors

$\left\langle{\tilde{n}}_{i\alpha}^{*}\left(f\right){\tilde{n}}_{j}\left(f^{\prime}\right)\right\rangle\equiv
0\;,$ (4.3)

where $i=1,\dots,M$ and $j=1,\dots,N$. For the L-shape detectors, we can then
write

$\left\langle{\tilde{n}}_{i}^{*}\left(f\right){\tilde{n}}_{j}\left(f^{\prime}\right)\right\rangle\equiv\frac{\delta_{ij}}{2}\delta_{D}\left(f-f^{\prime}\right)P_{i}\left(\left|f\right|\right)\;,$
(4.4)

where $i=1,\dots,N$ and $j=1,\dots,N$.

For the variance of the noise of the triangular-shaped detectors, accounting
for the fact that every interferometer shares one arm with each of the other
two interferometers (of the same detector), we write

$\left\langle{\tilde{n}}_{i\alpha}^{*}\left(f\right){\tilde{n}}_{j\beta}\left(f^{\prime}\right)\right\rangle\equiv\frac{\delta_{ij}}{2}\delta_{D}\left(f-f^{\prime}\right)N_{\alpha\beta,i}\left(\left|f\right|\right)\;,$
(4.5)

where $i=1,\dots,M$, $j=1,\dots,M$. If the $i$-th detector has an exact
equilateral configuration, with three identical instruments at the vertices,
the noise matrix at each site is a $3\times 3$ matrix having the same element
$N_{d,i}$ in all diagonal entries and the same element $N_{o,i}$ in all off-
diagonal entries. This matrix can be diagonalized in the three channels

$\displaystyle m_{iO}\equiv
c_{O\alpha}\,m_{i\alpha}\;\;,\;\;\;\;O=A,\,E,\,T\;\;,\;\;\alpha=X,\,Y,\,Z\,,$
(4.6)

(we note that the same combinations $c_{O\alpha}$ are taken in each detector
$i$) with

$c\equiv\left(\begin{array}[]{ccc}\frac{4}{3\sqrt{3}}&-\frac{2}{3\sqrt{3}}&-\frac{2}{3\sqrt{3}}\\\
0&-\frac{2}{3}&\frac{2}{3}\\\
\left(\frac{2}{3}\right)^{3/2}&\left(\frac{2}{3}\right)^{3/2}&\left(\frac{2}{3}\right)^{3/2}\end{array}\right)\;.$
(4.7)

These linear combinations, already introduced in [46] for the LISA experiment,
which also has an equilateral configuration (the explicit frequency dependence
of $N_{d,i}$ and $N_{o,i}$ for LISA can be found for example in ref. [47])
were used for ET in [28].

We use these channels as they diagonalize the noise matrix, which simplifies
the computation of the SNR that we perform below. In fact, combining eqs.
(4.5) and (4.6), one obtains

$\displaystyle\left\langle{\tilde{n}}_{iO}^{*}\left(f\right){\tilde{n}}_{jO^{\prime}}\left(f^{\prime}\right)\right\rangle=\frac{1}{2}\delta_{D}\left(f-f^{\prime}\right)\,\delta_{ij}\,\delta_{OO^{\prime}}\,N_{O,i}\left(\left|f\right|\right)\;,$
(4.8)

with

$N_{A,i}\left(f\right)=N_{E,i}\left(f\right)=\frac{8}{9}\left[N_{d,i}\left(f\right)-N_{o,i}\left(f\right)\right]\;\;,\;\;N_{T,i}\left(f\right)=\frac{8}{9}\left[N_{d,i}\left(f\right)+2N_{o,i}\left(f\right)\right]\;.$
(4.9)

As in [28], we consider this choice of channels in the present work since they
diagonalize the noise matrix. Furthermore, as explicitly verified in [28], the
$A-$ and $E-$channels behave as $90^{\circ}$ interferometers for small
frequencies (a limit that holds in the relevant frequency range for Einstein
Telescope and instruments of similar dimensions 888The condition for the low
frequency / short arm limit is $2\pi fL\ll 1$, where $L$ is the arm length. In
the diagonalization of eq. (4.7) we changed the normalization of each channel
with respect to [46], so that, in the low frequency / short arm limit, the
$A-$ and $E-$channels behave as $90^{\circ}$ degrees interferometers, with
unit-norm arm factors $d_{iA,E}^{ab}$ given in eq. 5.11.. In this regime, the
contribution of the $T-$channel to the signal vanishes, so we disregard it in
our computations [32].

##  5 Measurement by a network of detectors

Similarly to what has been done in [27] and [28], we start by introducing the
time-dependent Fourier transform of eq. (4.1) for the measurement in the
$i$-th L-shaped interferometer

${\tilde{m}}_{i}\left(f,\,t\right)\equiv\int_{t-\tau/2}^{t+\tau/2}dt^{\prime}{\rm
e}^{-2\pi ift^{\prime}}\,m_{i}\left(t^{\prime}\right)\;.$ (5.1)

We label the measurements in triangular-shaped interferometers with an
additional index $O$ accounting for the channel considered (namely, the A- or
the E-channel for that interferometer), and we also transform the measurements
according to

${\tilde{m}}_{iO}\left(f,\,t\right)\equiv\int_{t-\tau/2}^{t+\tau/2}dt^{\prime}{\rm
e}^{-2\pi ift^{\prime}}\,m_{iO}\left(t^{\prime}\right)\;.$ (5.2)

In these relations, the integration in time is performed over a timescale
$\tau$ much greater than the inverse of the smallest frequency of the GWs that
we want to study but sufficiently small that we can neglect the effect due to
the rotation of the Earth during this time. We introduce the Fourier
transforms of the signals ${\tilde{s}}_{i}$ and ${\tilde{s}}_{jO}$, and of the
noise ${\tilde{n}}_{i}$, and ${\tilde{n}}_{jO}$ contributions in an analogous
manner.

We can then define the set of estimators

$\displaystyle{\cal C}_{ij}\left(t\right)$
$\displaystyle\equiv\int_{-\infty}^{\infty}df{\tilde{m}}_{i}^{*}\left(f,\,t\right){\tilde{m}}_{j}\left(f,\,t\right)\,{\tilde{Q}}_{ij}\left(f\right)\;,$
$\displaystyle{\cal C}_{ij,O}\left(t\right)$
$\displaystyle\equiv\int_{-\infty}^{\infty}df{\tilde{m}}_{iO}^{*}\left(f,\,t\right){\tilde{m}}_{j}\left(f,\,t\right)\,{\tilde{Q}}_{ij,O}\left(f\right)\;,$
$\displaystyle{\cal C}_{ij,OO^{\prime}}\left(t\right)$
$\displaystyle\equiv\int_{-\infty}^{\infty}df{\tilde{m}}_{iO}^{*}\left(f,\,t\right){\tilde{m}}_{jO^{\prime}}\left(f,\,t\right)\,{\tilde{Q}}_{ij,OO^{\prime}}\left(f\right)\;,$
(5.3)

where the functions ${\tilde{Q}}_{ij,OO^{\prime}}\left(f\right)$,
${\tilde{Q}}_{ij,O}\left(f\right)$ and ${\tilde{Q}}_{ij}\left(f\right)$ are
weights (in the sum over channels, interferometers of the network, and
frequencies) that will be chosen later to maximize the effect of the signal
over the noise. Note that in the definition of the estimator (5) we do not
consider the correlations between instruments located at the same site. This
is done because it is difficult to characterize all the noise contributions
for collocated instruments properly.

From now on, to shorten the notation, we use the expression $(,OO^{\prime})$
to indicate all possible combinations of L-shaped and triangular-shaped
detectors collectively. For instance, the three expressions (5) can be
collectively written as

${\cal
C}_{ij(,OO^{\prime})}\left(t\right)\equiv\int_{-\infty}^{\infty}df{\tilde{m}}_{i(,O)}^{*}\left(f,\,t\right){\tilde{m}}_{j(,O^{\prime})}\left(f,\,t\right)\,{\tilde{Q}}_{ij,(,OO^{\prime})}\left(f\right)\;\;\;,\;\;\;i\neq
j\;.$ (5.4)

We assume that the statistical properties of the signal and the noise do not
change with time. Then for an anisotropic SGWB, the statistics of the signal
are periodic, with periodicity given by the rotation period
$T_{e}=\frac{2\pi}{\omega_{e}}$ of the Earth. It is then convenient to perform
the Fourier decomposition

${\cal C}_{ij(,OO^{\prime})}\left(t\right)=\sum_{m=-\infty}^{\infty}{\cal
C}_{ij(,OO^{\prime}),m}\,{\rm e}^{im\omega_{e}t}\;\;\;\;,\;\;\;\;{\cal
C}_{ij(,OO^{\prime}),m}\equiv\frac{1}{T}\int_{0}^{T}dt\,{\rm
e}^{-im\omega_{e}t}\,{\cal C}_{ij(,OO^{\prime})}\left(t\right)\;,$ (5.5)

where we take the observation time $T$ to be an integer multiple of one day
$T_{e}$.

Under the assumption of Gaussianity, we can define the Likelihood function of
the estimated parameters of the theory $\theta$ and their fiducial value
$\bar{\theta}$ for the set of measurements ${\cal C}_{m}$ provided by the
network:

$\displaystyle\mathcal{L}(\theta,\bar{\theta})\equiv\frac{1}{N(\bar{\theta})}e^{-\frac{1}{2}\chi^{2}(\theta,\bar{\theta})}\;,$
$\displaystyle\chi^{2}(\theta,\bar{\theta})\equiv\sum_{m^{\prime},m^{\prime\prime}}\sum_{\begin{subarray}{c}i\neq
j\\\ (,OO^{\prime})\end{subarray}}\sum_{\begin{subarray}{c}k\neq l\\\
(,O^{\prime\prime}O^{\prime\prime\prime})\end{subarray}}\left\langle
r^{*}_{ij(,OO^{\prime}),m^{\prime}}(\theta,\bar{\theta})\right\rangle\Sigma_{ij(,OO^{\prime}),kl(,O^{\prime\prime}O^{\prime\prime\prime}),m^{\prime}m^{\prime\prime}}^{-2}(\bar{\theta})\left\langle
r_{kl(,O^{\prime\prime}O^{\prime\prime\prime}),m^{\prime\prime}}(\theta,\bar{\theta})\right\rangle\;,$
(5.6)

where we introduced the quantities

$\displaystyle\Sigma_{ij(,OO^{\prime}),kl(,O^{\prime\prime}O^{\prime\prime\prime}),m^{\prime}m^{\prime\prime}}^{2}(\hat{\theta})\equiv\left\langle
r^{*}_{ij(,OO^{\prime}),m^{\prime}}(\bar{\theta},\bar{\theta})\,r_{kl(,O^{\prime\prime}O^{\prime\prime\prime}),m^{\prime\prime}}(\bar{\theta},\bar{\theta})\right\rangle\;,$
$\displaystyle
r_{ij(,OO^{\prime}),m}(\theta,\bar{\theta})\equiv\mathcal{C}_{ij(,OO^{\prime}),m}(\theta)-\langle\mathcal{C}_{ij(,OO^{\prime}),m}(\bar{\theta})\rangle\;.$
(5.7)

We stress that the sums over $ij(,OO^{\prime})$,
$kl(,O^{\prime\prime}O^{\prime\prime\prime})$ are meant over all the possible
combinations of the observables built-in eq. (5). The $\chi^{2}$ in eq. (5) is
a well-defined $\chi^{2}$ since it has a global minimum in
$\theta=\bar{\theta}$.

The set of parameters for our analysis is
$\theta\equiv(\Omega_{0},\delta^{GW}_{\ell m},p,\alpha)$, where we recall that
$\Omega_{0}$ is the amplitude of the monopolar power spectrum as in eq.
(3.15), $\delta^{GW}_{\ell m}$ is the ratio between the ($\ell,m$) multipole
and the monopole, as in eq. (3.10), $p$ is the polarization factor as defined
in eq. (3.12), and $\alpha$ is the tilt of the power spectrum in eq. (3.15).

###  5.1 Expectation value of the signal

In appendix A we show that the expectation value of eq. (5.5) is

$\displaystyle\left\langle
C_{ij(,OO^{\prime}),m}\right\rangle=\sum_{\lambda}\frac{8\pi\,\tau}{5}\int_{-\infty}^{+\infty}df\,H_{\lambda}\left(|f|\right)\sum_{\ell=|m|}^{\infty}\delta^{GW}_{\ell
m,\lambda}\gamma_{\ell
m,ij(,OO^{\prime}),\lambda}\left(f\right)Q_{ij(,OO^{\prime}),m}\left(f\right)\;,$
(5.8)

where, following [27], the overlap functions are defined as 999The factor
$\frac{5}{8\pi}$ is conventional, and it has the purpose of eliminating an
overall factor in the monopole term; see [28] for further details.

$\displaystyle\gamma_{\ell
m,ij(,OO^{\prime}),\lambda}(f)\equiv\frac{5}{8\pi}\int d^{2}{\hat{n}}{\rm
e}^{2\pi if{\hat{n}}\cdot\Delta\vec{x}_{ij}}\,Y_{\ell
m}\left({\hat{n}}\right)e_{ab}^{*\lambda}\left({\hat{n}}\right)e_{cd}^{\lambda}\left({\hat{n}}\right)\,d_{i(,O)}^{ab}\,d_{j(,O^{\prime})}^{cd}\;.$
(5.9)

It is useful to rewrite the $\gamma$ coefficients of eq. (5.9) as

$\displaystyle\gamma_{\ell m,ij(,OO^{\prime}),\lambda}(f)\equiv\gamma_{\ell
m,\lambda,abcd}(\kappa_{ij},\hat{s}_{ij})\times
d_{i(O)}^{ab}\,d_{j(O^{\prime})}^{cd}\;,$ (5.10)

where we define the tensors

$\displaystyle
d_{i(\alpha)}^{ab}\equiv\frac{\hat{u}_{i(\alpha)}^{a}\hat{u}_{i(\alpha)}^{b}-\hat{v}_{i(\alpha)}^{a}\hat{v}_{i(\alpha)}^{b}}{2}\;\;,\;\;d^{ab}_{iO}\equiv
c_{O\alpha}\,d_{i\alpha}^{ab}\;,$ (5.11)

where $\hat{u}_{i}$, $\hat{v}_{i}$ are the arm directions of the L-shaped
interferometers, while $\hat{u}_{i\alpha}$, $\hat{v}_{i\alpha}$ denote the arm
directions of the triangular-shaped interferometers (we recall that
$\alpha=X,\,Y,\,Z$, while $O=A,\,E,\,T$, and the matrix $c$ is given in eq.
(4.7).

In eq. (5.10) we also defined

$\kappa_{ij}\equiv 2\pi
f\left|\Delta\vec{x}_{ij}\right|\;\;\;,\;\;\;\hat{s}_{ij}\equiv\frac{\Delta\vec{x}_{ij}}{\left|\Delta\vec{x}_{ij}\right|}\;\;,\;\;{\rm
with}\;\;\Delta\vec{x}_{ij}=\vec{x}_{i}-\vec{x}_{j}\;.$ (5.12)

The decomposition of eq. (5.10) introduces the coefficients

$\displaystyle\gamma_{\ell
m,\lambda,abcd}(\kappa,\hat{s})\equiv\frac{5}{8\pi}\int d^{2}{\hat{n}}{\rm
e}^{i\kappa{\hat{n}}\cdot\hat{s}}\,Y_{\ell
m}\left({\hat{n}}\right)e_{ab}^{*\lambda}\left({\hat{n}}\right)e_{cd}^{\lambda}\left({\hat{n}}\right)\;,$
(5.13)

that can be evaluated just once for all for any multipole of interest,
independently of the pair of detectors under consideration. These coefficients
can be evaluated numerically or analytically. In appendix C, we provide the
explicit analytic expressions for the monopole, dipole, and quadrupole
contributions101010The expressions for the monopole and dipole can be found in
[29], and we report them here for completeness, while the expression for the
quadrupole is an original result of this paper. The explicit analytic
expressions for the non-vanishing coefficients in the zero frequency limit can
be found in [28]., and we explain how this computation can be extended to
higher multipoles.

###  5.2 Covariance matrix and optimal $\chi^{2}$

In appendix B, we evaluate the covariance matrix from the first of eqs. (5)
for an arbitrary level of signal and noise. We show that the covariance
assumes the structure (indices $OO^{\prime}$ accounting for the channels of
each triangular-shaped detector are omitted for brevity)

$\displaystyle\Sigma_{ij,kl,m^{\prime}m^{\prime\prime}}^{2}\simeq\delta_{m^{\prime}m^{\prime\prime}}(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk})\int_{-\infty}^{+\infty}dfV_{ij}(f)\left|Q_{ij,m^{\prime}}(f)\right|^{2}\;,$
(5.14)

where

$\displaystyle V_{ij}(f)$
$\displaystyle=\left(\frac{8\pi}{5}\right)^{2}\frac{\tau^{2}\delta_{00,GW}^{2}}{T}\Bigg{[}\sum_{\lambda}H_{\lambda}(f)\gamma_{ii,\lambda}(f)\sum_{\lambda^{\prime}}H_{\lambda^{\prime}}(f)\gamma_{jj,\lambda^{\prime}}(f)$
$\displaystyle+\sum_{\lambda}H_{\lambda}(f)\gamma_{ij,\lambda}(f)\sum_{\lambda^{\prime}}H_{\lambda^{\prime}}(f)\gamma_{ji,\lambda^{\prime}}(f)$
$\displaystyle+\frac{1}{\sqrt{4\pi}}\frac{5}{16\pi}\left(\sum_{\lambda}H_{\lambda}(f)\gamma_{ii,\lambda}(f)P_{j}(f)+\sum_{\lambda}H_{\lambda}(f)\gamma_{jj,\lambda}(f)P_{i}(f)\right)$
$\displaystyle+\frac{1}{4\pi}\left(\frac{5}{16\pi}\right)^{2}P_{i}(f)P_{j}(f)\Bigg{]}\;,$
(5.15)

where we defined $\gamma_{ij,\lambda}(f)\equiv\gamma_{ij,00,\lambda}(f)$ taken
from eq. (5.9). The terms proportional to $\gamma_{jj,\lambda}$ in this last
expression are due to the contribution of the dominant isotropic component of
the signal, while the terms proportional to $P_{i}$ (which are the Power
Spectral Density plotted in figure 1) are due to the noise.

The tensorial structure in eq. (5.14) is exact in the noise-dominated regime
since the noise is diagonal across different instruments and channels. On the
contrary, the variance of the signal has contributions that are diagonal in
each pair of detectors, and off-diagonal contributions that originate from
products of correlators taken at different pairs of detectors. Only the former
contributions respect the tensorial structure of eq. (5.14). We verified, by a
direct inspection of about $100$ correlators in the networks that we
considered in section 2, that the latter contributions are negligible. We
refer the reader to the discussion after eq. (B.26) for more details on this.

Hereafter, we consider an anisotropic power spectrum as in eq. (3.10) such
that

$P(\hat{n})=1+\delta^{GW}_{\ell m}Y_{lm}(\hat{n})\;,$ (5.16)

with $\ell$ and $m\neq 0$. In this case the only measurements
$C_{ij,m^{\prime}}$ in eq. (5.8) with a non-vanishing expectation value are
$C_{ij,0}$ and $C_{ij,m}$. Plugging this expression, along with eqs. (5.8) and
(5.14), into the second equation of eqs. (5) one obtains the following
expression for the $\chi^{2}$:

$\displaystyle\chi^{2}$
$\displaystyle=\sum_{ij}\sum_{\\{m^{\prime},m^{\prime\prime}\\}=\\{0,m\\}}\int
df\langle
r_{ij,m^{\prime}}^{*}(f)\rangle\,Q_{ij,m^{\prime}}^{*}(f)\Sigma_{ij,m^{\prime},m^{\prime\prime}}^{-2}\int
df\langle r_{ij,m^{\prime\prime}}(f)\rangle\,Q_{ij,m^{\prime\prime}}(f)=$
$\displaystyle=\sum_{ij}\left[\frac{\left|\int df\langle
r_{ij,0}(f)\rangle\,Q_{ij,0}(f)\right|^{2}}{\int
dfV_{ij}(f)\,\left|Q_{ij,0}(f)\right|^{2}}+\frac{\left|\int df\langle
r_{ij,m}(f)\rangle\,Q_{ij,m}(f)\right|^{2}}{\int
dfV_{ij}(f)\,\left|Q_{ij,m}(f)\right|^{2}}\right]\;,$ (5.17)

where we used the fact that the covariance matrix is diagonal in the indices
$(m^{\prime},m^{\prime\prime})$, as we show in appendix B.2 and we have
defined the functions $r_{ij,m}(f)$ through

$\displaystyle r_{ij,m}=\int_{-\infty}^{\infty}dfr_{ij,m}(f)Q_{ij,m}(f)\,,$
(5.18)

with the coefficients $r_{ij,m}$ introduced in eq. (5).

The optimization of $Q_{ij,0}(f)$ and $Q_{ij,m}(f)$ (namely, the choice of the
filter functions that maximize the signal-to-noise ratio) is performed by
equating to zero the first functional derivative of the $\chi^{2}$ over
$Q_{ij,0}(f)$ and $Q_{ij,m}(f)$. Doing so, the optimal filters are found to be

$\displaystyle Q^{*}_{ij,0}(f)=\frac{\langle
r_{ij,0}(f)\rangle}{V_{ij}(f)}\;\;\;,\;\;\;Q^{*}_{ij,m}(f)=\frac{\langle
r_{ij,m}(f)\rangle}{V_{ij}(f)}\;,$ (5.19)

leading to the optimal $\chi^{2}$:

$\displaystyle\chi^{2}_{\text{opt}}=\sum_{ij}\left[\int df\frac{|\langle
r_{ij,0}(f)\rangle|^{2}}{V_{ij}(f)}+\int df\frac{|\langle
r_{ij,m}(f)\rangle|^{2}}{V_{ij}(f)}\right]\;.$ (5.20)

This relation, with the numerator from eqs. (5) and (5.8), and the denominator
from eq. (5.15) is our main analytic result. It can be employed to study an
SGWB which is both anisotropic (with the anisotropy dominated by one
multipole) and polarized. In the following two sections, we study how this
expression simplifies in cases in which only one of these two properties is
present.

##  6 Detection of anisotropies

As the first application of the formalism developed in the previous section,
we quantify the ability to detect anisotropies in the SGWB with a network of
ground-based detectors. We do so through a Fisher forecast. To achieve this,
in subsection 6.1 we write explicitly how the expectation value of the
measurements depends on the relevant parameters for the study of the
anisotropies. On the other hand, in subsection 6.3 we use the $\chi^{2}$ to
obtain a well-posed posterior function for the parameters controlling the
anisotropy in the Gaussian limit.

###  6.1 Expectation value of the signal

We firstly consider an unpolarized power spectrum by setting $p=0$ in eq.
(3.12) and assuming that the tilt of the power spectrum $\alpha$ in eq. (3.15)
has a known value. Under this assumption, from eq. (5), (5.8), (3.15), and
(5.16) we obtain

$\displaystyle\langle r_{ij,0}(f)\rangle$
$\displaystyle=\frac{3H_{0}^{2}\tau}{20\pi^{2}}\delta^{GW}_{00}(\Omega_{0}-\bar{\Omega}_{0})\left(\frac{f}{f_{0}}\right)^{\alpha}f^{-3}\gamma_{ij}(f)\;,$
$\displaystyle\langle r_{ij,m}(f)\rangle$
$\displaystyle=\frac{3H_{0}^{2}\tau}{20\pi^{2}}(\Omega_{0}\delta^{GW}_{\ell
m}-\hat{\Omega}_{0}\hat{\delta}^{GW}_{\ell
m})\left(\frac{f}{f_{0}}\right)^{\alpha}f^{-3}\gamma_{ij,lm}(f)\;,$ (6.1)

where we introduced the shorthand notation for

$\displaystyle\gamma_{ij,\ell m}(f)\equiv\sum_{\lambda}\gamma_{\ell
m,ij,\lambda}(f)\;\;,\;\;\gamma_{ij}(f)\equiv\gamma_{ij,00}(f)\;.$ (6.2)

###  6.2 Covariance matrix and optimal $\chi^{2}$

With an explicit evaluation of eq. (5.14) under the assumption of an
unpolarized background $(p=0)$, the variance matrix, evaluated for a generic
fiducial value $\bar{\Omega}_{0}$ (i.e. without assuming either the signal- or
the noise-dominated regime), reads

$\displaystyle\Sigma$
${}_{ij,kl,m^{\prime}m^{\prime\prime}}^{2}\simeq\delta_{m^{\prime}m^{\prime\prime}}(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk})\frac{\tau^{2}}{T}\left(\frac{8\pi}{5}\right)^{2}\delta^{2}_{00,GW}\int_{-\infty}^{+\infty}dfH^{2}(f)\left|Q_{ij,m^{\prime}}(f)\right|^{2}$
$\displaystyle\;\;\;\;\;\;\;\;\times\Bigg{[}\gamma_{ii}(f)\gamma_{jj}(f)+\gamma_{ij}(f)\gamma_{ji}(f)+\frac{1}{\sqrt{4\pi}}\frac{5}{16\pi}\frac{\gamma_{ii}(f)P_{j}(f)+\gamma_{jj}(f)P_{i}(f)}{H(f)}$
$\displaystyle\;\;\;\;\;\;\;\;+\frac{1}{4\pi}\left(\frac{5}{16\pi}\right)^{2}\frac{P_{i}(f)P_{j}(f)}{H^{2}(f)}\Bigg{]}\;.$
(6.3)

Setting in this expression the filters to maximize the $\chi^{2}$ as in eq.
(5.19) (see appendix B for details), we obtain

$\displaystyle\chi^{2}_{\rm opt}$
$\displaystyle=\frac{T}{1\text{year}}\left[\left(\frac{\Omega_{0}}{\bar{\Omega}_{0}}-1\right)^{2}I_{00}+\frac{1}{4\pi}\left|\frac{\Omega_{0}}{\bar{\Omega}_{0}}\delta^{GW}_{\ell
m}-\hat{\delta}^{GW}_{\ell m}\right|^{2}I_{\ell m}\right]\;,$ (6.4)

where we defined (a power law frequency dependence as in eq. (3.15) is
assumed)

$\displaystyle I_{00/\ell m}$ $\displaystyle\equiv\sum_{i\neq
j}\int_{0}^{\infty}\frac{df}{\text{Hz}}6.31\times
10^{7}\frac{|\gamma_{ij,00/\ell m}(f)|^{2}}{D_{ij}(f)}\;,$ $\displaystyle
D_{ij}(f)$
$\displaystyle\equiv\gamma_{ii}(f)\gamma_{jj}(f)+\gamma_{ij}(f)\gamma_{ji}(f)+\frac{\gamma_{ii}(f)P_{j}(f)+\gamma_{jj}(f)P_{i}(f)}{\bar{\Omega}_{0}\frac{16\pi\sqrt{4\pi}}{5}\left(\frac{f}{f_{0}}\right)^{\alpha}f^{-3}\frac{3H_{0}^{2}}{32\pi^{3}}}+$
$\displaystyle\;\;\;\;+\frac{P_{i}(f)P_{j}(f)}{\bar{\Omega}_{0}^{2}\left(\frac{16\pi\sqrt{4\pi}}{5}\right)^{2}\left(\frac{f}{f_{0}}\right)^{2\alpha}f^{-6}\left(\frac{3H_{0}^{2}}{32\pi^{3}}\right)^{2}}\;.$
(6.5)

We note that the result depends linearly on the total observation time $T$. In
contrast, the “auxiliary time” $\tau$ (the interval duration over which the
time stream is Fourier transformed) cancels out from the combined expression.
Furthermore, in this regime, the coefficients $I_{00}$ and $I_{\ell m}$ do not
depend on the specific frequency shape of $\Omega_{GW}(f)$ but solely on the
geometry of the network (i.e. the baselines between the detectors and the
orientation of their arms). This is because the anisotropy coefficients do not
depend on frequency, cf. eq. (3.8), so the frequency shape cancels out in the
ratio between the expected signal and the variance.

We also note that the final $\chi^{2}$ can be rewritten just in terms of the
relative amplitude of $\Omega_{0}$ with respect to its fiducial value
$\bar{\Omega}_{0}$, namely on the ratio

$\displaystyle\omega\equiv\frac{\Omega_{0}}{\bar{\Omega}_{0}}\;,$ (6.6)

whose fiducial value is, by definition, $\bar{\omega}=1$. As a last remark, a
direct inspection of $\gamma_{ij,\ell m}(f)$ in eq. (6.3) shows that this
$\chi^{2}$ is symmetric under the redefinition $m\to-m$.

We conclude this subsection with two expressions showing how eq. (6.2)
simplifies in the noise-dominated and the signal-dominated regime. These
expressions are shown for clarity and completeness. However, they are not used
in our explicit evaluations, where the general expression for the covariance
is used to make a Fisher forecast for detecting the anisotropies. In the
noise-dominated regime, eq. (6.2) reduces to

$\displaystyle I^{N}_{00/\ell m}\equiv\bar{\Omega}_{0}^{2}\sum_{i\neq
j}\int_{0}^{\infty}\frac{df}{\text{Hz}}1.28\times
10^{-66-4\alpha}\,\left(\frac{f}{1\,{\rm
Hz}}\right)^{2\alpha-6}\frac{|\gamma_{ij,00/\ell
m}(f)|^{2}}{P_{i}(f)P_{j}(f)\text{Hz}^{2}}\;.$ (6.7)

while in the signal-dominated case, we have

$\displaystyle I_{00/\ell m}^{S}=\sum_{i\neq
j}\int_{0}^{\infty}\frac{df}{\text{Hz}}6.31\times
10^{7}\frac{|\gamma_{ij,00/\ell
m}(f)|^{2}}{\left[\gamma_{ii}(f)\gamma_{jj}(f)+\gamma_{ij}(f)\gamma_{ji}(f)\right]}\;.$
(6.8)

###  6.3 Fisher forecast

In eq. (6.4) we introduced the $\chi^{2}$ as a function of the parameters
($\omega$ and $\delta^{GW}_{\ell m}$) and of their fiducial values
($\bar{\omega}=1,\hat{\delta}^{GW}_{\ell m}$). We consider the null case of
$\hat{\delta}^{GW}_{\ell m}=0$ to obtain the range of measurements that is
compatible with this hypothesis at a given confidence level (C.L.) so that a
measurement outside this range can be considered as an evidence of an
anisotropic signal at that C.L.

The Fisher analysis consists of Taylor expanding the $\chi^{2}$ about an
extremum (namely, about the fiducial values) to quadratic order in the
parameters. One only needs the Fisher matrix

$\displaystyle\mathcal{F}_{\alpha\beta}\equiv\left.\frac{1}{2}\frac{\partial^{2}\chi^{2}_{\rm
opt}}{\partial\alpha\partial\beta}\right|_{\begin{subarray}{c}\alpha=\hat{\alpha}\\\
\beta=\hat{\beta}\end{subarray}}\;,$ (6.9)

which, for the set of parameters $(\omega,\delta^{GW}_{\ell m})$ of our
interest, evaluates to

$\displaystyle\mathcal{F}=\begin{pmatrix}\mathcal{F}_{\omega\omega}&\mathcal{F}_{\omega\,\delta_{\ell
m}^{GW}}\\\ \mathcal{F}_{\delta_{\ell m}^{GW}\omega}&\mathcal{F}_{\delta_{\ell
m}^{GW}\delta_{\ell
m}^{GW}}\end{pmatrix}=\frac{T}{1\text{year}}\begin{pmatrix}I_{00}&0\\\
0&\frac{1}{4\pi}I_{\ell m}\end{pmatrix}\;.$ (6.10)

Correspondingly, the expanded $\chi^{2}$ is

$\displaystyle\chi^{2}_{F}(\omega,\delta^{GW}_{\ell m})$
$\displaystyle=\begin{pmatrix}\omega-1\\\ \delta^{GW}_{\ell
m}-\hat{\delta}^{GW}_{\ell
m}\end{pmatrix}^{T}\cdot\begin{pmatrix}\mathcal{F}_{\omega\omega}&\mathcal{F}_{\omega\,\delta_{\ell
m}^{GW}}\\\ \mathcal{F}_{\delta_{\ell m}^{GW}\omega}&\mathcal{F}_{\delta_{\ell
m}^{GW}\delta_{\ell m}^{GW}}\end{pmatrix}\cdot\begin{pmatrix}\omega-1\\\
\delta^{GW}_{\ell m}-\hat{\delta}^{GW}_{\ell m}\end{pmatrix}=$
$\displaystyle=\frac{T}{1\text{year}}\left[\left(\omega-1\right)^{2}I_{00}+\frac{1}{4\pi}|\delta^{GW}_{\ell
m}|^{2}I_{\ell m}\right]\;,$ (6.11)

where the $I_{00}$ and $I_{\ell m}$ elements are given in eq. (6.2) for the
general case but can be approximated to eq. (6.7) in the noise-dominated case
and to eq. (6.8) in the signal-dominated case, depending on the specific
network of instruments used. From the $\chi^{2}$ one obtains the posterior
distribution on the parameters $(\omega,\delta^{GW}_{\ell m})$

$\displaystyle P(\omega,\delta^{GW}_{\ell
m})\propto\Pi(\omega,\delta^{GW}_{\ell
m})e^{-\frac{1}{2}\chi^{2}_{fish}(\omega,\delta^{GW}_{\ell m})}\;,$ (6.12)

where $\Pi(\omega,\delta^{GW}_{\ell m})$ is the prior function on the
parameters, and the proportionality is given up to a constant normalization
factor.

If one then wishes to marginalize over all the possible values of $\omega$,
the forecasted posterior function (assuming a flat prior) on the parameter
$\delta^{GW}_{\ell m}$ is

$\displaystyle P\left(\delta^{GW}_{\ell
m}\right)=\sqrt{\frac{\mathcal{F}_{\delta\delta}}{2\pi}}\exp(-\frac{\delta_{\ell
m,GW}^{2}}{2}\mathcal{F}_{\delta\delta})\;.$ (6.13)

Under these assumptions, we define the error of the measurement of the
anisotropic coefficient $(\ell,m)$ as $\sigma_{\delta_{\ell
m}}=\mathcal{F}_{\delta\delta}^{-1}$.

These expressions are used to obtain the numerical results of subsection 2.1.

##  7 Detection of a net polarization

In this section, we apply the above formalism to probe a net polarization of
the SGWB. Analogously to the previous section, this analysis aims to obtain a
Fisher forecast of the measurement of a net polarization with a network of
interferometers. To do so, in subsection 7.1, we derive the expectation value
of the measurement as a function of the theoretical parameters (the amplitude
of the signal and its polarization). We evaluate the covariance and Fisher
matrices in subsection 7.2. In subsection 7.3, we finally derive the explicit
form of a posterior function on the parameters used to describe the polarized
SGWB in the Gaussian limit.

###  7.1 Expectation value of the signal

We evaluate eqs. (5.8), (5.14), and (5.20) in the case of a power spectrum
described by eqs. (3.12) and eq. (3.15), for a generic value of net
polarization $p$ (ranging from $-1$ to $+1$) and assuming a completely
isotropic background, i.e. $\delta^{GW}_{\ell m}=0$ for every
$(\ell,m)\neq(0,0)$. In this case, the only coefficient
$\mathcal{C}_{ij,m}(f)$ which has a non-zero expectation value is
$\mathcal{C}_{ij,0}(f)$, therefore we are only interested in computing

$\displaystyle\langle
r_{ij,0}(f)\rangle=\frac{8\pi\tau}{5}\delta^{GW}_{00}\left(\frac{f}{f_{0}}\right)^{\alpha}f^{-3}\frac{3H_{0}^{2}}{32\pi^{3}}\left(\Omega_{0}\sum_{\lambda}\gamma_{ij,\lambda}(f)(1+p\lambda)-\bar{\Omega}_{0}\sum_{\lambda}\gamma_{ij,\lambda}(f)(1+\hat{p}\lambda)\right)\;,$
(7.1)

where $\gamma_{ij,\lambda}(f)=\gamma_{ij,00,\lambda}(f)$. Then, separating the
response functions in their parity invariant and parity-violating (introduced
in appendix C.1),

$\displaystyle\gamma_{ij,\lambda}(f)=\gamma_{ij,I}(f)+\lambda\gamma_{ij,V}(f)\;,$
(7.2)

the expectation value of the signal can then be rewritten as

$\displaystyle\langle
r_{ij,0}(f)\rangle=\frac{16\pi\tau}{5}\delta^{GW}_{00}\left(\frac{f}{f_{0}}\right)^{\alpha}f^{-3}\frac{3H_{0}^{2}}{32\pi^{3}}\left((\Omega_{0}-\bar{\Omega}_{0})\gamma_{ij,I}(f)+(p\,\Omega_{0}-\bar{p}\,\bar{\Omega}_{0})\gamma_{ij,V}(f)\right)\;.$
(7.3)

###  7.2 Covariance matrix and optimal $\chi^{2}$

The covariance evaluates to

$\displaystyle\Sigma_{ij,kl}^{2}$
$\displaystyle=(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk})\frac{\tau^{2}}{T}\left(\frac{8\pi}{5}\right)^{2}\delta^{2}_{00,GW}\int_{-\infty}^{+\infty}df\,H^{2}(f)\left|Q_{ij}(f)\right|^{2}\times$
$\displaystyle\times\Bigg{[}\gamma_{ii,I}(f)\gamma_{jj,I}(f)+\left(\gamma_{ij,I}(f)+\hat{p}\gamma_{ij,V}(f)\right)^{2}+$
$\displaystyle+\frac{1}{\sqrt{4\pi}}\frac{5}{32\pi}\frac{\gamma_{ii,I}(f)P_{j}(f)+\gamma_{jj,I}(f)P_{i}(f)}{H(f)}+\frac{1}{4\pi}\left(\frac{5}{32\pi}\right)^{2}\frac{P_{i}(f)P_{j}(f)}{H^{2}(f)}\Bigg{]}\;,$
(7.4)

since $\gamma_{ii,V}(f)=0$, and therefore

$\displaystyle\chi^{2}$ $\displaystyle=\frac{T}{1\text{year}}\sum_{i\neq
j}\int_{0}^{\infty}\frac{df}{\text{Hz}}6.31\times
10^{7}\frac{\left[(\omega-1)\gamma_{ij,I}(f)+(\omega
p-\bar{p})\gamma_{ij,V}(f)\right]^{2}}{D_{ij}(f)}\;,$ (7.5)

with

$\displaystyle D_{ij}(f)$
$\displaystyle=\Bigg{[}\gamma_{ii,I}(f)\gamma_{jj,I}(f)+\gamma_{ij,I}^{2}(f)+\bar{p}^{2}\gamma_{ij,V}^{2}(f)+2\bar{p}\gamma_{ij,I}(f)\gamma_{ij,V}(f)+$
$\displaystyle+\frac{\gamma_{ii,I}(f)P_{j}(f)+\gamma_{jj,I}(f)P_{i}(f)}{\bar{\Omega}_{0}\frac{32\pi\sqrt{4\pi}}{5}\left(\frac{f}{f_{0}}\right)^{\alpha}f^{-3}\frac{3H_{0}^{2}}{32\pi^{3}}}+\frac{P_{i}(f)P_{j}(f)}{\bar{\Omega}_{0}^{2}\left(\frac{32\pi\sqrt{4\pi}}{5}\right)^{2}\left(\frac{f}{f_{0}}\right)^{2\alpha}f^{-6}\left(\frac{3H_{0}^{2}}{32\pi^{3}}\right)^{2}}\Bigg{]}\;.$
(7.6)

The vector of parameters will be $\theta=\left(\begin{smallmatrix}\omega\\\
p\end{smallmatrix}\right)$, and we can then perform the standard Fisher
approximation

$\displaystyle\chi^{2}_{F}(\omega,p)$
$\displaystyle\simeq\begin{pmatrix}\omega-1\\\
p-\bar{p}\end{pmatrix}^{T}\mathcal{F}(\hat{p})\begin{pmatrix}\omega-1\\\
p-\bar{p}\end{pmatrix}\;\;\;,\;\;\;\mathcal{F}_{ab}\equiv\left.\frac{1}{2}\frac{\partial^{2}\chi^{2}(\omega,p)}{\partial
a\partial b}\right|_{\begin{subarray}{c}a=\bar{a}\\\
b=\bar{b}\end{subarray}}\;.$ (7.7)

The Fisher matrix evaluates to (recall that $\bar{\omega}=1$ by definition)

$\displaystyle\mathcal{F}_{\omega\omega}$
$\displaystyle=\frac{T}{1\text{year}}\sum_{i\neq
j}\int_{0}^{\infty}\frac{df}{\text{Hz}}6.31\times
10^{7}\frac{(\gamma_{ij,I}(f)+\bar{p}\gamma_{ij,V}(f))^{2}}{D_{ij}(f)}\;,$
$\displaystyle\mathcal{F}_{\omega p}$
$\displaystyle=\frac{T}{1\text{year}}\sum_{i\neq
j}\int_{0}^{\infty}\frac{df}{\text{Hz}}6.31\times
10^{7}\frac{(\gamma_{ij,I}(f)+\bar{p}\gamma_{ij,V}(f))\gamma_{12,V}(f)}{D_{ij}(f)}\;,$
$\displaystyle\mathcal{F}_{pp}$
$\displaystyle=\frac{T}{1\text{year}}\sum_{i\neq
j}\int_{0}^{\infty}\frac{df}{\text{Hz}}6.31\times
10^{7}\frac{\gamma_{ij,V}^{2}(f)}{D_{ij}(f)}\;.$ (7.8)

Analogously to what was done in the previous section, we conclude this
subsection with two parts, where we report the explicit evaluation of the
Fisher matrix of eq. (7.2) in the noise-dominated and the signal-dominated
regimes. These are presented for completeness even if the full formula in eq.
(7.2) is used to perform the numerical results in our analysis. In particular,
in the noise-dominated regime eq. (7.2) reduces to

$\displaystyle\mathcal{F}_{\omega\omega}^{N}$
$\displaystyle=\frac{T}{1\text{year}}\sum_{i\neq
j}\int_{0}^{\infty}\frac{df}{\text{Hz}}1.28\times
10^{-66-4\alpha}\,\left(\frac{f}{1\,{\rm
Hz}}\right)^{2\alpha-6}\frac{(\gamma_{ij,I}(f)+\bar{p}\gamma_{ij,V}(f))^{2}}{P_{i}(f)P_{j}(f)\text{Hz}^{2}}\;,$
$\displaystyle\mathcal{F}_{\omega p}^{N}$
$\displaystyle=\frac{T}{1\text{year}}\sum_{i\neq
j}\int_{0}^{\infty}\frac{df}{\text{Hz}}1.28\times
10^{-66-4\alpha}\,\left(\frac{f}{1\,{\rm
Hz}}\right)^{2\alpha-6}\frac{(\gamma_{ij,I}(f)+\bar{p}\gamma_{ij,V}(f))\gamma_{12,V}(f)}{P_{i}(f)P_{j}(f)\text{Hz}^{2}}\;,$
$\displaystyle\mathcal{F}_{pp}^{N}$
$\displaystyle=\frac{T}{1\text{year}}\sum_{i\neq
j}\int_{0}^{\infty}\frac{df}{\text{Hz}}1.28\times
10^{-66-4\alpha}\,\left(\frac{f}{1\,{\rm
Hz}}\right)^{2\alpha-6}\frac{\gamma_{ij,V}^{2}(f)}{P_{i}(f)P_{j}(f)\text{Hz}^{2}}\,,$
(7.9)

while in the signal-dominated case, we have

$\displaystyle\mathcal{F}_{\omega\omega}^{S}$
$\displaystyle=\frac{T}{1\text{year}}\sum_{i\neq
j}\int_{0}^{\infty}\frac{df}{\text{Hz}}6.31\times
10^{7}\frac{(\gamma_{ij,I}(f)+\bar{p}\gamma_{ij,V}(f))^{2}}{\gamma_{ii,I}(f)\gamma_{jj,I}(f)+\left(\gamma_{ij,I}(f)+\bar{p}\gamma_{ij,V}(f)\right)^{2}}\;,$
$\displaystyle\mathcal{F}_{\omega p}^{S}$
$\displaystyle=\frac{T}{1\text{year}}\sum_{i\neq
j}\int_{0}^{\infty}\frac{df}{\text{Hz}}6.31\times
10^{7}\frac{(\gamma_{ij,I}(f)+\bar{p}\gamma_{ij,V}(f))\gamma_{12,V}(f)}{\gamma_{ii,I}(f)\gamma_{jj,I}(f)+\left(\gamma_{ij,I}(f)+\bar{p}\gamma_{ij,V}(f)\right)^{2}}\;,$
$\displaystyle\mathcal{F}_{pp}^{S}$
$\displaystyle=\frac{T}{1\text{year}}\sum_{i\neq
j}\int_{0}^{\infty}\frac{df}{\text{Hz}}6.31\times
10^{7}\frac{\gamma_{ij,V}^{2}(f)}{\gamma_{ii,I}(f)\gamma_{jj,I}(f)+\left(\gamma_{ij,I}(f)+\bar{p}\gamma_{ij,V}(f)\right)^{2}}\;.$
(7.10)

###  7.3 Fisher forecast

Finally, we can compute the posterior function of the parameters $(\omega,p)$,
obtaining an expression analogous to eq. (6.12)

$\displaystyle
P(\omega,p)\propto\Pi(\omega,p)e^{-\frac{1}{2}\chi^{2}_{F}(\omega,p)}\;.$
(7.11)

Under the assumption of a uniform prior on the parameter $p$, after
marginalizing over all the possible values of $\omega$, the posterior
distribution on the single parameter p will be

$\displaystyle
P(p)\propto\exp(-\frac{p^{2}}{2}\frac{\mathcal{F}_{\omega\omega}\mathcal{F}_{pp}-\mathcal{F}_{\omega
p}^{2}}{\mathcal{F}_{\omega\omega}})\;.$ (7.12)

Note that since the Fisher matrix depends on $\bar{p}$, the posterior
functions do as well. As a last remark, we note that thanks to the Schwarz
inequality, we have

$\displaystyle\det(F)=\mathcal{F}_{\omega\omega}\mathcal{F}_{pp}-\mathcal{F}_{p\omega}^{2}\geq
0\;,$ (7.13)

which means that the Fisher matrix is always definite positive, and therefore,
the quadratic approximation for the $\chi^{2}$ is always well defined. In this
case, the error in the measurement of the polarization factor $p$,
marginalized over the possible values of the amplitude of the logarithmic
energy density $\Omega_{0}$, is given by

$\displaystyle\sigma_{p}=\left(\frac{\mathcal{F}_{\omega\omega}\mathcal{F}_{pp}-\mathcal{F}_{\omega
p}^{2}}{\mathcal{F}_{\omega\omega}}\right)^{-\frac{1}{2}}\;.$ (7.14)

These expressions are used to obtain the numerical results of subsection 2.2.

## Acknowledgments

We thank Vuk Mandic and Angelo Ricciardone for useful discussions. G.M.
acknowledges support from the Imperial College London Schrödinger Scholarship
scheme. C.R.C. acknowledges support under a UKRI Consolidated Grant
ST/T000791/1. M.P. is supported by Istituto Nazionale di Fisica Nucleare
(INFN) through the Theoretical Astroparticle Physics (TAsP) and the Inflation,
Dark Matter and the Large-Scale Structure of the Universe (InDark) project.

##  Appendix A Signal expectation value

As explained in section 5, we work under the assumption that the measurements
at different sites are noise-uncorrelated. For this reason, the expectation
value of the estimator ${\cal C}_{ij,m}$ in eq. (5.5) (the ones for ${\cal
C}_{ij,O,m}$ and ${\cal C}_{ij,OO^{\prime},m}$ are defined analogously)
contains the contribution coming from the correlator of the signal component
in the data stream

$\displaystyle\left\langle{\cal
C}_{ij,m}\right\rangle=\frac{1}{T}\int_{0}^{T}dt\,{\rm
e}^{-im\omega_{e}t}\,\int_{-\infty}^{+\infty}df\left\langle{\tilde{s}}_{i}^{*}\left(f,\,t\right){\tilde{s}}_{j}\left(f,\,t\right)\right\rangle{\tilde{Q}}_{ij}(f)\;,$
(A.1)

between two different sites, $i\neq j$.

To show explicitly how this component depends on the amplitude of the SGWB, we
recall that, to first order in the GW perturbation, a light beam starting from
$\vec{x}$ at the unperturbed time $t-2L$, arriving at $\vec{x}+L\,{\hat{l}}$,
and returning back to $\vec{x}$ at the unperturbed time $t$ covers the
distance in a time

$T_{\rm
return}=2L+\frac{{\hat{l}}^{a}\,{\hat{l}}^{b}}{2}\int_{0}^{L}ds\,h_{ab}\left(t-2L+s,\,\vec{x}+s\,{\hat{l}}\right)+\frac{{\hat{l}}^{a}\,{\hat{l}}^{b}}{2}\int_{0}^{L}ds\,h_{ab}\left(t-L+s,\,\vec{x}+L\,{\hat{l}}-s\,{\hat{l}}\right)\;,$
(A.2)

where $T_{0}=2L$ is the return travel time in unperturbed Minkowski spacetime.
Since our ground-based interferometers (existing and planned) are in the short
arm regime, namely $2\pi\,f\,L\ll 1$, we can approximate the GW appearing in
eq. (A.2) as $h_{ab}\left(t,\,\vec{x}\right)$, which is constant along the
integral, and therefore

$T_{\rm ret}=2L+\frac{{\hat{l}}^{a}\,{\hat{l}}^{b}}{2}\times
2\,L\,h_{ab}\left(t,\,\vec{x}\right)\;.$ (A.3)

Here we denote by $\vec{x}_{i\alpha}\left(t\right)$ the location of the vertex
$\alpha$ of the $i-$th triangular-shaped detector and by
$\vec{x}_{j}\left(t\right)$ the location of the $j-$th L-shaped detector of
the network. In a similar manner, ${\hat{u}}_{i\alpha}\left(t\right)$ and
${\hat{v}}_{i\alpha}\left(t\right)$ are the unit vectors in the directions of
the two arms that start from this vertex of each triangular-shaped detector,
as defined in [28]. These quantities are time-dependent due to the rotation of
the Earth about its axis. Under these conventions, the signal at the time $t$
at the vertex $\alpha$ of the $i$-th triangular-shaped detector, located at
$\vec{x}_{i\alpha}\left(t\right)$, is

$s_{i\alpha}(t)=d_{i\alpha}^{ab}(t)h_{ab}(t,\vec{x}_{i}(t))\;\;,\;\;d_{i\alpha}^{ab}(t)\equiv\frac{\hat{u}_{i\alpha}^{a}(t)\hat{u}_{i\alpha}^{b}(t)-\hat{v}_{i\alpha}^{a}(t)\hat{v}_{i\alpha}^{b}(t)}{2}\;.$
(A.4)

Analogously, since ${\hat{u}}_{j\alpha}\left(t\right)$ and
${\hat{v}}_{j\alpha}\left(t\right)$ are the unit vectors in the directions of
the two arms that start from the $j-$th L-shaped instrument, the signal at the
time $t$ at the $j$-th L-shaped detector, located at
$\vec{x}_{j}\left(t\right)$, is

$s_{j}(t)=d_{j}^{ab}(t)h_{ab}(t,\vec{x}_{i}(t))\;\;,\;\;d_{j}^{ab}(t)\equiv\frac{\hat{u}_{j}^{a}(t)\hat{u}_{j}^{b}(t)-\hat{v}_{j}^{a}(t)\hat{v}_{j}^{b}(t)}{2}\;,$
(A.5)

With eqs. (A.4), (A.5) and (3.6) one finds (using a compact notation as in eq.
(5.4)

$\displaystyle{\tilde{s}}_{i(\alpha)}\left(f,\,t\right)$
$\displaystyle=\sum_{\lambda=R,L}\int
d^{2}{\hat{n}}\int_{-\infty}^{+\infty}df^{\prime}{\rm e}^{-2\pi
i\left(f-f^{\prime}\right)t}\delta_{\tau}\left(f-f^{\prime}\right)$
$\displaystyle\times{\rm e}^{-2\pi
if^{\prime}{\hat{n}}\cdot\vec{x}_{i(\alpha)}\left(t\right)}h_{\lambda}\left(f^{\prime},\,{\hat{n}}\right)e_{ab}^{\lambda}\left({\hat{n}}\right)d_{i(\alpha)}^{ab}\left(t\right)\;,$
(A.6)

where the interferometer location and arm directions can be treated as
constant in the time interval $\tau$, provided that its value is sufficiently
small. In analogy with [27] we have introduced the notation

$\delta_{\tau}\left(f\right)\equiv\frac{\sin\left(\pi\tau f\right)}{\pi
f}\;\;\;\;,\;\;\;\;\lim_{\tau\to\infty}\delta_{\tau}\left(f\right)=\delta_{D}\left(f\right)\;.$
(A.7)

We can now compute the correlators

$\displaystyle\left\langle{\tilde{s}}_{i(\alpha)}^{*}\left(f,t\right){\tilde{s}}_{j(\beta)}\left(f^{\prime},t\right)\right\rangle=$
$\displaystyle\sum_{\lambda=R,L}\int_{-\infty}^{+\infty}df^{\prime}\,\delta_{\tau}^{2}\left(f-f^{\prime}\right)H_{\lambda}\left(\left|f^{\prime}\right|\right)\int
d^{2}{\hat{n}}{\rm e}^{2\pi
if^{\prime}{\hat{n}}\cdot\left(\vec{x}_{i(\alpha)}\left(t\right)-\vec{x}_{j(\beta)}\left(t\right)\right)}$
$\displaystyle\times\sum_{\ell m}\delta^{GW}_{\ell m,\lambda}Y_{\ell
m}\left({\hat{n}}\right)\,e_{ab}^{*\lambda}\left({\hat{n}}\right)e_{cd}^{\lambda}\left({\hat{n}}\right)d_{i(\alpha)}^{ab}\left(t\right)d_{j(\beta)}^{cd}\left(t\right)\;.$
(A.8)

Since the integration time $\tau$ is set to be much greater than the inverse
of the typical measured GW frequencies, one of the two $\delta_{\tau}$ in this
expression can be replaced by a Dirac $\delta-$function, while the other one
evaluates to the integration time. We obtain

$\displaystyle\left\langle{\tilde{s}}_{i(\alpha)}^{*}\left(f,\,t\right){\tilde{s}}_{j(\beta)}\left(f,\,t\right)\right\rangle=$
$\displaystyle\tau\sum_{\lambda=R,L}H_{\lambda}\left(\left|f\right|\right)\int
d^{2}{\hat{n}}\,{\rm e}^{2\pi
if{\hat{n}}\cdot\left(\vec{x}_{i(\alpha)}\left(t\right)-\vec{x}_{j(\beta)}\left(t\right)\right)}$
$\displaystyle\times\sum_{\ell m}\delta^{GW}_{\ell m,\lambda}Y_{\ell
m}\left({\hat{n}}\right)\,e_{ab}^{*\lambda}\left({\hat{n}}\right)e_{cd}^{\lambda}\left({\hat{n}}\right)d_{i(\alpha)}^{ab}\left(t\right)d_{j(\beta)}^{cd}\left(t\right)\;.$
(A.9)

At this point, we can go from the correlation of the signal at the three
vertices to the correlator of the signal in the three channels and write the
first equation of eqs. (A.9) as

$\displaystyle\left\langle{\tilde{s}}_{iO}^{*}\left(f,\,t\right){\tilde{s}}_{jO^{\prime}}\left(f,\,t\right)\right\rangle=$
$\displaystyle\tau\sum_{\lambda=R,L}H_{\lambda}\left(\left|f\right|\right)\int
d^{2}{\hat{n}}{\rm e}^{2\pi
if{\hat{n}}\cdot\left(\vec{x}_{i}\left(t\right)-\vec{x}_{j}\left(t\right)\right)}$
$\displaystyle\times$ $\displaystyle\sum_{\ell m}\delta^{GW}_{\ell
m,\lambda}Y_{\ell
m}\left({\hat{n}}\right)e_{ab}^{*\lambda}\left({\hat{n}}\right)e_{cd}^{\lambda}\left({\hat{n}}\right)d_{iO}^{ab}\left(t\right)d_{jO^{\prime}}^{cd}\left(t\right)\;,$
(A.10)

where $d_{iO}^{ab}\equiv c_{O\alpha}d_{i\alpha}^{ab}$. By evaluating these
linear combinations, we find

$d_{iA}^{ab}=\frac{2}{\sqrt{3}}\,d_{iX}^{ab}\;\;\;\;\;,\;\;\;\;\;d_{iE}^{ab}=-\frac{2}{3}\left(d_{iX}^{ab}+2d_{iY}^{ab}\right)\;\;\;\;\;,\;\;\;\;\;d_{iT}^{ab}=0\;,$
(A.11)

namely, only the two channels $A$ and $E$ are non-vanishing in the short-arm
approximation. Inserting all this in eq. (A.1), we obtain (again trivial
extension for $\langle{\cal C}_{ij,O,m}\rangle$ and $\langle{\cal
C}_{ij,OO^{\prime},m}\rangle$)

$\displaystyle\left\langle C_{ij,m}\right\rangle=$
$\displaystyle\frac{\tau}{T}\int_{0}^{T}dt\;{\rm
e}^{-im\omega_{e}t}\sum_{\lambda=R,L}\int_{-\infty}^{+\infty}df\,H_{\lambda}\left(|f|\right)\int
d^{2}{\hat{n}}\,e_{ab}^{*\lambda}\left({\hat{n}}\right)e_{cd}^{\lambda}\left({\hat{n}}\right)\,\sum_{\ell
m^{\prime}}p_{\ell m^{\prime},\lambda}\,Y_{\ell
m^{\prime}}\left({\hat{n}}\right)$ $\displaystyle\times$ $\displaystyle
d^{ab}_{i}\left(t\right)\,d^{cd}_{j}\left(t\right)Q_{ij}\left(f\right){\rm
e}^{2\pi
if{\hat{n}}\cdot\left(\vec{x}_{i}\left(t\right)-\vec{x}_{j}\left(t\right)\right)}\;.$
(A.12)

The quantities $d^{ab}_{i,A/E}\left(t\right)$ and $d^{ab}_{i}\left(t\right)$
are time-dependent because of the rotation of the Earth in the frame of the
fixed stars. Therefore we denote by $d^{ab}_{i,A/E}$ and $d^{ab}_{i}$ the same
quantities in a frame that is fixed with respect to the Earth,

$\displaystyle
d_{i,A/E}^{ab}\left(t\right)=R_{aa^{\prime}}\left(t\right)R_{bb^{\prime}}\left(t\right)\,d_{i,A/E}^{a^{\prime}b^{\prime}}\;,$
$\displaystyle
d_{i}^{ab}\left(t\right)=R_{aa^{\prime}}\left(t\right)R_{bb^{\prime}}\left(t\right)\,d_{i}^{a^{\prime}b^{\prime}}\;,$
(A.13)

with $R\left(t\right)$ the rotation matrix of period $T_{e}$ around the
$z-$axis. By changing integration variable ${\hat{n}}\to R{\hat{n}}$ in eq.
(A.12), we can reabsorb this rotation and take advantage of the properties of
the basis tensors. We see that the rotation matrix does not appear anywhere
apart from the argument of the spherical harmonic, where it produces $Y_{\ell
m^{\prime}}\left(R{\hat{n}}\right)={\rm e}^{im^{\prime}\omega_{e}t}Y_{\ell
m^{\prime}}\left({\hat{n}}\right)$. At this point, we define the position of
the $i$-th interferometer in this reference frame as $\vec{x}_{i}$, in such a
way that we can set, as in eq. (5.12),

$\Delta\vec{x}_{ij}=\vec{x}_{i}-\vec{x}_{j}\;.$ (A.14)

It then follows immediately that the integration in the time variable forces
$m^{\prime}=m$, and

$\displaystyle\left\langle C_{ij,m}\right\rangle$
$\displaystyle=\tau\sum_{\lambda=R,L}\int_{-\infty}^{+\infty}df\,H_{\lambda}\left(|f|\right)\int
d^{2}{\hat{n}}\,e_{ab}^{*\lambda}\left({\hat{n}}\right)e_{cd}^{\lambda}\left({\hat{n}}\right)$
$\displaystyle\times\sum_{\ell=|m|}^{\infty}\delta^{GW}_{\ell
m,\lambda}\,Y_{\ell
m}\left({\hat{n}}\right)\,d^{ab}_{i}d^{cd}_{j}Q_{ij}\left(f\right){\rm
e}^{2\pi if{\hat{n}}\cdot\Delta\vec{x}_{ij}}\;,$ (A.15)

which leads to eq. (5.8) of the main text.

##  Appendix B Covariance matrix

We compute the elements of the covariance matrix
$\Sigma_{ij,kl,m^{\prime}m^{\prime\prime}}^{2}$ as defined in eq. (5),
considering initial contributions coming from both the instrumental noise and
the variance of the signal111111As remarked in the main text, for the sake of
simplicity, we perform the computation assuming a single non-vanishing (and
subdominant) multipole in addition to the monopole. For brevity, we do not
report the details of the computation of
$\Sigma_{ij(OO^{\prime}),kl(O^{\prime\prime}O^{\prime\prime\prime}),m^{\prime}m^{\prime\prime}}^{2}$,
where triangular-shaped instruments are involved, as it follows analogous
steps to the computation presented here.:

$\displaystyle\Sigma$
${}_{ij,kl,m^{\prime}m^{\prime\prime}}^{2}=\langle\mathcal{C}^{*}_{ij,m^{\prime}}\mathcal{C}_{kl,m^{\prime\prime}}\rangle+\langle\mathcal{C}^{*}_{ij,m^{\prime}}\rangle\langle\mathcal{C}_{kl,m^{\prime\prime}}\rangle-2\langle\mathcal{C}^{*}_{ij,m^{\prime}}\rangle\langle\mathcal{C}_{kl,m^{\prime\prime}}\rangle=$
$\displaystyle=\frac{1}{T^{2}}\int_{0}^{T}dtdt^{\prime}\mathrm{e}^{im^{\prime}\omega_{e}t}\mathrm{e}^{-im^{\prime\prime}\omega_{e}t^{\prime}}\int_{-\infty}^{\infty}df\int_{-\infty}^{\infty}df^{\prime}\Bigg{[}\left\langle
m_{i}(f,t)m_{k}^{*}(f^{\prime},t^{\prime})\right\rangle\left\langle
m_{j}^{*}(f,t)m_{l}(f^{\prime},t^{\prime})\right\rangle$
$\displaystyle\quad\quad\quad\quad\quad\quad+\left\langle
m_{i}(f,t)m_{l}(f^{\prime},t^{\prime})\right\rangle\left\langle
m_{j}^{*}(f,t)m_{k}^{*}(f^{\prime},t^{\prime})\right\rangle\Bigg{]}Q^{*}_{ij,m^{\prime}}(f)Q_{kl,m^{\prime\prime}}(f^{\prime})\;,$
(B.1)

which simplifies to

$\displaystyle\Sigma$
${}_{ij,kl,m^{\prime}m^{\prime\prime}}^{2}=\frac{1}{T^{2}}\int_{0}^{T}dtdt^{\prime}\mathrm{e}^{im^{\prime}\omega_{e}t}\mathrm{e}^{-im^{\prime\prime}\omega_{e}t^{\prime}}\int_{-\infty}^{\infty}df\int_{-\infty}^{\infty}df^{\prime}Q^{*}_{ij,m^{\prime}}(f)Q_{kl,m^{\prime\prime}}(f^{\prime})$
$\displaystyle\times\Bigg{[}\left\langle
s_{i}(f,t)s_{k}^{*}(f^{\prime},t^{\prime})\right\rangle\left\langle
s_{j}^{*}(f,t)s_{l}(f^{\prime},t^{\prime})\right\rangle+\left\langle
s_{i}(f,t)s_{l}(f^{\prime},t^{\prime})\right\rangle\left\langle
s_{j}^{*}(f,t)s_{k}^{*}(f^{\prime},t^{\prime})\right\rangle$
$\displaystyle+\left\langle
s_{i}(f,t)s_{k}^{*}(f^{\prime},t^{\prime})\right\rangle\left\langle
n_{j}^{*}(f,t)n_{l}(f^{\prime},t^{\prime})\right\rangle+\left\langle
s_{i}(f,t)s_{l}(f^{\prime},t^{\prime})\right\rangle\left\langle
n_{j}^{*}(f,t)n_{k}^{*}(f^{\prime},t^{\prime})\right\rangle$
$\displaystyle+\left\langle
n_{i}(f,t)n_{k}^{*}(f^{\prime},t^{\prime})\right\rangle\left\langle
s_{j}^{*}(f,t)s_{l}(f^{\prime},t^{\prime})\right\rangle+\left\langle
n_{i}(f,t)n_{l}(f^{\prime},t^{\prime})\right\rangle\left\langle
s_{j}^{*}(f,t)s_{k}^{*}(f^{\prime},t^{\prime})\right\rangle$
$\displaystyle+\left\langle
n_{i}(f,t)n_{k}^{*}(f^{\prime},t^{\prime})\right\rangle\left\langle
n_{j}^{*}(f,t)n_{l}(f^{\prime},t^{\prime})\right\rangle+\left\langle
n_{i}(f,t)n_{l}(f^{\prime},t^{\prime})\right\rangle\left\langle
n_{j}^{*}(f,t)n_{k}^{*}(f^{\prime},t^{\prime})\right\rangle\Bigg{]}\;.$ (B.2)

In terms of the quantities introduced in the main text, this rewrites as

$\displaystyle\Sigma_{ij,kl,m^{\prime}m^{\prime\prime}}^{2}=\left(\frac{8\pi}{5}\right)^{2}\frac{1}{T^{2}}\int_{0}^{T}dtdt^{\prime}\mathrm{e}^{im^{\prime}\omega_{e}t}\mathrm{e}^{-im^{\prime\prime}\omega_{e}t^{\prime}}\int_{-\infty}^{\infty}dfdf^{\prime}Q^{*}_{ij,m^{\prime}}(f)Q_{kl,m^{\prime\prime}}(f^{\prime})$
$\displaystyle\;\;\;\;\times\int_{-\infty}^{\infty}df_{1}df_{2}df_{3}df_{4}\delta_{\tau}(f-f_{1})\delta_{\tau}(f-f_{3})\delta(f_{1}-f_{2})\delta(f_{3}-f_{4})$
$\displaystyle\;\;\;\;\times\Bigg{[}e^{2\pi it(f-f_{1})}e^{-2\pi
it^{\prime}(f^{\prime}-f_{2})}e^{-2\pi it(f-f_{3})}e^{2\pi
it^{\prime}(f^{\prime}-f_{4})}\delta_{\tau}(f^{\prime}-f_{2})\delta_{\tau}(f^{\prime}-f_{4})\,\tilde{\cal
T}_{ik}^{*}\left(f_{1},f_{2}\right)\,\tilde{\cal
T}_{lm}\left(f_{3},f_{4}\right)$ $\displaystyle\;\;\;\;+e^{2\pi
it(f-f_{1})}e^{-2\pi it^{\prime}(-f^{\prime}-f_{2})}e^{-2\pi
it(f-f_{3})}e^{2\pi
it^{\prime}(-f^{\prime}-f_{4})}\delta_{\tau}(f^{\prime}+f_{2})\delta_{\tau}(f^{\prime}+f_{4})\,\tilde{\cal
T}_{il}^{*}\left(f_{1},f_{2}\right)\,\tilde{\cal
T}_{jk}\left(f_{3},f_{4}\right)\Bigg{]}\;,$ (B.3)

where we defined

$\displaystyle\mathcal{\tilde{T}}_{ij}\left(f,g,t,t^{\prime}\right)$
$\displaystyle\equiv\sum_{\lambda}H_{\lambda}(f)(\delta^{GW}_{00}\tilde{\gamma}_{ij,00,\lambda}(f,g,t,t^{\prime})+\delta^{GW}_{\ell
m}\;\tilde{\gamma}_{ij,\ell
m,\lambda}(f,g,t,t^{\prime}))+\frac{5}{16\pi}\delta_{ij}P_{i}(f)\;,$
$\displaystyle\tilde{\gamma}_{\ell m,ij,\lambda}(f,g,t,t^{\prime})$
$\displaystyle\equiv\frac{5}{8\pi}\int d^{2}\hat{n}e^{2\pi
i\hat{n}(f\textbf{x}_{i}(t)-g\textbf{x}_{j}(t^{\prime}))}Y_{\ell,m}(\hat{n})e^{*\lambda}_{ab}(\hat{n})e^{\lambda}_{bc}(\hat{n})d_{i}^{ab}(t)d_{j}^{cd}(t^{\prime})\;.$
(B.4)

After performing the $f_{2}$, $f_{4}$ integrals we obtain

$\displaystyle\Sigma$
${}_{ij,kl,m^{\prime}m^{\prime\prime}}^{2}=\left(\frac{8\pi}{5}\right)^{2}\frac{1}{T^{2}}\int_{0}^{T}dtdt^{\prime}\mathrm{e}^{im^{\prime}\omega_{e}t}\mathrm{e}^{-im^{\prime\prime}\omega_{e}t^{\prime}}$
$\displaystyle\times\int_{-\infty}^{+\infty}dfdf^{\prime}df_{1}df_{3}\delta_{\tau}(f-f_{1})\delta_{\tau}(f-f_{3})e^{2\pi
it(f_{3}-f_{1})}e^{-2\pi
it^{\prime}(f_{3}-f_{1})}Q^{*}_{ij,m^{\prime}}(f)Q_{kl,m^{\prime\prime}}(f^{\prime})$
$\displaystyle\times\Bigg{[}\delta_{\tau}(f^{\prime}-f_{1})\delta_{\tau}(f^{\prime}-f_{3})\tilde{\cal
T}_{ik}^{*}\left(f_{1},f_{1},t,t^{\prime}\right)\tilde{\cal
T}_{jl}\left(f_{3},f_{3},t,t^{\prime}\right)$
$\displaystyle\;\;\;\;+\delta_{\tau}(f^{\prime}+f_{1})\delta_{\tau}(f^{\prime}+f_{3})\tilde{\cal
T}_{il}^{*}\left(f_{1},f_{1},t,t^{\prime}\right)\tilde{\cal
T}_{jk}\left(f_{3},f_{3},t,t^{\prime}\right)\Bigg{]}\;,$ (B.5)

and with the change of variables $a=t-t^{\prime}$ and $b=t+t^{\prime}$, this
simplifies to

$\displaystyle\Sigma_{ij,kl,m^{\prime}m^{\prime\prime}}^{2}=\left(\frac{8\pi}{5}\right)^{2}\frac{1}{2T^{2}}\int_{0}^{2T}db\,\mathrm{e}^{i(m^{\prime}-m^{\prime\prime})\omega_{e}b}\int_{-T}^{T}da\,\mathrm{e}^{i(m^{\prime}+m^{\prime\prime})\omega_{e}a}$
$\displaystyle\times\int_{-\infty}^{+\infty}dfdf^{\prime}df_{1}df_{3}\delta_{\tau}(f-f_{1})\delta_{\tau}(f-f_{3})e^{2\pi
ia(f_{3}-f_{1})}Q^{*}_{ij,m^{\prime}}(f)Q_{kl,m^{\prime\prime}}(f^{\prime})$
$\displaystyle\times\Bigg{[}\delta_{\tau}(f^{\prime}-f_{1})\delta_{\tau}(f^{\prime}-f_{3})\mathcal{\tilde{T}}_{ik}^{*}\left(f_{1},a\right)\mathcal{\tilde{T}}_{jl}\left(f_{3},a\right)+\delta_{\tau}(f^{\prime}+f_{1})\delta_{\tau}(f^{\prime}+f_{3})\mathcal{\tilde{T}}_{il}^{*}\left(f_{1},a\right)\mathcal{\tilde{T}}_{jk}\left(f_{3},a\right)\Bigg{]}\;,$
(B.6)

where we have introduced121212This can be shown to be true after a change of
variables as a rotation $R(\theta)$ by an angle $\theta=\omega_{e}t^{\prime}$.

$\displaystyle\mathcal{\tilde{T}}_{ij}\left(f,a\right)$
$\displaystyle\equiv\sum_{\lambda}H_{\lambda}(f)(\delta^{GW}_{00}\tilde{\gamma}_{ij,00,\lambda}(f,a)+\delta^{GW}_{\ell
m}\;\tilde{\gamma}_{ij,\ell
m,\lambda}(f,a))+\frac{5}{16\pi}\delta_{ij}P_{i}(f)\;,$
$\displaystyle\tilde{\gamma}_{\ell m,ij,\lambda}(f,a)$
$\displaystyle\equiv\frac{5}{8\pi}\int d^{2}\hat{n}e^{2\pi
if\hat{n}(\textbf{x}_{i}(a)-\textbf{x}_{j})}Y_{\ell,m}(\hat{n})e^{*\lambda}_{ab}(\hat{n})e^{\lambda}_{bc}(\hat{n})d_{i}^{ab}(a)d_{j}^{cd}\;.$
(B.7)

The integration in $db$ then becomes trivial. This is because the covariance
is built by products of two correlations taken at two different times $t$ and
$t^{\prime}$ between different pairs. As we are assuming stationary statistics
of both the signal and the noise, the covariance is therefore insensitive to
the “absolute” time $\frac{b}{2}=\frac{t+t^{\prime}}{2}$, so that only the
dependence in $a=t-t^{\prime}$ is nontrivial. This leads to

$\displaystyle\Sigma_{ij,kl,m^{\prime}m^{\prime\prime}}^{2}=\delta_{m^{\prime}m^{\prime\prime}}\left(\frac{8\pi}{5}\right)^{2}\frac{1}{T}\int_{-T}^{T}da\,\mathrm{e}^{2im^{\prime}\omega_{e}a}$
$\displaystyle\times\int_{-\infty}^{+\infty}dfdf^{\prime}df_{1}df_{3}\delta_{\tau}(f-f_{1})\delta_{\tau}(f-f_{3})e^{2\pi
ia(f_{3}-f_{1})}Q^{*}_{ij,m^{\prime}}(f)Q_{kl,m^{\prime}}(f^{\prime})$
$\displaystyle\times\Bigg{[}\delta_{\tau}(f^{\prime}-f_{1})\delta_{\tau}(f^{\prime}-f_{3})\mathcal{\tilde{T}}_{ik}^{*}\left(f_{1},a\right)\mathcal{\tilde{T}}_{jl}\left(f_{3},a\right)+\delta_{\tau}(f^{\prime}+f_{1})\delta_{\tau}(f^{\prime}+f_{3})\mathcal{\tilde{T}}_{il}^{*}\left(f_{1},a\right)\mathcal{\tilde{T}}_{jk}\left(f_{3},a\right)\Bigg{]}\;,$
(B.8)

where it is now manifest that the covariance matrix is diagonal in
$m^{\prime},m^{\prime\prime}$.

The last expression is evaluated in the remainder of this appendix.
Specifically, for clarity of exposition, in subsection B.1 we evaluate it for
the simplest case of just two L-shaped instruments, to outline the procedure
adopted to perform all the integrations. In subsection B.2 we then extend this
computation to a generic network.

###  B.1 The simplest network: two L-shaped detectors

In the case of two L-shaped detectors, we can set $i=k=1$ and $j=l=2$ in eq.
(B). We recall that we are assuming small anisotropies, namely
$\delta^{GW}_{\ell m}\ll\delta^{GW}_{00}=\sqrt{4\pi}$. In this case

$\displaystyle\Sigma_{m}^{2}\equiv\left(\frac{8\pi}{5}\right)^{2}\frac{1}{T}\int_{-T}^{T}da\,\mathrm{e}^{im\omega_{e}a}\int_{-\infty}^{+\infty}dfdf^{\prime}df_{1}df_{3}\delta_{\tau}(f-f_{1})\delta_{\tau}(f-f_{3})e^{2\pi
ia(f_{3}-f_{1})}Q^{*}_{m}(f)Q_{m}(f^{\prime})$
$\displaystyle\times\Bigg{[}\delta_{\tau}(f^{\prime}-f_{1})\delta_{\tau}(f^{\prime}-f_{3})\mathcal{\tilde{T}}_{11}^{*}\left(f_{1},a\right)\mathcal{\tilde{T}}_{22}\left(f_{3},a\right)+\delta_{\tau}(f^{\prime}+f_{1})\delta_{\tau}(f^{\prime}+f_{3})\mathcal{\tilde{T}}_{12}^{*}\left(f_{1},a\right)\mathcal{\tilde{T}}_{21}\left(f_{3},a\right)\Bigg{]}\;,$
$\displaystyle\mathcal{\tilde{T}}_{ij}\left(f,a\right)\simeq\sum_{\lambda}H_{\lambda}(f)\delta^{GW}_{00}\tilde{\gamma}_{ij,\lambda}(f,a)+\frac{5}{16\pi}\delta_{ij}P_{i}(f)\;,$
(B.9)

where we used the shorthand notation
$\tilde{\gamma}_{ij,\lambda}=\tilde{\gamma}_{ij,00,\lambda}$ and
$Q_{m}(f)\equiv Q_{12,m}(f)$. Introducing the dimensionless integration
variable $u\equiv\frac{a}{T_{e}}$ and, making use of the periodicity of the
integrand over the time $T_{e}$, and of the fact that $T=nT_{e}$ where $n$ is
a large integer number, we have

$\displaystyle\Sigma_{m}^{2}$
$\displaystyle=\left(\frac{8\pi}{5}\right)^{2}\frac{1}{n}\int_{-n}^{n}du\,\int_{-\infty}^{+\infty}dfdf^{\prime}df_{1}df_{3}\delta_{\tau}(f-f_{1})\delta_{\tau}(f-f_{3})e^{2\pi
iu[m+T_{e}(f_{3}-f_{1})]}Q^{*}_{m}(f)Q_{m}(f^{\prime})$
$\displaystyle\times\Bigg{[}\delta_{\tau}(f^{\prime}-f_{1})\delta_{\tau}(f^{\prime}-f_{3})\mathcal{\tilde{T}}_{11}\left(f_{1},a\right)\mathcal{\tilde{T}}_{22}\left(f_{3},a\right)+\delta_{\tau}(f^{\prime}+f_{1})\delta_{\tau}(f^{\prime}+f_{3})\mathcal{\tilde{T}}_{12}\left(f_{1},a\right)\mathcal{\tilde{T}}_{21}\left(f_{3},a\right)\Bigg{]}\;,$
(B.10)

where we also used the fact that
$\mathcal{\tilde{T}}_{ij}^{*}\left(f,a\right)=\mathcal{\tilde{T}}_{ij}\left(f,a\right)$.
We note that the functions $H_{\lambda}(f)$ and $\tilde{\gamma}(f)$ vary on a
much slower scale than the $\delta_{\tau}(f)$ functions, which can therefore
be treated as Dirac $\delta$-functions, enforcing $f^{\prime}\to f$ in the
filter function when it multiplies the first term, and $f^{\prime}\to-f$ when
it multiplies the second term. Therefore, we can write

$\displaystyle\Sigma_{m}^{2}$
$\displaystyle=\left(\frac{8\pi}{5}\right)^{2}\frac{1}{n}\int_{-n}^{n}du\,\int_{-\infty}^{+\infty}dfdf^{\prime}df_{1}df_{3}\delta_{\tau}(f-f_{1})\delta_{\tau}(f-f_{3})e^{2\pi
iu[m+T_{e}(f_{3}-f_{1})]}$
$\displaystyle\times\Bigg{[}\delta_{\tau}(f^{\prime}-f_{1})\delta_{\tau}(f^{\prime}-f_{3})\mathcal{\tilde{T}}_{11}\left(f,u\right)\mathcal{\tilde{T}}_{22}\left(f,u\right)\left|Q_{m}(f)\right|^{2}$
$\displaystyle\;\;\;\;\;\;+\delta_{\tau}(f^{\prime}+f_{1})\delta_{\tau}(f^{\prime}+f_{3})\mathcal{\tilde{T}}_{12}\left(f,u\right)\mathcal{\tilde{T}}_{21}\left(f,u\right)Q^{*}_{m}(f)Q_{m}(-f)\Bigg{]}\;.$
(B.11)

Now let $g=f_{3}-f$, $h=f_{1}-f$ and $k=f^{\prime}-f$, noting that the second
term of the sum is symmetric in $f^{\prime}\to-f^{\prime}$, giving

$\displaystyle\Sigma_{m}^{2}$
$\displaystyle=\left(\frac{8\pi}{5}\right)^{2}\frac{1}{n}\int_{-n}^{n}du\,\int_{-\infty}^{+\infty}df\Bigg{[}\mathcal{\tilde{T}}_{11}\left(f,u\right)\mathcal{\tilde{T}}_{22}\left(f,u\right)\left|Q_{m}(f)\right|^{2}+\mathcal{\tilde{T}}_{12}\left(f,u\right)\mathcal{\tilde{T}}_{21}\left(f,u\right)Q^{*}_{m}(f)Q_{m}(-f)\Bigg{]}$
$\displaystyle\quad\quad\quad\quad\times\int_{-\infty}^{+\infty}dk\,dh\,dg\,e^{2\pi
iu[m+T_{e}(g-h)]}\delta_{\tau}(h)\delta_{\tau}(g)\delta_{\tau}(k-h)\delta_{\tau}(k-g)$
$\displaystyle=\left(\frac{8\pi}{5}\right)\frac{\delta^{2}_{00,GW}}{n}\int_{-n}^{n}du\,e^{2\pi
ium}\int_{-\infty}^{+\infty}df\Bigg{[}\mathcal{\tilde{T}}_{11}\left(f,u\right)\mathcal{\tilde{T}}_{22}\left(f,u\right)\left|Q_{m}(f)\right|^{2}$
$\displaystyle\quad\quad\quad\quad+\mathcal{\tilde{T}}_{12}\left(f,u\right)\mathcal{\tilde{T}}_{21}\left(f,u\right)Q^{*}_{m}(f)Q_{m}(-f)\Bigg{]}\int_{-\infty}^{+\infty}dk\,I(k,u)\,I(k,-u)\;,$
(B.12)

where

$\displaystyle I(k,u)\equiv\int_{-\infty}^{+\infty}dg\,e^{2\pi
iuT_{e}g}\delta_{\tau}(g)\delta_{\tau}(k-g)\;.$ (B.13)

To compute this quantity explicitly we expand the $\delta_{\tau}(f)$ and
$\sin(x)$ functions to obtain

$\displaystyle I(k,u)$ $\displaystyle=\frac{2\cos(\pi\tau k)}{(2\pi
i)^{2}}\int_{-\infty}^{+\infty}dg\,\frac{e^{2\pi
iuT_{e}g}}{g(k-g)}-\int_{-\infty}^{+\infty}dg\frac{e^{-i\pi\tau k}}{(2\pi
i)^{2}}\frac{e^{2i\pi\tau g(\frac{uTe}{\tau}+1)}}{g(k-g)}$
$\displaystyle-\int_{-\infty}^{+\infty}dg\frac{e^{i\pi\tau k}}{(2\pi
i)^{2}}\frac{e^{2i\pi\tau g(\frac{uTe}{\tau}-1)}}{g(k-g)}\equiv
I_{1}(k,u)+I_{2}(k,u)+I_{3}(k,u)\;.$ (B.14)

In the evaluation of the first term, we separate the two factors at the
denominator, $\frac{1}{g(k-g)}=\frac{1}{kg}-\frac{1}{k(g-k)}$, we write the
two terms in separate integrals, and we change variable $g\to g+k$ in the
second one, to obtain

$\displaystyle I_{1}(k,u)$ $\displaystyle=\frac{2\cos(\pi\tau k)}{k(2\pi
i)^{2}}(1-e^{2\pi iuT_{e}k})\int_{-\infty}^{+\infty}dg\,\frac{e^{2\pi
iuT_{e}g}}{g}\;.$ (B.15)

In the Euler expansion of the exponential, only the sine term provides a non-
vanishing contribution to the integral as $\frac{\cos(x)}{x}$ is odd

$\displaystyle I_{1}(k,u)$ $\displaystyle=\frac{2i\cos(\pi\tau k)}{k(2\pi
i)^{2}}(1-e^{2\pi iuT_{e}k})\int_{-\infty}^{+\infty}dg\,\frac{\sin(2\pi
uT_{e}g)}{g}\;,$ $\displaystyle=\frac{\cos(\pi\tau k)}{2\pi ik}(1-e^{2\pi
iuT_{e}k})\,\text{sgn}(u)\;,$ (B.16)

with $\text{sgn}(u)$ being the sign function. In an analogous manner, we find

$\displaystyle I_{2}(k,u)$ $\displaystyle=-\frac{1}{2}\frac{e^{-i\pi\tau
k}}{2\pi ik}(1-e^{2\pi
i\tau(\frac{uTe}{\tau}+1)k})\text{sgn}\left(\frac{uTe}{\tau}+1\right)\;,$
$\displaystyle I_{3}(k,u)$ $\displaystyle=-\frac{1}{2}\frac{e^{i\pi\tau
k}}{2\pi ik}(1-e^{2\pi
i\tau(\frac{uTe}{\tau}-1)k})\text{sgn}\left(\frac{uTe}{\tau}-1\right)\;,$
(B.17)

Adding the three terms, one can verify that

$\displaystyle\int_{-\infty}^{+\infty}dk\,I(k,u)I(k,-u)$
$\displaystyle=\int_{-\infty}^{+\infty}dk\,\frac{\sin^{2}(\pi k(\tau-
T_{e}\left|u\right|))}{\pi^{2}k^{2}}\,\theta\left(\frac{\tau}{T_{e}}-\left|u\right|\right)$
$\displaystyle=\left(\tau-
T_{e}\,\left|u\right|\right)\,\theta\left(\frac{\tau}{T_{e}}-\left|u\right|\right)\;,$
(B.18)

where $\theta$ is the Heaviside $\theta$ function. Therefore, eq. (B.1)
becomes

$\displaystyle\Sigma_{m}^{2}$
$\displaystyle=\left(\frac{8\pi}{5}\right)^{2}\frac{1}{n}\int_{-\frac{\tau}{T_{e}}}^{\frac{\tau}{T_{e}}}du\,(\tau-
T_{e}\left|u\right|)e^{2\pi ium}$
$\displaystyle\times\int_{-\infty}^{+\infty}df\Bigg{[}\mathcal{\tilde{T}}_{11}\left(f,u\right)\mathcal{\tilde{T}}_{22}\left(f,u\right)\left|Q_{m}(f)\right|^{2}+\mathcal{\tilde{T}}_{12}\left(f,u\right)\mathcal{\tilde{T}}_{21}\left(f,u\right)Q^{*}_{m}(f)Q_{m}(-f)\Bigg{]}$
$\displaystyle=\left(\frac{8\pi}{5}\right)^{2}\frac{\tau^{2}}{n\,T_{e}}\int_{-\infty}^{+\infty}df\,V_{m}(f)\;,$
(B.19)

where we used the fact that the functions $\mathcal{\tilde{T}}$ are even in
$u$, and where we introduced the quantity

$\displaystyle V_{m}(f)\equiv
2\int_{0}^{1}du^{\prime}\,\left(1-u^{\prime}\right)\cos\left(2\pi
u^{\prime}\frac{\tau}{T_{e}}m\right)$
$\displaystyle\times\Bigg{[}\mathcal{\tilde{T}}_{11}\left(f,u^{\prime}\frac{\tau}{T_{e}}\right)\mathcal{\tilde{T}}_{22}\left(f,u^{\prime}\frac{\tau}{T_{e}}\right)\left|Q_{m}(f)\right|^{2}+\mathcal{\tilde{T}}_{12}\left(f,u^{\prime}\frac{\tau}{T_{e}}\right)\mathcal{\tilde{T}}_{21}\left(f,u^{\prime}\frac{\tau}{T_{e}}\right)Q^{*}_{m}(f)Q_{m}(-f)\Bigg{]}\;.$
(B.20)

Next, we use the fact that the cosine term and the $\mathcal{\tilde{T}}$
functions vary very slowly in the domain of integration for $u^{\prime}$. As a
consequence, $V_{m}(f)$ can be well approximated by

$\displaystyle V_{m}(f)$
$\displaystyle\simeq\Bigg{[}\mathcal{\tilde{T}}_{11}\left(f,0\right)\mathcal{\tilde{T}}_{22}\left(f,0\right)\left|Q(f)\right|^{2}+\mathcal{\tilde{T}}_{12}\left(f,0\right)\mathcal{\tilde{T}}_{21}\left(f,0\right)Q^{*}_{m}(f)Q_{m}(-f)\Bigg{]}2\int_{0}^{1}du^{\prime}\,\left(1-u^{\prime}\right)=$
$\displaystyle=\Bigg{[}\mathcal{T}_{11}\left(f\right)\mathcal{T}_{22}\left(f\right)\left|Q_{m}(f)\right|^{2}+\mathcal{T}_{12}\left(f\right)\mathcal{T}_{21}\left(f\right)Q^{*}_{m}(f)Q_{m}(-f)\Bigg{]}\;,$
$\displaystyle\mathcal{T}_{ij}\left(f\right)$
$\displaystyle\equiv\sum_{\lambda}H_{\lambda}(f)\delta^{GW}_{00}\gamma_{ij,\lambda}(f)+\frac{5}{16\pi}\delta_{ij}P_{i}(f)\;,$
(B.21)

where we have used that fact that the function
${\tilde{\gamma}}_{ij,\lambda}\left(f,\,0\right)$, defined in eq. (B.7), and
evaluated at $a=0$, is the monopole response function
$\gamma_{ij,\lambda}\left(f\right)$ (which we evaluate in eq. (C)).

Finally, we can note that the filters $Q_{m}(f)$ must satisfy the condition
$Q_{m}^{*}(-f)=Q_{m}(f)$ so that the correlators in eq. (5) are real.
Therefore, we can write

$\displaystyle\Sigma_{m}^{2}$
$\displaystyle=\left(\frac{8\pi}{5}\right)^{2}\frac{\tau^{2}}{T}\int_{-\infty}^{+\infty}df\,\left|Q_{m}(f)\right|^{2}\Bigg{[}\mathcal{T}_{11}\left(f\right)\mathcal{T}_{22}\left(f\right)+\mathcal{T}_{12}\left(f\right)\mathcal{T}_{21}\left(f\right)\Bigg{]}\;.$
(B.22)

where $T=nT_{e}$ is the total time of observation.

###  B.2 The most general network

Extending the above computation to a generic network of $N$ instruments
(labelled by $i=1\dots N$) is straightforward. The only nontrivial extension
is related to the definition of optimal filters. Proceeding as in B.1, we
arrive to the expression

$\displaystyle\Sigma_{ij,kl,m}^{2}$
$\displaystyle=\left(\frac{8\pi}{5}\right)^{2}\frac{\tau^{2}}{T}\int_{-\infty}^{+\infty}df\,Q^{*}_{ij,m}(f)Q_{kl,m}(f)\Bigg{[}\mathcal{T}_{ik}\left(f\right)\mathcal{T}_{jl}\left(f\right)+\mathcal{T}_{il}\left(f\right)\mathcal{T}_{jk}\left(f\right)\Bigg{]}$
$\displaystyle=\int_{-\infty}^{+\infty}dfV^{\prime}_{ij,kl,m}Q^{*}_{ij,m}(f)Q_{kl,m}(f)\;,$
(B.23)

where we have introduced the quantities

$\displaystyle V_{ij,kl,m}^{\prime}(f)$
$\displaystyle\equiv\left(\frac{8\pi}{5}\right)^{2}\frac{\tau^{2}}{T}\int_{-\infty}^{+\infty}df\Bigg{[}\mathcal{T}_{ik}\left(f\right)\mathcal{T}_{jl}\left(f\right)+\mathcal{T}_{il}\left(f\right)\mathcal{T}_{jk}\left(f\right)\Bigg{]}\;,$
$\displaystyle\mathcal{T}_{ij}\left(f\right)$
$\displaystyle\equiv\sum_{\lambda}H_{\lambda}(f)\delta^{GW}_{00}\gamma_{ij,\lambda}(f)+\frac{5}{16\pi}\delta_{ij}P_{i}(f)\;.$
(B.24)

We use the shorthand notation to write all the possible pairs $a\equiv(i,j)$
(where $a$ ranges from $1$ to the combinatorial factor $\binom{N}{2}$). We
study the $\chi^{2}$ (which is diagonal in $m$ and therefore we can study the
problem at a fixed value of $m$). With this notation, the $\chi^{2}$ is

$\chi^{2}=\sum_{m^{\prime}=\\{0,m\\}}\sum_{a,b}\int
df\,r_{a,m^{\prime}}^{*}(f)\,Q_{a,m^{\prime}}^{*}(f)\;\Sigma_{ab,m^{\prime}}^{-2}\;\int
df\,r_{b,m^{\prime}}(f)\,Q_{b,m^{\prime}}(f)\;,$ (B.25)

where we recall that
$r_{a,m}\equiv\mathcal{C}_{a,m}-\left\langle\mathcal{C}_{a,m}\right\rangle$,
where the expectation value is obtained with the fiducial parameters, see eq.
(5). Moreover,

$\displaystyle\Sigma^{2}_{ab,m^{\prime}}\equiv\begin{pmatrix}\int
dfV_{aa}^{\prime}(f)\,\left|Q_{a,m^{\prime}}(f)\right|^{2}&\int
dfV_{ab}^{\prime}(f)\,Q_{a,m^{\prime}}(f)Q_{b,m^{\prime}}^{*}(f)\\\ \int
dfV_{ba}^{\prime}(f)\,Q_{a,m^{\prime}}^{*}(f)Q_{b,m^{\prime}}(f)&\int
dfV_{bb}^{\prime}(f)\,\left|Q_{b,m^{\prime}}(f)\right|^{2}\end{pmatrix}\;.$
(B.26)

To find the optimal filters, we set the functional derivatives of the
$\chi^{2}$ with respect to the filters $Q_{a,m^{\prime}}(f)$ to zero, and we
solve the corresponding linear system. We verified through inspection of
$V^{\prime}_{ab}(f)$ and $V^{\prime}_{ba}(f)$ that the off-diagonal terms of
the covariance matrix are much smaller than the diagonal ones. Therefore, we
can approximately invert eq. (B.26) to linear order in the off-diagonal
$V_{ab}^{\prime}(f)$ (with $a\neq b$) as

$\displaystyle\Sigma_{ab,m^{\prime}}^{-2}\equiv\begin{pmatrix}\frac{1}{\int
dfV_{aa}^{\prime}(f)\,\left|Q_{a,m^{\prime}}(f)\right|^{2}}&-\frac{\int
dfV_{ab}^{\prime}(f)\,Q_{a,m^{\prime}}(f)Q_{b,m^{\prime}}^{*}(f)}{\int
dfV_{aa}^{\prime}(f)\,\left|Q_{a,m^{\prime}}(f)\right|^{2}\int
dfV_{bb}^{\prime}(f)\,\left|Q_{b,m^{\prime}}(f)\right|^{2}}\\\ -\frac{\int
dfV_{ba}^{\prime}(f)\,Q_{a,m^{\prime}}^{*}(f)Q_{b,m^{\prime}}(f)}{\int
dfV_{aa}^{\prime}(f)\,\left|Q_{a,m^{\prime}}(f)\right|^{2}\int
dfV_{bb}^{\prime}(f)\,\left|Q_{b,m^{\prime}}(f)\right|^{2}}&\frac{1}{\int
dfV_{bb}^{\prime}(f)\,\left|Q_{b,m^{\prime}}(f)\right|^{2}}\end{pmatrix}\;,$
(B.27)

where we note that repeated indices such as in $V_{aa}^{\prime}(f)$ do not
imply a summation. We also note that
$V_{ba}^{\prime}(f)=V_{ab}^{*^{\prime}}(f)$ and therefore the covariance
matrix is hermitian.

Inserting this expression into eq. (B.25), and setting to zero the functional
derivatives, $\frac{\delta\chi^{2}}{\delta Q_{c,m^{\prime\prime}}}=0$, we
obtain

$\displaystyle
Q_{a,m^{\prime}}^{*}(f)=\frac{r_{a,m^{\prime}}(f)}{V_{aa}^{\prime}(f)}-\sum_{b\neq
a}\left[\frac{V_{ab}^{\prime}(f)\,r_{b,m^{\prime}}(f)}{V_{aa}^{\prime}(f)V_{bb}^{\prime}(f)}+{\cal
O}\left(V_{ab}^{{}^{\prime}2}\left(f\right)\right)\right]\;,$ (B.28)

Inserting this back into the $\chi^{2}$, we obtain the optimal chi-squared

$\chi^{2}_{\rm opt}=\sum_{a,m^{\prime}}\left\\{\int
df\frac{\left|r_{a,m^{\prime}}(f)\right|^{2}}{V^{\prime}_{aa}(f)}-\sum_{b\neq
a}\int
df\left[\frac{V^{\prime}_{ab}(f)}{V^{\prime}_{aa}(f)V^{\prime}_{bb}(f)}r^{*}_{a,m^{\prime}}(f)r_{b,m^{\prime}}(f)+{\cal
O}\left(V_{ab}^{{}^{\prime}2}\left(f\right)\right)\right]\right\\}\;.$ (B.29)

A numerical inspection of this equation shows that the contributions from the
linear terms in the off-diagonal $V^{\prime}_{ab}(f)$ are $O(10^{-4})$
suppressed with respect to the first ($b-$independent) integral, and therefore
we can write, with high accuracy

$\chi^{2}_{\rm
opt}\simeq\sum_{a,m^{\prime}}\int_{-\infty}^{+\infty}df\frac{\left|r_{a,m^{\prime}}(f)\right|^{2}}{V^{\prime}_{aa}(f)}\;.$
(B.30)

This shows that the covariance matrix in eq. (B.2) can be well approximated by
only its diagonal part in the $(i,j),(k,l)$ pairs, namely

$\displaystyle\Sigma_{ij,ij,m^{\prime}}^{2}=\left(\frac{8\pi}{5}\right)^{2}\frac{\tau^{2}}{T}\int_{-\infty}^{+\infty}df\,\left|Q_{ij,m^{\prime}}(f)\right|^{2}\Bigg{[}\mathcal{T}_{ii}\left(f\right)\mathcal{T}_{jj}\left(f\right)+\mathcal{T}_{ij}\left(f\right)\mathcal{T}_{ji}\left(f\right)\Bigg{]}\;,$
(B.31)

leading to eq. (5.14) of the main text.

##  Appendix C Evaluation of the response functions

In this appendix, we provide exact analytic expressions for the detector-
independent response functions of eq. (5.13) for the monopole, dipole, and
quadrupole. Although the response functions can be evaluated numerically,
employing the explicit analytical expressions considerably speed up the
computation of the correlators of eq. (5.8).

The expressions for the monopole and dipole can be found in [29], and we
report them here for completeness, while the expression for the quadrupole is
an original result of this paper. The extension to higher multipoles can be
performed proceeding as we do in eqs. (C.1), (C.3), and (C.10). From the
expression of the spherical harmonics, we can write

$\displaystyle\gamma_{00,abcd,\lambda}(\kappa,\hat{s})=\frac{5}{2\sqrt{\pi}}\Gamma_{abcd,\lambda}^{M}(\kappa,\hat{s})\;,$
$\displaystyle\gamma_{10,abcd,\lambda}(\kappa,\hat{s})=-\frac{5i}{2}\sqrt{\frac{3}{\pi}}\Gamma_{abcd,\lambda}^{D}(\kappa,\hat{s},\hat{z})\;,$
$\displaystyle\gamma_{1\pm
1,abcd,\lambda}(\kappa,\hat{s})=\frac{5i}{2}\sqrt{\frac{3}{2\pi}}\left[\pm\Gamma_{abcd,\lambda}^{D}(\kappa,\hat{s},\hat{x})+i\Gamma_{abcd,\lambda}^{D}(\kappa,\hat{s},\hat{y})\right]\;,$
$\displaystyle\gamma_{20,abcd,\lambda}(\kappa,\hat{s})=\frac{5}{4}\sqrt{\frac{5}{\pi}}\left[-3\Gamma_{abcd,\lambda}^{Q}(\kappa,\hat{s},\hat{z},\hat{z})-\Gamma_{abcd,\lambda}^{M}(\kappa,\hat{s})\right]\;,$
$\displaystyle\gamma_{2\pm
1,abcd,\lambda}(\kappa,\hat{s})=\frac{5}{2}\sqrt{\frac{15}{2\pi}}\left[\pm\Gamma_{abcd,\lambda}^{Q}(\kappa,\hat{s},\hat{x},\hat{z})+i\Gamma_{abcd,\lambda}^{Q}(\kappa,\hat{s},\hat{y},\hat{z})\right]\;,$
$\displaystyle\gamma_{2\pm
2,abcd,\lambda}(\kappa,\hat{s})=\frac{5}{4}\sqrt{\frac{15}{2\pi}}\left[-\Gamma_{abcd,\lambda}^{Q}(\kappa,\hat{s},\hat{x},\hat{x})\mp
2i\Gamma_{abcd,\lambda}^{Q}(\kappa,\hat{s},\hat{x},\hat{y})+\Gamma_{abcd,\lambda}^{Q}(\kappa,\hat{s},\hat{y},\hat{y})\right]\;,$
(C.1)

where we introduced the following functions:

$\displaystyle\Gamma_{abcd,\lambda}^{M}(\kappa,\hat{s})=\frac{1}{8\pi}\int
d^{2}{\hat{n}}\,{\rm
e}^{i\kappa{\hat{n}}\cdot\hat{s}}\,e_{ab}^{*\lambda}\left({\hat{n}}\right)e_{cd}^{\lambda}\left({\hat{n}}\right)\;,$
$\displaystyle\Gamma_{abcd,\lambda}^{D}(\kappa,\hat{s},\hat{v})=\frac{1}{8\pi}\int
d^{2}{\hat{n}}\,{\rm
e}^{i\kappa{\hat{n}}\cdot\hat{s}}\,e_{ab}^{*\lambda}\left({\hat{n}}\right)e_{cd}^{\lambda}\left({\hat{n}}\right)\,(i\hat{v}\cdot\hat{n})\;,$
$\displaystyle\Gamma_{ab,cd,\lambda}^{Q}(\kappa,\hat{s},\hat{v},\hat{w})=\frac{1}{8\pi}\int
d^{2}{\hat{n}}\,{\rm
e}^{i\kappa{\hat{n}}\cdot\hat{s}}\,e_{ab}^{*\lambda}\left({\hat{n}}\right)e_{cd}^{\lambda}\left({\hat{n}}\right)\,(i\hat{v}\cdot\hat{n})\,(i\hat{w}\cdot\hat{n})\;.$
(C.2)

In [29] it has been proven that the first two expressions can be rewritten in
the form

$\displaystyle\Gamma_{abcd,\lambda}^{M}(\kappa,\hat{s})=\Gamma_{abcd,\lambda}(\kappa,\vec{s})_{\rvert
s=1}\;,$
$\displaystyle\Gamma_{abcd,\lambda}^{D}(\kappa,\hat{s},\hat{v})=\left[\frac{1}{\kappa}\hat{v}_{i}\frac{\partial}{\partial
s_{i}}\Gamma_{abcd,\lambda}(\kappa,\vec{s})\right]_{\rvert s=1}\;,$ (C.3)

where it was found that

$\displaystyle\Gamma_{abcd,\lambda}(\kappa,\vec{s})$
$\displaystyle=f_{A}(\kappa
s)\frac{\delta_{ac}\delta_{bd}+\delta_{ad}\delta_{bc}}{2}+\frac{f_{B}(\kappa
s)}{s^{2}}\frac{\delta_{ac}s_{b}s_{d}+\delta_{bd}s_{a}s_{c}+\delta_{ad}s_{b}s_{c}+\delta_{bc}s_{a}s_{d}}{4}+$
$\displaystyle+\frac{f_{C}(\kappa
s)}{s^{4}}s_{a}s_{b}s_{c}s_{d}+\lambda\frac{f_{D}(\kappa
s)}{s}\frac{\delta_{ac}s_{e}\epsilon_{bde}+\delta_{bd}s_{e}\epsilon_{ace}+\delta_{ad}s_{e}\epsilon_{bce}+\delta_{bc}s_{e}\epsilon_{ade}}{4}+$
$\displaystyle+\lambda\frac{f_{E}(\kappa
s)}{s^{3}}\frac{s_{a}s_{c}s_{e}\epsilon_{bde}+s_{b}s_{d}s_{e}\epsilon_{ace}+s_{a}s_{d}s_{e}\epsilon_{bce}+s_{b}s_{c}s_{e}\epsilon_{ade}}{4}\;,$
(C.4)

and

$\displaystyle
f_{A}(y)=\frac{j_{1}(y)}{2y}+\frac{1-y^{2}}{2y^{2}}j_{2}(y)\;\;,\;\;f_{B}(y)=\frac{j_{1}(y)}{y}+\frac{5-y^{2}}{y^{2}}j_{2}(y)\;\;,\;\;f_{C}(y)=-7\frac{j_{1}(y)}{4y}+\frac{35-y^{2}}{y^{2}}j_{2}(y)$
$\displaystyle
f_{D}(y)=\frac{j_{1}(y)}{2}-\frac{j_{2}(y)}{2y}\;\;,\;\;f_{E}(y)=-\frac{j_{1}(y)}{2}+\frac{5j_{2}(y)}{2y}\;,$
(C.5)

where $j_{1}(y)$ and $j_{2}(y)$ are the spherical Bessel Functions, of,
respectively, the first and the second type. On being contracted with the
detector geometry, the response function acquires the form 131313We denote the
two response functions as $\mathcal{M}_{ij}^{\lambda}$ and
$\mathcal{D}_{ij}^{\lambda}$ in the notation of [29]. We note that the former
quantity coincides with the monopole response function that we have introduced
in eq. (5.9) up to a constant factor,
$\gamma_{00,ij,\lambda}=\frac{5}{2\sqrt{\pi}}\,{\cal M}_{ij}^{\lambda}$, while
$\mathcal{D}_{ij}^{\lambda}$ is a linear combination of
$\gamma_{1m,ij,\lambda}$.

$\displaystyle\mathcal{M}_{ij}^{\lambda}(k,\hat{s})=\Gamma_{abcd,\lambda}^{M}(\kappa,\hat{s})d_{i}^{ab}d_{j}^{cd}\;\;\;,\;\;\;\mathcal{D}_{ij}^{\lambda}(k,\hat{s},\hat{v})=\Gamma_{abcd,\lambda}^{D}(\kappa,\hat{s},\hat{v})d_{i}^{ab}d_{j}^{cd}\;,$
(C.6)

where [29]

$\displaystyle\mathcal{M}_{ij}^{\lambda}(k)$
$\displaystyle=f_{A}(\kappa)\operatorname{tr}\left[d_{i}d_{j}\right]+f_{B}(\kappa)\left(d_{i}\hat{s}_{ij}\right)^{a}\left(d_{j}\hat{s}_{ij}\right)^{a}+f_{C}(\kappa)\left(d_{i}\hat{s}_{ij}\hat{s}_{ij}\right)\left(d_{j}\hat{s}_{ij}\hat{s}_{ij}\right)$
$\displaystyle+\lambda
f_{D}(\kappa)\left[d_{i}d_{j}\right]^{ab}\epsilon_{abc}\hat{s}_{ij}^{c}+\lambda
f_{E}(\kappa)\left(d_{i}\hat{s}_{ij}\right)^{a}\left(d_{j}\hat{s}_{ij}\right)^{b}\epsilon_{abc}\hat{s}_{ij}^{c}\;,$
(C.7)

and

$\displaystyle\mathcal{D}_{ij,\lambda}$
$\displaystyle(\kappa,\hat{s},\hat{v})=f_{A}^{\prime}(\kappa)(\hat{v}\hat{s})\left(d_{i}d_{j}\right)^{aa}$
$\displaystyle+\left[f_{B}^{\prime}(\kappa)-\frac{2f_{B}(\kappa)}{\kappa}\right](\hat{v}\hat{s})\left(d_{i}\hat{s}\right)^{a}\left(d_{j}\hat{s}\right)^{a}+\frac{f_{B}(\kappa)}{\kappa}\left[\left(d_{i}\hat{v}\right)^{a}\left(d_{j}\hat{s}\right)^{a}+\left(d_{i}\hat{s}\right)^{a}\left(d_{j}\hat{v}\right)^{a}\right]$
$\displaystyle+\left[f_{C}^{\prime}(\kappa)-\frac{4f_{C}(\kappa)}{\kappa}\right](\hat{v}\hat{s})\left(d_{i}\hat{s}\hat{s}\right)\left(d_{j}\hat{s}\hat{s}\right)+\frac{2f_{C}(\kappa)}{\kappa}\left[\left(d_{i}\hat{v}\hat{s}\right)\left(d_{j}\hat{s}\hat{s}\right)+\left(d_{i}\hat{s}\hat{s}\right)\left(d_{j}\hat{v}\hat{s}\right)\right]$
$\displaystyle+\lambda(\hat{v}\hat{s})\left[f_{D}^{\prime}(\kappa)-\frac{f_{D}(\kappa)}{\kappa}\right]\left(d_{i}d_{j}\times\hat{s}\right)+\lambda\frac{f_{D}(\kappa)}{\kappa}\left(d_{i}d_{j}\times\hat{v}\right)$
$\displaystyle+\lambda\left[f_{E}^{\prime}(\kappa)-\frac{3f_{E}(\kappa)}{\kappa}\right](\hat{v}\hat{s})\left[\left(d_{i}\hat{s}\right)\left(d_{j}\hat{s}\right)\times\hat{s}\right]$
$\displaystyle+\lambda\frac{f_{E}(\kappa)}{\kappa}\left\\{\left[\left(d_{i}\hat{v}\right)\left(d_{j}\hat{s}\right)\times\hat{s}\right]+\left[\left(d_{i}\hat{s}\right)\left(d_{j}\hat{v}\right)\times\hat{s}\right]+\left[\left(d_{i}\hat{s}\right)\left(d_{j}\hat{s}\right)\times\hat{v}\right]\right\\}\;.$
(C.8)

In these expressions, we defined

$\displaystyle(d_{i}d_{j})^{ab}$ $\displaystyle\equiv
d_{i}^{ac}d_{j}^{cb}\;\;,\;\;(d_{i}\hat{s})^{a}\equiv
d_{i}^{ab}\hat{s}^{b}\;\;,\;\;(d_{i}\hat{s}\hat{s})\equiv
d_{i}^{ab}\hat{s}^{a}\hat{s}^{b}\;,$
$\displaystyle\left(d_{i}d_{j}\times\hat{s}\right)$ $\displaystyle\equiv
d_{i}^{ae}d_{j}^{eb}\epsilon_{abc}\hat{s}_{c}\;\;,\;\;\left[\left(d_{i}\hat{s}\right)\left(d_{j}\hat{s}\right)\times\hat{s}\right]\equiv\left(d_{i}\hat{s}\right)^{a}\left(d_{j}\hat{s}\right)^{b}\epsilon_{abc}\hat{s}_{c}\;,$
$\displaystyle(\hat{v}\hat{s})$
$\displaystyle\equiv\hat{v}^{i}\hat{s}^{i}\;\;,\;\;(d_{i}\hat{v})^{a}\equiv
d_{i}^{ab}\hat{v}^{b}\;\;,\;\;(d_{i}\hat{v}\hat{s})\equiv
d_{i}^{ab}\hat{v}^{a}\hat{s}^{b}\;,$
$\displaystyle\left(d_{i}d_{j}\times\hat{v}\right)$ $\displaystyle\equiv
d_{i}^{ae}d_{j}^{eb}\epsilon_{abc}\hat{v}_{c}\;\;,\;\;\left[\left(d_{i}\hat{s}\right)\left(d_{j}\hat{s}\right)\times\hat{v}\right]\equiv\left(d_{i}\hat{s}\right)^{a}\left(d_{j}\hat{s}\right)^{b}\epsilon_{abc}\hat{v}_{c}\;,$
$\displaystyle\left[\left(d_{i}\hat{v}\right)\left(d_{j}\hat{s}\right)\times\hat{s}\right]$
$\displaystyle\equiv\left(d_{i}\hat{v}\right)^{a}\left(d_{j}\hat{s}\right)^{b}\epsilon_{abc}\hat{s}_{c}\;\;,\;\;\left[\left(d_{i}\hat{s}\right)\left(d_{j}\hat{v}\right)\times\hat{s}\right]\equiv\left(d_{i}\hat{s}\right)^{a}\left(d_{j}\hat{v}\right)^{b}\epsilon_{abc}\hat{s}_{c}\;.$
(C.9)

We extend the procedure of eq. (C.3) to the quadrupole contribution
141414Analytic response functions to higher multipoles can be obtained from
higher derivatives of $\Gamma_{abcd,\lambda}$ with respect to $\vec{s}$.

$\displaystyle\Gamma_{ab,cd,\lambda}^{Q}(\kappa,\hat{s},\hat{v},\hat{w})=\left[\frac{1}{\kappa^{2}}\hat{v}_{i}\hat{w}_{j}\frac{\partial^{2}}{\partial
s_{i}\partial s_{j}}\Gamma_{abcd,\lambda}(\kappa,\vec{s})\right]_{\rvert
s=1}\;,$ (C.10)

obtaining, after projecting on the geometry of the pair of detectors

$\displaystyle\mathcal{Q}_{ij,\lambda}(\kappa,\hat{s},\hat{v},\hat{w})$
$\displaystyle=\Gamma_{abcd,\lambda}^{Q}(\kappa,\hat{s},\hat{v},\hat{w})d_{i}^{ab}d_{j}^{cd}$
$\displaystyle=\left.\mathcal{Q}_{ij,\lambda}(\kappa,\hat{s},\hat{v},\hat{w})\right|_{\text{Parity
invariant
}}+\left.\mathcal{Q}_{ij,\lambda}(\kappa,\hat{s},\hat{v},\hat{w})\right|_{\text{Parity
violating}}\;,$ (C.11)

where the parity-invariant and parity-violating terms are given, respectively,
by

$\displaystyle\left.\mathcal{Q}_{ij,\lambda}(\kappa,\hat{s},\hat{v},\hat{w})\right|_{\text{Parity
invariant }}=$
$\displaystyle\left\\{\frac{f_{A}^{\prime}(\kappa)}{\kappa}(\hat{v}\hat{w})+\left[f_{A}^{\prime\prime}(\kappa)-\frac{f_{A}^{\prime}(\kappa)}{\kappa}\right](\hat{v}\hat{s})(\hat{w}\hat{s})\right\\}\left(d_{i}d_{j}\right)^{aa}$
$\displaystyle+$
$\displaystyle\left\\{\left[\frac{f_{B}^{\prime}(\kappa)}{\kappa}-\frac{2f_{B}(\kappa)}{\kappa^{2}}\right](\hat{v}\hat{w})+\left[f_{B}^{\prime\prime}(\kappa)-\frac{5f_{B}^{\prime}(\kappa)}{\kappa}+\frac{8f_{B}(\kappa)}{\kappa^{2}}\right](\hat{v}\hat{s})(\hat{w}\hat{s})\right\\}\left(d_{i}\hat{s}\right)^{a}\left(d_{j}\hat{s}\right)^{a}$
$\displaystyle+$
$\displaystyle\left[\frac{f_{B}^{\prime}(\kappa)}{\kappa}-\frac{2f_{B}(\kappa)}{\kappa^{2}}\right](\hat{v}\hat{s})\left[\left(d_{i}\hat{w}\right)^{a}\left(d_{j}\hat{s}\right)^{a}+\left(d_{i}\hat{s}\right)^{a}\left(d_{j}\hat{w}\right)^{a}\right]$
$\displaystyle+$
$\displaystyle\left[\frac{f_{B}^{\prime}(\kappa)}{\kappa}-\frac{2f_{B}(\kappa)}{\kappa^{2}}\right](\hat{w}\hat{s})\left[\left(d_{i}\hat{v}\right)^{a}\left(d_{j}\hat{s}\right)^{a}+\left(d_{i}\hat{s}\right)^{a}\left(d_{j}\hat{v}\right)^{a}\right]$
$\displaystyle+$
$\displaystyle\frac{f_{B}(\kappa)}{\kappa^{2}}\left[\left(d_{i}\hat{v}\right)^{a}\left(d_{j}\hat{w}\right)^{a}+\left(d_{i}\hat{w}\right)^{a}\left(d_{j}\hat{v}\right)^{a}\right]$
$\displaystyle+$
$\displaystyle\left\\{\left[\frac{f_{C}^{\prime}(\kappa)}{\kappa}-\frac{4f_{C}(\kappa)}{\kappa^{2}}\right](\hat{v}\hat{w})+\left[f_{C}^{\prime\prime}(\kappa)-\frac{9f_{C}^{\prime}(\kappa)}{\kappa}+\frac{24f_{C}(\kappa)}{\kappa^{2}}\right](\hat{v}\hat{s})(\hat{w}\hat{s})\right\\}\left(d_{i}\hat{s}\hat{s}\right)\left(d_{j}\hat{s}\hat{s}\right)$
$\displaystyle+2\left[\frac{f_{C}^{\prime}(\kappa)}{\kappa}-\frac{4f_{C}(\kappa)}{\kappa^{2}}\right](\hat{v}\hat{s})\left[\left(d_{i}\hat{w}\hat{s}\right)\left(d_{j}\hat{s}\hat{s}\right)+\left(d_{i}\hat{s}\hat{s}\right)\left(d_{j}\hat{w}\hat{s}\right)\right]$
$\displaystyle+2\left[\frac{f_{C}^{\prime}(\kappa)}{\kappa}-\frac{4f_{C}(\kappa)}{\kappa^{2}}\right](\hat{w}\hat{s})\left[\left(d_{i}\hat{v}\hat{s}\right)\left(d_{j}\hat{s}\hat{s}\right)+\left(d_{i}\hat{s}\hat{s}\right)\left(d_{j}\hat{v}\hat{s}\right)\right]$
$\displaystyle\quad+\frac{2f_{C}(\kappa)}{\kappa^{2}}\left[\left(d_{i}\hat{v}\hat{w}\right)\left(d_{j}\hat{s}\hat{s}\right)+\left(d_{i}\hat{s}\hat{s}\right)\left(d_{j}\hat{v}\hat{w}\right)+2\left(d_{i}\hat{v}\hat{s}\right)\left(d_{j}\hat{w}\hat{s}\right)+2\left(d_{i}\hat{w}\hat{s}\right)\left(d_{j}\hat{v}\hat{s}\right)\right]\;,$
(C.12)

and by

$\displaystyle\left.\mathcal{Q}_{ij,\lambda}(\kappa,\hat{s},\hat{v},\hat{w})\right|_{\text{Parity
violating }}=$
$\displaystyle\lambda\left\\{\left[\frac{f_{D}^{\prime}(\kappa)}{\kappa}-\frac{f_{D}(\kappa)}{\kappa^{2}}\right](\hat{v}\hat{w})+\left[f_{D}^{\prime\prime}(\kappa)-\frac{3f_{D}^{\prime}(\kappa)}{\kappa}+\frac{3f_{D}(\kappa)}{\kappa^{2}}\right](\hat{v}\hat{s})(\hat{w}\hat{s})\right\\}\left(d_{i}d_{j}\times\hat{s}\right)$
$\displaystyle+\lambda\left[\frac{f_{D}^{\prime}(\kappa)}{\kappa}-\frac{f_{D}(\kappa)}{\kappa^{2}}\right](\hat{v}\hat{s})\left(d_{i}d_{j}\times\hat{w}\right)+\lambda\left[\frac{f_{D}^{\prime}(\kappa)}{\kappa}-\frac{f_{D}(\kappa)}{\kappa^{2}}\right](\hat{w}\hat{s})\left(d_{i}d_{j}\times\hat{v}\right)$
$\displaystyle+\lambda\left\\{\left[\frac{f_{E}^{\prime}(\kappa)}{\kappa}-\frac{3f_{E}(\kappa)}{\kappa^{2}}\right](\hat{v}\hat{w})+\left[f_{E}^{\prime\prime}(\kappa)-\frac{7f_{E}^{\prime}(\kappa)}{\kappa}+\frac{15f_{E}(\kappa)}{\kappa^{2}}\right](\hat{v}\hat{s})(\hat{w}\hat{s})\right\\}\left[\left(d_{i}\hat{s}\right)\left(d_{j}\hat{s}\right)\times\hat{s}\right]$
$\displaystyle+\lambda\left[\frac{f_{E}^{\prime}(\kappa)}{\kappa}-\frac{3f_{E}(\kappa)}{\kappa^{2}}\right](\hat{v}\hat{s})\left\\{\left[\left(d_{i}\hat{w}\right)\left(d_{j}\hat{s}\right)\times\hat{s}\right]+\left[\left(d_{i}\hat{s}\right)\left(d_{j}\hat{w}\right)\times\hat{s}\right]+\left[\left(d_{i}\hat{s}\right)\left(d_{j}\hat{s}\right)\times\hat{w}\right]\right\\}$
$\displaystyle+\lambda\left[\frac{f_{E}^{\prime}(\kappa)}{\kappa}-\frac{3f_{E}(\kappa)}{\kappa^{2}}\right](\hat{w}\hat{s})\left\\{\left[\left(d_{i}\hat{v}\right)\left(d_{j}\hat{s}\right)\times\hat{s}\right]+\left[\left(d_{i}\hat{s}\right)\left(d_{j}\hat{v}\right)\times\hat{s}\right]+\left[\left(d_{i}\hat{s}\right)\left(d_{j}\hat{s}\right)\times\hat{v}\right]\right\\}$
$\displaystyle+\lambda\frac{f_{E}^{\prime}(\kappa)}{\kappa^{2}}\Bigg{\\{}\left[\left(d_{i}\hat{v}\right)\left(d_{j}\hat{s}\right)\times\hat{w}\right]+\left[\left(d_{i}\hat{s}\right)\left(d_{j}\hat{v}\right)\times\hat{w}\right]+\left[\left(d_{i}\hat{w}\right)\left(d_{j}\hat{s}\right)\times\hat{v}\right]+\left[\left(d_{i}\hat{s}\right)\left(d_{j}\hat{w}\right)\times\hat{v}\right]$
$\displaystyle+\left[\left(d_{i}\hat{v}\right)\left(d_{j}\hat{w}\right)\times\hat{s}\right]+\left[\left(d_{i}\hat{w}\right)\left(d_{j}\hat{v}\right)\times\hat{s}\right]\Bigg{\\}}\;.$
(C.13)

Beside the quantities defined in eq. (C.9), we have also introduced

$\displaystyle(\hat{v}\hat{w})$
$\displaystyle\equiv\hat{v}^{i}\hat{w}^{i}\;\;,\;\;(d_{i}\hat{v}\hat{w})\equiv
d_{i}^{ab}\hat{v}^{a}\hat{w}^{b}\;\;,\;\;\left[\left(d_{i}\hat{v}\right)\left(d_{j}\hat{s}\right)\times\hat{w}\right]\equiv\left(d_{i}\hat{v}\right)^{a}\left(d_{j}\hat{s}\right)^{b}\epsilon_{abc}\hat{w}_{c}\;,$
$\displaystyle\left[\left(d_{i}\hat{s}\right)\left(d_{j}\hat{v}\right)\times\hat{w}\right]$
$\displaystyle\equiv\left(d_{i}\hat{s}\right)^{a}\left(d_{j}\hat{v}\right)^{b}\epsilon_{abc}\hat{w}_{c}\;\;,\;\;\left[\left(d_{i}\hat{v}\right)\left(d_{j}\hat{w}\right)\times\hat{s}\right]\equiv\left(d_{i}\hat{v}\right)^{a}\left(d_{j}\hat{w}\right)^{b}\epsilon_{abc}\hat{s}_{c}\;.$
(C.14)

For a network that also contains triangular-shaped interferometers, similar to
eqs. (C.6) and (C.11) we introduce

$\displaystyle\mathcal{M}_{ij,O}^{\lambda}(k,\hat{s})=\Gamma_{abcd,\lambda}^{M}(\kappa,\hat{s})d_{i,O}^{ab}d_{j}^{cd}\;\;,\;\;\mathcal{M}_{ij,OO^{\prime}}^{\lambda}(k,\hat{s})=\Gamma_{abcd,\lambda}^{M}(\kappa,\hat{s})d_{i,O}^{ab}d_{j,O^{\prime}}^{cd}\;,$
(C.15)

and analogously for the dipole and quadrupole terms.

###  C.1 Monopole response function for a chiral SGWB

We conclude this appendix with the explicit expression for the response
function to a chiral isotropic SGWB, which we separate in a parity even and a
parity odd contribution

$\displaystyle\gamma_{ij,00,\lambda}(f)$
$\displaystyle\equiv\gamma_{ij,I}(f)+\lambda\gamma_{ij,V}(f)\;,$ (C.16)

From the first of eq. (C.1), the first of eq. (C.6), and eq. (C), the two
terms acquire the form

$\displaystyle\gamma_{ij,I}(f)$
$\displaystyle\equiv\frac{5}{2\sqrt{\pi}}\Bigg{[}f_{A}(\kappa)\operatorname{tr}\left[d_{i}d_{j}\right]+f_{B}(\kappa)\left(d_{i}\hat{s}_{ij}\right)^{a}\left(d_{j}\hat{s}_{ij}\right)^{a}+f_{C}(\kappa)\left(d_{i}\hat{s}_{ij}\hat{s}_{ij}\right)\left(d_{j}\hat{s}_{ij}\hat{s}_{ij}\right)\Bigg{]}\;,$
$\displaystyle\gamma_{ij,I}(f)$
$\displaystyle\equiv\frac{5}{2\sqrt{\pi}}\Bigg{[}f_{D}(\kappa)\left[d_{i}d_{j}\right]^{ab}\epsilon_{abc}\hat{s}_{ij}^{c}+f_{E}(\kappa)\left(d_{i}\hat{s}_{ij}\right)^{a}\left(d_{j}\hat{s}_{ij}\right)^{b}\epsilon_{abc}\hat{s}_{ij}^{c}\Bigg{]}\;,$
(C.17)

where we recall that $\hat{s}_{ij}$ is the unit vector in the direction of the
baseline between the two interferometers, see eq. (5.12).

##  Appendix D Geometry of the detectors

In this appendix, we specify the position and orientation of the detectors
considered in this work. Let $\theta_{i}$, and $\phi_{i}$ be the latitude and
longitude of the $i$-th detector, respectively. We define latitude so that the
north pole is at $\theta=0$, while the equator is at $\theta=\pi/2$, while we
use the standard notation for longitude, with $\phi=0$ corresponding to
Greenwich meridian. We then consider a Cartesian system centred at the centre
of the Earth (assumed to be a sphere), with the $x-$ axis (resp., the
$y-$axis) points toward the location of $0$ longitude (resp., $\pi/2$
longitude) on the equator, and with the $z-$axis pointing toward the north
pole. In these Cartesian coordinates, the north and east directions for a
location specified by latitude and longitude $\left\\{\theta,\,\phi\right\\}$
on the surface of the Earth are given by

$\displaystyle{\hat{v}}_{{\rm
north}}=\left\\{-\sin\theta\,\cos\phi,\,-\sin\theta\,\sin\phi,\,\cos\theta\right\\}\;\;,\;\;{\hat{v}}_{{\rm
east}}=\left\\{-\sin\phi,\,\cos\phi,\,0\right\\}\;.$ (D.1)

We can then specify the orientation of the detector arm ${\hat{u}}_{X}$
through an angle $\beta$ that it forms with the north direction,

${\hat{u}}_{X}=\cos\left(\beta\right)\,{\hat{v}}_{\rm
north}+\sin\left(\beta\right)\,{\hat{v}}_{\rm east}\;,$ (D.2)

and for what concerns the L-shaped detectors, the other arm direction will be

${\hat{u}}_{Y}=-\sin\left(\beta\right)\,{\hat{v}}_{\rm
north}+\cos\left(\beta\right)\,{\hat{v}}_{\rm east}\;,$ (D.3)

while for triangular-shaped interferometers, we follow the conventions
indicated in figure 7 to define the other two arm directions ${\hat{u}}_{Y}$
and ${\hat{u}}_{Z}$.

Figure 7: Our convention for the orientation of the three ET arms, with the
labels “N” and “E” indicating, respectively, the north and east direction at
the location of the interferometer, and with $\beta$ being the angle formed by
the direction ${\hat{u}}_{X}$ of the “first” arm and the north direction.

Using this convention, in table 8, we display the actual locations of the
existing interferometers LIGO, Virgo and KAGRA, with the ET in the proposed
site in Sardinia, Italy ([48]). In the table, we set $\beta=0$, with $\beta$
defined in eq. (D.2).

| Central location | $\\{-0.338,-0.600,0.725\\}$
---|---|---
LIGO Hanford (or CE) | First Arm | $\\{-0.224,0.799,0.557\\}$
| Second Arm | $\\{-0.914,0.0261,-0.405\\}$
| Central location | $\\{-0.0116,-0.861,0.508\\}$
LIGO Livingston | First Arm | $\\{-0.953,-0.144,-0.266\\}$
| Second Arm | $\\{0.302,-0.488,-0.819\\}$
| Central location | $\\{0.712,0.132,0.690\\}$
Virgo | First Arm | $\\{-0.701,0.201,0.684\\}$
| Second Arm | $\\{-0.0485,-0.971,0.236\\}$
| Central location | $\\{-0.591,0.546,0.594\\}$
KAGRA | First Arm | $\\{-0.390,-0.838,0.382\\}$
| Second Arm | $\\{0.706,-0.00580,0.709\\}$
| Central location | $\\{0.751,0.125,0.649\\}$
ET - Sardinia | First Arm | $\\{-0.640,-0.106,0.761\\}$
| Second Arm | $\\{0.178,0.908,-0.381\\}$
| Third Arm | $\\{0.462,-0.801,-0.381\\}$
Figure 8: Cartesian coordinates of the unit-vectors specifying the positions
of the interferometers and the direction of their arms in the coordinate
system described in this appendix. For ET, we consider the configuration with
$\beta=0$, with $\beta$ defined in eq. (D.2), located in the proposed site in
Sardinia, Italy ([48]).

## References

  * Abbott et al. [2016] B. P. Abbott, R. Abbott, T. D. Abbott, M. R. Abernathy, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso, R. X. Adhikari, et al. (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett. 116, 061102 (2016), URL https://link.aps.org/doi/10.1103/PhysRevLett.116.061102.
  * Abbott et al. [2021] R. Abbott, T. D. Abbott, S. Abraham, F. Acernese, K. Ackley, A. Adams, C. Adams, R. X. Adhikari, V. B. Adya, C. Affeldt, et al. (LIGO Scientific Collaboration, Virgo Collaboration, and KAGRA Collaboration), Phys. Rev. D 104, 022004 (2021), URL https://link.aps.org/doi/10.1103/PhysRevD.104.022004.