# Mapping the Gravitational-wave Sky with LISA: A Bayesian Spherical Harmonic
Approach

Sharan Banagiri1, Alexander Criswell1, Tommy Kuan1, Vuk Mandic1, Joseph D.
Romano2, Stephen R. Taylor3
1 School of Physics and Astronomy, University of Minnesota, Minneapolis,
Minnesota 55455, USA
2 Texas Tech University, Physics & Astronomy Department, Box 41051, Lubbock,
Texas 79409, USA
3 Department of Physics & Astronomy, Vanderbilt University, 2301 Vanderbilt
Place, Nashville, TN 37235, USA
E-mail<EMAIL_ADDRESS>
###### Abstract

The millihertz gravitational-wave frequency band is expected to contain a rich
symphony of signals with sources ranging from galactic white dwarf binaries to
extreme mass ratio inspirals. Many of these gravitational-wave signals will
not be individually resolvable. Instead, they will incoherently add to produce
stochastic gravitational-wave confusion noise whose frequency content will be
governed by the dynamics of the sources. The angular structure of the power of
the confusion noise will be modulated by the distribution of the sources
across the sky. Measurement of this structure can yield important information
about the distribution of sources on galactic and extra-galactic scales, their
astrophysics and their evolution over cosmic timescales. Moreover, since the
confusion noise is part of the noise budget of LISA, mapping it will also be
essential for studying resolvable signals. In this paper, we present a
Bayesian algorithm to probe the angular distribution of the stochastic
gravitational-wave confusion noise with LISA using a spherical harmonic basis.
We develop a technique based on Clebsch–Gordan coefficients to mathematically
constrain the spherical harmonics to yield a non-negative distribution, making
them optimal for expanding the gravitational-wave power and amenable to
Bayesian inference. We demonstrate these techniques using a series of
simulations and analyses, including recovery of simulated distributed and
localized sources of gravitational-wave power. We also apply this method to
map the gravitational-wave foreground from galactic white-dwarfs using a
simplified model of the galactic white dwarf distribution.

###### keywords:

gravitational waves, methods:statistical, methods: data analysis

††pubyear: 2021††pagerange: Mapping the Gravitational-wave Sky with LISA: A
Bayesian Spherical Harmonic Approach–References

##  1 Introduction

The upcoming space-based Laser Interferometer Space Antenna (LISA) (Amaro-
Seoane. et al., 2017) promises access to the millihertz gravitational-wave
(GW) frequency band, which is inaccessible to terrestrial ground-based
detectors like the Laser Interferometer Gravitational-wave Observatory (LIGO).
A rich collection of galactic and extra-galactic sources emit GWs at these
frequencies. Among them are double white dwarfs (DWDs) from both within the
Milky Way (Marsh, 2011; Korol et al., 2017) and neighboring satellite galaxies
(Korol et al., 2018, 2020), extreme mass ratio inspirals (Amaro-Seoane et al.,
2007; Gair et al., 2017; Babak et al., 2017), supermassive blackhole binaries
(Barausse et al., 2015; Klein et al., 2016), extragalactic stellar-mass binary
black hole (BBH) and binary neutron star inspirals (Lau et al., 2020), and
even exoplanets orbiting white dwarfs (Tamanini & Danielski, 2019; Danielski
et al., 2019). Cosmological sources like cosmic strings (Auclair et al.,
2020), cosmological phase transitions (Caprini et al., 2016) and primordial GW
backgrounds (Bartolo et al., 2016) from the early universe are also
potentially accessible by LISA, not to mention the prospect of multi-
wavelength GW science in conjunction with next-generation ground-based
detectors (see for example Fitz Axen et al. (2018))

The GWs from many astrophysical sources will not be individually detectable
and will overlap to form a confusion noise in the detector, usually called the
stochastic gravitational-wave background (Regimbau, 2011; Romano & Cornish,
2017). The stochastic GW background from binary inspirals is a major
scientific target for both current ground-based detectors like Advanced LIGO
(Aasi et al., 2015) and Advanced Virgo (Acernese et al., 2015), and for pulsar
timing arrays (Brazier et al., 2019; Manchester et al., 2013; Kramer &
Champion, 2013). The current limits from ground-based detectors are
$\Omega_{\text{GW}}(f=25\text{ Hz})<4.8\times 10^{-8}$ (Abbott et al., 2019a)
for an isotropic background and $h_{0}<(3.6\textendash 4.7)\times 10^{-25}$
for point sources (Abbott et al., 2019b). The latest results from pulsar
timing arrays are from the NANOGrav collaboration which see strong evidence
for a common-spectrum red-noise process with a median strain of $1.92\times
10^{-15}$ at a frequency of $1\text{ yr}^{-1}$ (Arzoumanian et al., 2020).

In the case of LISA, there is also a GW foreground from galactic DWDs that
stands above the instrumental noise for a part of the LISA band while still
being stochastic in nature. While this foreground is considered an
inconvenient noise source for resolvable signals, it is of astrophysical
interest in its own right and contains useful information about the physical
and spectral distribution of the DWDs (Benacquista & Holley-Bockelmann, 2006;
Breivik et al., 2020b). For the rest of the paper, we will use the term
stochastic gravitational-wave confusion noise (SGCN) to refer to both
stochastic backgrounds and foregrounds collectively.

The angular structure of the GW power from a SGCN will directly follow the
distribution of the sources which generates it. The antenna patterns of the
detectors are not isotropic and also change throughout LISA’s orbit, which
ensures that the angular structure can, in principle, be measured with enough
integration time. Spherical harmonic functions are a natural basis to describe
the distribution of power on the sky and have been frequently used in
algorithms developed to measure these anisotropies, both for LIGO and LISA,
albeit in a frequentist manner (Thrane et al., 2009; Ungarelli & Vecchio,
2001; Kudoh & Taruya, 2005; Taruya & Kudoh, 2005; Taruya, 2006; Renzini &
Contaldi, 2018). In particular, recently Contaldi et al. (2020) developed a
frequentist maximum likelihood method to map GW power with LISA.

In this paper, we present a Bayesian algorithm to map the power of an SGCN
using a spherical harmonic basis. There are several advantages to developing a
Bayesian version of this method, especially in the case of LISA, where the
galactic foreground dominates. First, the Bayesian version can be better
integrated with global analyses designed to extract multiple resolvable
signals (Littenberg et al., 2020), in order to map the foreground
simultaneously along with them. Accounting for the foreground in this way can
be crucial for accurately inferring the properties of the resolvable signals.
Frequentist searches also generally require the inversion of a Fisher matrix
connecting different sky directions or harmonics. The poor angular sensitivity
of GW detectors, LISA included, creates degeneracies that make this inversion
mathematically ill-conditioned, necessitating the use of techniques like
singular value decomposition (see for example Thrane et al. (2009); Contaldi
et al. (2020)). These degeneracies are better accommodated with a Bayesian
approach, which requires no such inversion. Additionally, the angular
sensitivity of the spherical harmonic expansion is set by cutting off the
expansion in $\ell$ at some $\ell_{\text{max}}$ parameter value. This
parameter value is chosen in a somewhat ad hoc way in frequentist searches but
can be much more naturally accommodated in Bayesian searches by allowing the
data to determine it.

Historically, one hindrance of a Bayesian spherical harmonic implementation
has been that the generic expansion describes a complex field on the sky,
while GW power is non-negative by definition. We demonstrate a way to
mathematically impose this constraint using Clebsch–Gordan coefficients in a
Bayesian spherical harmonic analysis. This mathematical technique was used
recently to measure the angular distribution of GW detections by LIGO-Virgo
(Payne et al., 2020), while a similar method was also recently used in the
pulsar timing array band (Taylor et al., 2020). We also introduce the Bayesian
LISA Pipeline (BLIP) designed to simulate LISA data, perform the spherical
harmonic analysis on the simulated data, and conduct Bayesian inference to
recover the simulated parameters.

The rest of the paper is structured as follows. In Sec.2 we review SGCNs and
the spherical harmonic basis and calculate the detector response function of
LISA to GW power in the spherical harmonic basis. In Sec. 3 we calculate the
Clebsch–Gordan decomposition for non-negative fields and develop the
parameterization necessary for Bayesian inference. Section 4 introduces the
BLIP pipeline and discusses the likelihood function and the configuration used
for the analyses in this paper. Section 5 demonstrates measurement of
anisotropies in simulated LISA data in the spherical harmonic basis with the
Clebsch–Gordan decomposition. Section 6 discusses applications of the
technique to the galactic foreground followed by a discussion and conclusion
in Sec. 7.

##  2 Stochastic gravitational-wave confusion noise

Astrophysical SGCNs result from an incoherent superposition of GWs from many
disparate sources which are not individually resolvable (Regimbau, 2011)
111which implies that the strength of the background depends on the
sensitivity of the detector.. Appealing to the central limit theorem, one can
characterize the SGCN as colored Gaussian noise in the detectors (Allen &
Romano, 1999). This is usually a good assumption for sources in the LISA band
that overlap with one other. The metric perturbation at $(t,\textbf{x})$
corresponding to an SGCN can be written in the Fourier basis as,

$\begin{split}h_{ij}(t,\textbf{x})=\sum_{A}\int^{\infty}_{-\infty}df&\int
d^{2}n\,\tilde{h}_{A}(f,\textbf{n})\,e^{A}_{ij}(\textbf{n})\\\
&\times\,\exp\\{-2\pi if(t-\textbf{n}\cdot\textbf{x}/c)\\},\end{split}$ (2.1)

where $A=\\{+,\times\\}$ denotes polarization, $c$ is the speed of light, n is
the directional unit vector and $e^{A}_{ij}(\textbf{n})$ are the polarization
tensors 222A point on notation is appropriate here. Generally, most paper use
subscripts of $ab$ rather than $ij$ for the perturbation and polarization
tensors. Conventionally $a,b\in\\{1,2\\}$ to signify the transverse nature of
the GWs, in a coordinate system where n is aligned with the $z$ axis. This is
not a suitable choice for an SGCN which is an incoherent superposition of GWs
from all directions. So we use the subscript $i,j\in\\{1,2,3\\}$ to clarify
this point. See Sec. 7.8 of Maggiore (2008) for a more detailed description..
The frequency of the SGCN is represented by $f$ and
$\tilde{h}_{A}(f,\textbf{n})$ are the Fourier components of the perturbations
which satisfy the following,

$\begin{split}\langle\tilde{h}_{A}(f,\textbf{n})\rangle=&\,0,\\\
\langle\tilde{h}_{A}(f,\textbf{n})\,\tilde{h}^{*}_{A^{\prime}}(f^{\prime},\textbf{n}^{\prime})\rangle=&\frac{1}{2}\,S_{A}(f,\textbf{n})\,\delta_{A,A^{\prime}}\,\delta(f-f^{\prime})\,\delta^{2}(\textbf{n},\textbf{n}^{\prime}).\end{split}$
(2.2)

Under the assumption of Gaussianity the power spectrum $S_{A}(f,\textbf{n})$
is the main measurable quantity of an SGCN. For the rest of this paper we will
also assume that the SGCN is unpolarized i.e,
$S_{+}(f,\textbf{n})=S_{\times}(f,\textbf{n})=1/2\,S_{\text{GW}}(f,\textbf{n})$.
The power spectrum is conventionally characterized by the dimensionless energy
density $\Omega_{\text{GW}}(f,\textbf{n})$ per logarithmic frequency bin
(Allen & Romano, 1999), related to $S_{\text{GW}}(f,\textbf{n})$ via,

$\Omega_{\text{GW}}(f,\textbf{n})=\frac{2\pi^{2}f^{3}}{3H^{2}_{0}}S_{\text{GW}}(f,\textbf{n}),$
(2.3)

where $H_{0}$ is the Hubble constant. In general, the distribution of the
power on the sky will not be isotropic but will rather trace the distribution
of its sources. To describe the angular structure of the SGCN we will first
assume that its frequency and directional dependence can be factorized as

$\Omega_{\text{GW}}(f,\textbf{n})=\Omega(f)\mathcal{P}(\textbf{n}).$ (2.4)

The spectral shape of the SGCN is given by $\Omega(f)$ while
$\mathcal{P}(\textbf{n})$ describes the angular distribution of the
background, normalized so that:

$\int_{\mathcal{S}^{2}}d^{2}n\,\mathcal{P}(\textbf{n})=1.$ (2.5)

A common way to parameterize the spectral shape $\Omega(f)$ is to write it as
a power law:

$\Omega(f)=\Omega_{0}\left(\frac{f}{f_{0}}\right)^{\alpha},$ (2.6)

where $\alpha$ is the spectral index of the power law, $f_{0}$ is some
reference frequency and $\Omega_{0}=\Omega(f=f_{0})$. In particular stochastic
backgrounds and foregrounds from compact binaries are expected to follow a
power law with $\alpha=2/3$ (Phinney, 2001). For the rest of this paper, we
will assume the power-law spectral shape to hold.

Spherical harmonics provide a general orthonormal basis to parameterize an
arbitrary continuous and differentiable function on the two sphere:

$\mathcal{P}(\textbf{n})=\frac{1}{\sqrt{4\pi
a_{0,0}}}\sum_{\ell,m}a_{\ell,m}Y_{\ell,m}(\textbf{n}),$ (2.7)

where $\\{Y_{\ell,m}\\}$ are the spherical harmonic functions and the
normalization factor of $\sqrt{4\pi a_{0,0}}$ ensures that the expansion
satisfies Eq. 2.5. The coefficients $\\{a_{\ell,m}\\}$ are in general complex
numbers and characterize the distribution of the field on the sky. The
harmonics of positive and negative $m$ are related by:

$Y_{\ell,-m}(\textbf{n})=(-1)^{m}Y_{\ell,m}^{*}(\textbf{n}),$ (2.8)

The spherical harmonic functions form an orthonormal basis on the two sphere,

$\int
d^{2}n\,Y_{\ell,m}(\textbf{n})Y^{*}_{\ell^{\prime},m^{\prime}}(\textbf{n})=\delta_{\ell,\ell^{\prime}}\delta_{m,m^{\prime}}$
(2.9)

###  2.1 Detector response in the spherical harmonic basis

Since we decompose GW power in the orthonormal spherical harmonic basis, it is
also useful to calculate the detector response to each spherical harmonic
mode. First, we define the antenna pattern function of channel $I$ of the
detector as (Cornish & Larson, 2001; Cornish, 2002),

$F^{A}_{I}(f,t,\textbf{n})=D_{I}^{\text{ij}}(f,t,\textbf{n}):e^{A}_{\text{ij}}(\textbf{n}),$
(2.10)

where $D_{I}^{\text{ij}}(f,t,\textbf{n})$ is the detector response tensor. For
Michelson channels with arm orientations given by unit vectors u and v, the
response tensor is given by,

$\begin{split}D_{\text{mich}}(f,t,\textbf{n})=\frac{1}{2}\bigg{[}&\left(\textbf{u}\,\otimes\,\textbf{u}\right)\mathcal{T}(f,\textbf{u}\cdot\textbf{n})\\\
&-\left(\textbf{v}\,\otimes\,\textbf{v}\right)\mathcal{T}(f,\textbf{v}\cdot\textbf{n})\bigg{]},\end{split}$
(2.11)

assuming that the satellite motion is negligible during the round-trip light-
travel time between the satellites. The response tensor is a function of time
owing to the temporal variation of u and v as the satellites move in their
orbits. Here $\mathcal{T}(f,\textbf{u}\cdot\textbf{n})$ is the timing transfer
function of interferometric detectors to GWs, which for an equal arm detector
is given by (Schilling, 1997; Cornish & Larson, 2001):

$\begin{split}\mathcal{T}(f,\textbf{u}\cdot\textbf{n})&=\frac{1}{2}\bigg{[}\text{sinc}\left(\frac{f}{f_{*}}(1-\textbf{n}\cdot\textbf{u})\right)\exp\left(-i\frac{f}{f_{*}}(3+\textbf{n}\cdot\textbf{u})\right)\\\
+&\,\text{sinc}\left(\frac{f}{f_{*}}(1+\textbf{n}\cdot\textbf{u})\right)\exp\left(-i\frac{f}{f_{*}}(1+\textbf{n}\cdot\textbf{u})\right)\bigg{]},\end{split}$
(2.12)

where $L=2.5\times 10^{9}$ m is the arm length of LISA and $f_{*}\equiv c/\pi
L$. Since the SGCN is measured as the exces GW power in a detector, to detect
it, we will need to correlate two channels $I$ and $J$ 333Note that $I$ and
$J$ can be the same channel of LISA. In the case of LIGO and Virgo, we usually
only correlate distinct spatially-seperated interferometers, so the
correlation response function is called the overlap reduction function. The
response function to power distributed as $\mathcal{P}(\textbf{n})$ on the sky
is then given by,

$R^{IJ}(f,t)=\int\frac{d^{2}n}{4\pi}\,\mathcal{P}(\textbf{n})\left(\sum_{A}F^{A}_{I}(f,\textbf{n}){F^{A}_{J}}^{*}(f,\textbf{n})\right).$
(2.13)

Expanding $\mathcal{P}(\textbf{n})$ in the spherical harmonic basis as in Eq.
2.7 we can define the response function to the $Y_{\ell,m}$ mode to be,

$\mathcal{R}_{\ell,m}^{IJ}(f,t)=\int\frac{d^{2}n}{4\pi}\,Y_{\ell,m}(\textbf{n})\left(\sum_{A}F^{A}_{I}(f,\textbf{n}){F^{A}_{J}}^{*}(f,\textbf{n})\right).$
(2.14)

Finally combining Eqs. 2.3, 2.4, 2.7 and 2.14, the SGCN signal in the
correlation between channels $I$ and $J$ is,

$S_{IJ}^{\text{GW}}(f,t)=\frac{3H^{2}_{0}}{2\pi^{2}f^{3}}\frac{\Omega(f)}{\sqrt{4\pi
a_{0,0}}}\sum_{\ell,m}a_{\ell,m}\mathcal{R}_{\ell,m}^{IJ}(f,t).$ (2.15)

##  3 Clebsch-Gordan decomposition

The general spherical harmonic expansion describes a complex field. We can
constrain it to be real everywhere on the sky with the condition:

$a_{\ell,-m}=(-1)^{m}a_{\ell,m}^{*}.$ (3.1)

However, the decomposition of GW power needs to be not only real but also non-
negative for any direction on the sky, i.e., $\Omega(f,\textbf{n})\geq 0$.
This is especially important for LISA, as the GW power distribution will be
highly anisotropic due to the foreground from galactic binaries. Implementing
this constraint is also necessary for Bayesian inference with the spherical
harmonic basis. This is because the posterior should be zero for any set of
$a_{\ell,m}$’s that contains even the tiniest spot on the sky with negative GW
power. Previous work in the PTA band attempted to solve this problem by
numerically checking the sign of the GW power on the sky using a grid (Taylor
& Gair, 2013; Taylor et al., 2015) and assigning a probability of zero for a
given set of $a_{\ell,m}$’s if any of the pixels have negative power. However,
one can always use finer and finer grids to check this which makes this
solution computationally untenable. A solution that mathematically guarantees
non-negative power is much preferable. In this section, we describe a solution
to this problem using Clebsch–Gordan coefficients.

First we define a function $\mathcal{S}(\textbf{n})$ which is the square root
of the spherical harmonic expansion, i.e,

$\mathcal{S}(\textbf{n})=\left[\sum_{\ell,m}a_{\ell,m}Y_{\ell,m}(\textbf{n})\right]^{1/2}.$
(3.2)

We then expand $\mathcal{S}(\textbf{n})$ through its own spherical harmonic
expansion as,

$\mathcal{S}(\textbf{n})=\sum_{\ell,m}b_{\ell,m}Y_{\ell,m}(\textbf{n}).$ (3.3)

The necessary and sufficient condition for the GW power to be non-negative
then becomes $\mathcal{S}(\textbf{n})\in\mathbb{R}$ which implies that
$b_{\ell,-m}=(-1)^{m}b_{\ell,m}^{*}$. From Eqs. 3.2 and 3.3 we get,

$\sum_{L,M}a_{L,M}Y_{L,M}=\left(\sum_{\ell,m}b_{\ell,m}Y_{\ell,m}(\textbf{n})\right)^{2}.$
(3.4)

Expanding the right-hand side gives us:

$\sum_{L,M}a_{L,M}Y_{L,M}=\sum_{\ell,m}\sum_{\ell^{\prime},m^{\prime}}b_{\ell,m}b_{\ell^{\prime},m^{\prime}}Y_{\ell,m}(\textbf{n})Y_{\ell^{\prime},m^{\prime}}(\textbf{n}).$
(3.5)

The Clebsch-Gordan coefficients $C^{LM}_{\ell m,\ell^{\prime}m^{\prime}}$
provide the means to write products of spherical harmonics as a sum over
spherical harmonics 444Alternatively one can use Wigner-3j symbols, i.e.;

$\begin{split}Y_{\ell,m}(\textbf{n})Y_{\ell^{\prime},m^{\prime}}(\textbf{n})=\sum_{L=L_{\text{min}}}^{L_{\text{max}}}&\sqrt{\frac{(2\ell+1)(2\ell^{\prime}+1)}{4\pi(2L+1)}}\\\
&\times C^{LM}_{\ell
m,\ell^{\prime}m^{\prime}}\,C^{L0}_{\ell\,0,\ell^{\prime}\,0}Y_{L,M}(\textbf{n}).\end{split}$
(3.6)

The expansion obeys selection rules related to the symmetries of the rotation
group SO(3). We point to Ch. 16 of Arfken (2012) for a detailed discussion of
the mathematics of spherical harmonics and Clebsch-Gordan coefficients. The
selection rules can be listed as:

  * •

$M=m+m^{\prime}$

  * •

$L_{\text{min}}=\text{min}(|\ell-\ell^{\prime}|,|m+m^{\prime}|)$ and
$L_{\text{max}}=\ell+\ell^{\prime}$

  * •

$L$ is an integer

For compactness, let us define $\beta^{L,M}_{\ell m,\ell^{\prime}m^{\prime}}$
such that:

$\beta^{L,M}_{\ell
m,\ell^{\prime}m^{\prime}}=\sqrt{\frac{(2\ell+1)(2\ell^{\prime}+1)}{4\pi(2L+1)}}C^{LM}_{\ell
m,\ell^{\prime}m^{\prime}}\,C^{L0}_{\ell\,0,\ell^{\prime}\,0},$ (3.7)

when the selection rules are satisfied, but $\beta^{L,M}_{\ell
m,\ell^{\prime}m^{\prime}}=0$ otherwise. We can then write Eq. 3.6 as,

$Y_{\ell,m}(\textbf{n})Y_{\ell^{\prime},m^{\prime}}(\textbf{n})=\sum_{L,M}\beta^{L,M}_{\ell
m,\ell^{\prime}m^{\prime}}Y_{L,M}(\textbf{n}).$ (3.8)

Combining this with equation Eq. 3.5 gives:

$\begin{split}\sum_{L,M}a_{L,M}Y_{L,M}(\textbf{n})=\sum_{L,M}&\left(\sum_{\ell
m}\sum_{\ell^{\prime}m^{\prime}}b_{\ell,m}b_{\ell^{\prime},m^{\prime}}\beta_{L,M}^{\ell
m,\ell^{\prime}m^{\prime}}\right)\\\ &\times Y_{L,M}(\textbf{n}).\end{split}$
(3.9)

Since the set of $Y_{LM}$ form an orthonormal basis, this provides the recipe
for converting between $a_{\ell,m}$ and $b_{\ell,m}$:

$a_{L,M}=\sum_{\ell,m}\sum_{\ell^{\prime},m^{\prime}}b_{\ell,m}b_{\ell^{\prime},m^{\prime}}\beta_{L,M}^{\ell
m,\ell^{\prime}m^{\prime}}.$ (3.10)

Often, we want to impose an artificial cutoff on angular sensitivity of some
$\ell^{a}_{\text{max}}$ on the expansion in $a_{\ell,m}$’s. This cutoff can
correspond to estimated resolution limits of the detector itself, or it might
be astrophysically motivated. The corresponding cutoff of the expansion on the
$b_{\ell,m}$’s is taken to be $\ell^{b}_{\text{max}}=\ell^{a}_{\text{max}}/2$
assuming that $\ell^{a}_{\text{max}}$ is an even number. This is a consequence
of the second selection rule that $L_{\text{max}}=\ell+\ell^{\prime}$, which
implies that the cutoff on the expansion in the $b_{\ell,m}$’s should be half
that of the $a_{\ell,m}$’s if we want all terms higher than
$\ell^{a}_{\text{max}}$ in the latter to be zero. For the rest of this paper
we will use this relation and assume $\ell^{a}_{\text{max}}$ is even.

Since $\Omega(f,\textbf{n})$ is proportional to $\mathcal{S}^{2}(\textbf{n})$,
it is invariant under any transformation which leaves the latter invariant.
The constraint that $\mathcal{S}(\textbf{n})$ be real 555In fact we could have
chosen $\mathcal{S}$ such that $\Omega=|\mathcal{S}|^{2}$ in which case
$\mathcal{S}$ could be any general complex number. The degree of degeneracies
would be higher, so it would be a less economical choice than what we actually
make. makes $\Omega(f,\textbf{n})$ invariant under a parity transformation
$\\{b_{\ell,m}\\}\rightarrow\\{-b_{\ell,m}\\}$. However, this leftover
symmetry introduces degeneracies that induce multiple modes in the posterior
distribution. Moreover, since the $\\{a_{\ell,m}\\}$ expansion is normalized
as in Eq. 2.5, there is also a scale invariance within the $\\{b_{\ell,m}\\}$
space, i.e with the transformation
$\\{b_{\ell,m}\\}\rightarrow\\{\kappa\,b_{\ell,m}\\}$, where $\kappa$ is some
constant. One can break both these symmetries by fixing the value of one of
the $b_{\ell,m}$ coefficients. In this paper we choose to fix $b_{0,0}=1$.

A similar method to what is used here can also be applied to expand a
probability distribution or any other non-negative function on the two sphere.
Most of the results in this section will generalize to that case with the
additional requirement that the normalization of the distributions be one.
Indeed, the Clebsch-Gordan based spherical harmonic decomposition was recently
used to constrain anisotropies in the distribution of BBH progenitors using
events detected by LIGO-Virgo up-to-the second observing run in Payne et al.
(2020), and also to develop an optimized anisotropic pipeline for PTAs in
Taylor et al. (2020).

Figure 1: Skymaps for the simulation and analysis described in Sec. 5 with the
analysis time scale of one year. The map on the left is the simulated skymap
while the right shows the posterior median recovery skymap. Both maps show the
distribution of $\Omega(f=1$ mHz) in the solar system barycentric frame. The
full posteriors corresponding to these maps are shown in Fig. 3 Figure 2: The
variance of $a_{\ell,m}$’s sampled from the prior distribution is shown here,
without the multiplicative factor of $\Omega(f)$. The nearly isotropic
distribution of the variance across the sky demonstrates that the priors have
broad support for multiple spherical harmonic modes.

##  4 BLIP pipeline

This section briefly introduces the BLIP pipeline
666https://github.com/sharanbngr/blip, which is an independent Python-based
implementation for LISA data analysis and the details of it. The pipeline is
written to make it easy to add new GW signal models and likelihood models
along with simulating instrumental Gaussian noise in the time domain. The
instrumental noise is simulated through the acceleration and position noises
using the spectral form described in the LISA proposal (Amaro-Seoane. et al.,
2017). The functional forms for the power spectrum of the acceleration and
position noise are given by

$\begin{split}S_{p}(f)=&N_{p}\left[1+\left(\frac{2\,\text{mHz}}{f}\right)^{4}\right]\text{Hz}^{-1},\\\
S_{a}(f)=&\left[1+\left(\frac{0.4\,\text{mHz}}{f}\right)^{2}\right]\left[1+\left(\frac{f}{8\,\text{mHz}}\right)^{4}\right]\\\
&\times\frac{N_{a}}{(2\pi f)^{4}}\text{Hz}^{-1}\end{split}$ (4.1)

The noise levels can be set by the end-user by modifying $N_{p}$ and $N_{a}$,
but the current implementation assumes that they are the same in all satellite
links. In this paper we set $N_{p}=9\times 10^{-42}\text{ and }N_{a}=3.6\times
10^{-49}~{}\text{Hz}^{-4}$ to match the instrumental noise levels described in
the LISA proposal. The code also implements time-delay interferometry (TDI)
with Michelson, $X-Y-Z$ and $A-E-T$ (Tinto & Dhurandhar, 2005; Adams &
Cornish, 2010) channels and heliocentric rigid-body orbits of the LISA
satellites. This is implemented in an adiabatic manner by modeling the
satellites to be stationary for small segments of time ($\ll 1\text{ year}$)
both for signal simulation and recovery, and allowing them to move between the
time segments. The current implementation of orbits neglects the differential
time delay for laser light on the round trip between two satellites – i.e the
travel time difference from satellite A to B and from B to A – and also
neglects their breathing modes. The satellites’ orbital motion is especially
important in partially breaking the degeneracies of the antenna patterns when
detecting an anisotropic SGCN.

The strain data $d_{I}(t)$ from any channel $I$ is the sum of the instrumental
noise $n_{I}(t)$ and the GW signal in the channel $h^{GW}_{I}(t)$,

$d_{I}(t)=n_{I}(t)+h^{GW}_{I}(t).$ (4.2)

Due to the linear nature of the Fourier transform this relationship carries
over to the Fourier domain as well,

$\tilde{d}_{I}(f)=\tilde{n}_{I}(f)+\tilde{h}^{GW}_{I}(f).$ (4.3)

In the case of SGCNs, due to the Gaussian nature of both the GW strain and the
instrumental noise and under the assumption that they are uncorrelated, the
PSD of the data is the sum of the data and the noise PSDs,

$S_{II}(f)=S^{n}_{II}(f)+S^{GW}_{II}(f).$ (4.4)

A similar relations will also hold for cross correlation between channels $I$
and $J$ too,

$S_{IJ}(f)=S^{n}_{IJ}(f)+S^{GW}_{IJ}(f).$ (4.5)

Finally, the BLIP pipeline is built to facilitate Bayesian inference and
supports emcee (Foreman-Mackey et al., 2013) and dynesty (Speagle, 2020)
samplers. All results in this paper were made through the dynesty sampler. The
Clebsch–Gordan coefficients are implemented with the help of the Wigner module
of SymPy (Meurer et al., 2017)

Figure 3: Posteriors corresponding to skymaps in Fig. 1 with a signal
amplitude $\Omega_{0}=2\times 10^{-7}$ for the duration of $1$ year. This
corresponds to a single channel theoretical SNR $\approx 149$. The shaded
region in the one-dimensional posteriors are $95\%$ confidence levels while
the light and dark regions in the two-dimensional posteriors are one and two
sigma confidence levels respectively. The parameters $N_{p}$ and $N_{a}$ are
the posterior measurements for the position and acceleration noises
respectively using the functional forms described in Eq. 4.1. The parameters
$\alpha$ and $\Omega_{0}$ measure the spectral shape of the SGCN while the
rest of the parameters are measurements of $b_{\ell,m}$’s which describe the
distribution of GW power on the sky. The dashed green lines are the true
values of these parameters used when simulating the data.

###  4.1 Analysis configuration

For the remainder of this paper we will use $X-Y-Z$ TDI channels with a rigid
body orbiting configuration of LISA. To analyze the simulated data, we employ
Fourier transforms with a duration of $T_{\text{seg}}=1\times 10^{5}$ s with
the aforementioned adiabatic approximation within each segment. The sampling
frequency of the data is $f_{s}=0.25$ Hz. We also want to approximate the
covariance matrix of the data in this frequency-time analysis to be diagonal
across frequency and time, for which we require the auto-correlation time
scale to be much smaller than $T_{\text{seg}}$ (Banagiri et al., 2020). The
value of $T_{\text{seg}}=1\times 10^{5}$ s is thus chosen as a compromise
between the auto-correlation time scale of LISA and the motion of the
satellites. Each time segment is Hann-windowed before Fourier transformation.
We only consider the Fourier components to be between $f_{\text{min}}=2\times
10^{-4}$ Hz and $f_{\text{max}}=2\times 10^{-2}$ Hz.

###  4.2 Likelihood function

The power spectral density (PSD) and the cross-spectral density (CSD) of the
data are combinations of the GW power from the SGCN and of the instrumental
noise power (Cornish & Romano, 2013). Assuming that both are Gaussian, the
Fourier domain likelihood is based on the multi-dimensional complex Gaussian
distribution (Adams & Cornish, 2010):

$\begin{split}\mathcal{L}(\tilde{d}|N_{p},N_{a},\\{b_{\ell,m}\\})&=\prod_{t,f}\frac{1}{2\pi
T_{\text{seg}}|C(t,f)|}\\\
&\times\exp\left(-\frac{2\,\tilde{d}^{*}_{t,f}\,C(t,f)^{-1}\,\tilde{d}_{t,f}}{T_{\text{seg}}}\right).\end{split}$
(4.6)

Here
$\tilde{d}_{t,f}=[\tilde{d}_{X}(t,f),\tilde{d}_{Y}(t,f),\tilde{d}_{Z}(t,f)]$
is the array of data in the Fourier domain for the three channels measured in
the time segment labeled by $t$ and at frequency $f$. As previously mentioned,
the data is Fourier transformed in segments of duration $T_{\text{seg}}$, and
the product is across all frequency bins and time-segments. Here $C(t,f)$ is
the $3\times 3$ covariance matrix across the three channels, and as seen in
Eqs. 4.4 and 4.5, its elements are the sum of the signal spectral densities
(defined in Eq. 2.15) and the instrumental noise spectral densities
$S^{n}_{IJ}(f)$. We follow the derivation in Adams & Cornish (2010) of the
actual expressions for the noise spectral densities.

The term $|C(t,f)|$ is the determinant of the covariance matrix. In many
application of Bayesian inference to GW data, one usually can ignore the
overall normalization because it is constant and model independent. This is
not true in this case and correctly modeling the normalization is essential to
the the problem of Bayesian inference of the SGCN.

Figure 4: Posteriors for an analysis with a signal amplitude
$\Omega_{0}=2\times 10^{-7}$ for the duration of $2$ months. This corresponds
to a single channel theoretical SNR $\approx 59$. The shaded region in the
one-dimensional posteriors are $95\%$ confidence levels while the light and
dark regions in the two-dimensional posteriors are one and two sigma
confidence levels respectively. The bimodalities seen in the one-dimensional
posteriors are due to the parity symmetries described in Sec. 3 which are only
approximately broken in the limit of a weak signal or short integration time.
With a stronger signal, the breaking of the symmetry becomes more complete and
the degenerate modes go away as the sampler finds the right mode, as can be
seen in the posterior for the $1$ year run (Fig. 3).

Figure 5: Recovered skymaps for simulated localized sources using three months
of simulated data. The red star indicates the true position of the source. The
map on the left uses a cutoff of $\ell_{max}^{a}=4$ while the map on the right
uses $\ell_{max}^{a}=6$. Consequently, the power is smeared to a larger extent
on the sky for the former compared to the latter.

##  5 Simulation and detection

###  5.1 Validation

We first validate the Clebsch–Gordan technique by recovering an ad hoc
distribution of power, simulated with a power-law spectral-index SGCN with
$\ell_{\text{max}}^{a}=4$. The simulated $\\{b_{\ell,m}\\}$ coefficients are
777in healpix order ($1.0,\,0.75,\,0.5,\,0.7j,\,0.7-0.3j,\,1.1j)$ which yield
the GW power distribution on the sky shown in Fig. 1. These values are chosen
only to validate the ability of the algorithm to recover an arbitrary
distribution of power on the sky. The spectral index of the power law is
chosen to be consistent with binary inspiral at $\alpha=2/3$ with
$\Omega(f=25\text{Hz})=2\times 10^{-7}$.

The spherical harmonic coefficients $\\{b_{\ell,m}\\}$ are complex if $m\neq
0$ and thus have two degrees of freedom. We parameterize them by their
amplitude $|b_{\ell,m}|$ and phase $\phi_{\ell,m}$. We set uniform priors
between $[0,3]$ on the amplitude and uniform priors between $[-\pi,\pi]$ on
the phase. For the modes with $m=0$, i.e $b_{\ell,0}$’s which are real, we set
uniform priors between $[-3,3]$. The variance of the prior sky map due to
these choices of priors is shown in Fig. 2. This shows broad support for many
modes and implies that the priors are not peaked at any particular region in
the space.

We use the single channel signal-to-noise ratio (SNR) as a metric to
characterize the strength of the signal, defined as

$\text{SNR}=\sqrt{T_{\text{seg}}\sum_{t,f}\delta
f\,\frac{S^{GW}_{XX}}{S^{n}_{XX}}}$ (5.1)

where the summation is over time and frequency band used for analyzing the
data. Note that while the SNR is a metric of the strength of the SGCN, because
of the non-isotropic distribution of the SGCN, it does not grow as $\sim
T^{1/2}$ as would be expected in the isotropic case.

We ran two separate analyses with a duration of one year and two months.
Figures 3 and 4 show the posterior corner plots for the former and the latter
respectively. We witness bimodalities in the posteriors for the $b_{\ell,m}$
in the two-month run, which are related to the parity symmetry described in
Sec. 3. While our choice of fixing $b_{0,0}=1$ breaks this symmetry in
principle, in practice, this breaking is only approximate and can fail in the
limit of a weak signal or a small amount of observation time. This is because
while two modes may be related by a sign change due to parity, we are limited
by the ability of LISA to resolve this relative sign (with respect to
$b_{0,0}$). However, with enough time or a stronger signal, the breaking of
the symmetry becomes complete, and the sampler finds the right mode as seen in
the posterior for the one year run. The corresponding injected skymap, and the
recovered median posterior skymap (for 1 year long run) are shown in Fig. 1.
The median posterior map is the skymap generated from median values of
posterior samples.

###  5.2 Localized sources

Next, we test the ability to recover signals from localized sources of GW
power, often referred to as point sources in the literature. Examples of such
sources could be combined gravitational radiation from a large globular
cluster or a nearby galaxy. The GW power from these sources would be an
incoherent superposition of all the GWs emitted by individual sources within
them. Such sources are generally localized to much smaller angular scales than
is possible to resolve with LISA. Thus when we map them with spherical
harmonic methods, the power is usually smeared on larger angular scale as
determined by the strength of the signal, the integration time and the
$\ell_{\text{max}}^{a}$ scale where we cutoff the expansion. The dependence on
the latter is demonstrated in Fig. 5 where three months of simulated data
containing a point source signal is analyzed using
$\ell_{\text{max}}^{a}=4~{}\text{and}~{}6$, with the same $\Omega(f)$ in both
cases. As expected the smearing of the power is smaller in the case of
$\ell_{\text{max}}^{a}=6$. In both cases, we still assumed $\alpha=2/3$ power-
law spectrum.

##  6 Galactic Foreground

The galactic white dwarf foreground is expected to be one of the primary
sources of confusion noise seen by LISA. Due to its strength and strong
anisotropy, characterizing this foreground will be a primary application of
the method described here. Accurately mapping the foreground will also be
important for searches of resolvable sources as it is a major source of noise.
The confusion noise in the detectors due to the galactic foreground will
exhibit a time modulation tied to its anisotropy as LISA orbits the sun.

We simulate the galactic foreground by modeling it as stochastic colored-
Gaussian noise with a power-law spectrum given by Eq. 2.6. Following Breivik
et al. (2020a, b) and McMillan (2011), we model the spatial distribution in
the galaxy using a disk+bulge model. However, instead of simulating an
ensemble of individual DWD systems, we instead distribute the DWD density in
the galactocentric frame using this model. The disk is modeled as,

$\rho_{d}(r,z)\propto\exp\left(-r/r_{h}\right)\exp\left(-z/z_{h}\right),$
(6.1)

where the radial scale height is assumed to be $r_{h}=2.9$ kpc, and the
vertical scale height $z_{h}=0.3$ kpc following the thin disk model of Breivik
et al. (2020a). The galactic bulge is modeled to be azimuthally symmetric,

$\rho_{b}(r)\propto\frac{\exp\left(-(r/r_{\text{cut}})^{2}\,\right)}{(1+r^{\prime}/r_{0})^{\gamma}},$
(6.2)

where $\gamma=1.8$, $r_{0}=0.075$ kpc, $r_{\text{cut}}=2.1$ kpc,
$r^{\prime}=\sqrt{r^{2}+(z/q)^{2}}$ and $q=0.5$.

The distribution is then transformed to the Solar System barycenter
coordinates. The DWD density at each point in this grid is modulated by the
distance to it from Earth and the grid is projected on to the sky to produce
an unnormalized distribution on the sky. The skymap is then normalized to one
and multiplied by the power spectrum of the foreground, modeled as an
$\alpha=2/3$ power law, producing the simulated GW power skymap seen on the
top panel of Fig. 6. The GW power distribution is then convolved with the
time-dependent LISA detector response to calculate the GW power in the
detector from which the simulated time-domain data is generated.

Here, we have simulated a year of galactic foreground data in the manner
described above, which we then analyzed using the Clebsch–Gordan method with
$\ell_{max}^{a}=4$. In the analysis, we have used $f_{\text{max}}$ of
$10^{-3}$ Hz to simulate the drop in the foreground due to mass transfer DWDs
in a simplistic manner. The bottom panel of Fig. 6 shows the recovered skymap.
We see that the overall shape is qualitatively well recovered, especially the
galactic bulge. This demonstrates the Clebsch–Gordan method’s efficacy in
mapping the galactic DWD foreground. We leave a more sophistical analysis of
recovering the foreground’s actual shape, and the exact frequency cut off
caused by mass transfer to a future paper.

##  7 Discussion & Conclusion

Figure 6: Skymaps for the simulation and recovery of the galactic DWD
foregrounds described in Sec.6 with one year of data. The top panel shows the
simulated skymap while the bottom shows posterior median recovery skymap. Both
maps show the distribution of $\Omega(f=1\times 10^{-3}$ Hz) in the solar
system barycentric frame. The bright spots in the map corresponds to the
galactic central bulge.

We have developed a Bayesian mapping algorithm using the spherical harmonic
basis that can optimally recover any arbitrary distribution of GW power on the
sky using LISA data, while imposing the physical constraint that the GW power
is non-negative in all directions on the sky. This method was validated
through a series of end-to-end simulations of different types of SGCNs.

While this paper developed mapping tools and validates them, there are several
directions to push this forward and apply them for astrophysical use. The
galactic foreground simulated in this paper is ultimately simplistic in that
it assumes that it is Gaussian and has the spectral shape of a power law with
$\alpha=2/3$. It will be essential to relax these assumptions and validate
this method on a realistic simulated foreground formed from catalogs of DWDs
generated with population synthesis codes. This might also require relaxing
the assumptions of Gaussianity. Additionally, as Breivik et al. (2020b) shows,
mapping the foreground on the sky can help constrain galactic structure in a
manner complementary to resolvable DWDs while also probing the stellar
evolution history of the galactic white dwarfs. It would be interesting to
study what mapping the foreground could tell us about the astrophysics of
galactic foreground with the optimal analysis developed in this work and using
realistic simulations. One can also expand this to study the astrophysics of
the galactic DWD population in a more model-independent manner by constraining
the properties of the foreground in narrow frequency bands.

It will also be essential to study how the algorithm’s angular sensitivity,
characterized by the $\ell_{max}^{a}$ parameter, will scale with the strength
of the stochastic confusion noise and the duration of the analysis. An optimal
way to do this could be by making $\ell_{max}^{a}$ an independent parameter to
be determined by data. Since the dimensionality of the spherical harmonic
parameter space itself depends on $\ell_{max}^{a}$, such an analysis could use
Reversible-Jump Markov-chain Monte Carlo approaches such as those developed in
Cornish & Littenberg (2015).

Finally, an important application would be a joint analysis with both
resolvable signals and the galactic DWD foreground. One of the core strengths
of mapping the foreground in a Bayesian manner is the possibility of
simultaneous inference of resolvable signals along with the galactic
foreground. This will allow an unbiased estimation of the properties of the
resolvable signals while accounting for the temporal modulation of the noise
due to the galactic foreground.

##  8 Acknowledgments

We are grateful to Katelyn Breivik for helpful discussions about the galactic
white dwarf foreground, and to Ethan Payne for discussion about the
Clebsch–Gordan decomposition. We are also thankful to Arianna Renzini for
useful discussion about using the spherical harmonic basis with LISA. All
corner plots were made with ChainConsumer (Hinton, 2016). The authors
acknowledge the Minnesota Supercomputing Institute (MSI) at the University of
Minnesota for providing computational resources that contributed to the
research reported within this paper. This work is supported by the NASA grant
80NSSC19K0318. SB acknowledges support by the Doctoral Dissertation Fellowship
at the University of Minnesota. AC acknowledges support by the National
Science Foundation under grant No. 1922512. SRT acknowledge support from NSF
AST-2007993, and a Dean’s Faculty Fellowship of Vanderbilt University’s
College of Arts & Science.

## References

  * Aasi et al. (2015) Aasi J., et al., 2015, Class. Quant. Grav., 32, 074001
  * Abbott et al. (2019a) Abbott B., et al., 2019a, Phys. Rev. D, 100, 061101
  * Abbott et al. (2019b) Abbott B., et al., 2019b, Phys. Rev. D, 100, 062001
  * Acernese et al. (2015) Acernese F., et al., 2015, Class. Quant. Grav., 32, 024001
  * Adams & Cornish (2010) Adams M. R., Cornish N. J., 2010, Phys. Rev. D, 82, 022002
  * Allen & Romano (1999) Allen B., Romano J. D., 1999, Phys. Rev. D, 59, 102001
  * Amaro-Seoane et al. (2007) Amaro-Seoane P., Gair J. R., Freitag M., Coleman Miller M., Mandel I., Cutler C. J., Babak S., 2007, Class. Quant. Grav., 24, R113
  * Amaro-Seoane. et al. (2017) Amaro-Seoane. P., et al., 2017, arXiv e-prints
  * Arfken (2012) Arfken G., 2012, Mathematical Methods for Physicists. Elsevier, doi:https://doi.org/10.1016/C2009-0-30629-7
  * Arzoumanian et al. (2020) Arzoumanian Z., et al., 2020, Astrophys. J. Lett., 905, L34
  * Auclair et al. (2020) Auclair P., et al., 2020, JCAP, 04, 034
  * Babak et al. (2017) Babak S., et al., 2017, Phys. Rev. D, 95, 103012
  * Banagiri et al. (2020) Banagiri S., Coughlin M. W., Clark J., Lasky P. D., Bizouard M., Talbot C., Thrane E., Mandic V., 2020, Mon. Not. Roy. Astron. Soc., 492, 4945
  * Barausse et al. (2015) Barausse E., Bellovary J., Berti E., Holley-Bockelmann K., Farris B., Sathyaprakash B., Sesana A., 2015, J. Phys. Conf. Ser., 610, 012001
  * Bartolo et al. (2016) Bartolo N., et al., 2016, JCAP, 12, 026
  * Benacquista & Holley-Bockelmann (2006) Benacquista M., Holley-Bockelmann K., 2006, Astrophys. J., 645, 589
  * Brazier et al. (2019) Brazier A., et al., 2019
  * Breivik et al. (2020a) Breivik K., et al., 2020a, Astrophys. J., 898, 71
  * Breivik et al. (2020b) Breivik K., Mingarelli C. M. F., Larson S. L., 2020b, Astrophys. J., 901, 4
  * Caprini et al. (2016) Caprini C., et al., 2016, JCAP, 04, 001
  * Contaldi et al. (2020) Contaldi C. R., Pieroni M., Renzini A. I., Cusin G., Karnesis N., Peloso M., Ricciardone A., Tasinato G., 2020, Phys. Rev. D, 102, 043502
  * Cornish (2002) Cornish N. J., 2002, Phys. Rev. D, 65, 022004
  * Cornish & Larson (2001) Cornish N. J., Larson S. L., 2001, Class. Quant. Grav., 18, 3473
  * Cornish & Littenberg (2015) Cornish N. J., Littenberg T. B., 2015, Class. Quant. Grav., 32, 135012
  * Cornish & Romano (2013) Cornish N. J., Romano J. D., 2013, Phys. Rev. D, 87, 122003
  * Danielski et al. (2019) Danielski C., Korol V., Tamanini N., Rossi E., 2019, Astron. Astrophys., 632, A113
  * Fitz Axen et al. (2018) Fitz Axen M., Banagiri S., Matas A., Caprini C., Mandic V., 2018, Phys. Rev. D, 98, 103508
  * Foreman-Mackey et al. (2013) Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, Publications of the Astronomical Society of the Pacific, 125, 306
  * Gair et al. (2017) Gair J. R., Babak S., Sesana A., Amaro-Seoane P., Barausse E., Berry C. P., Berti E., Sopuerta C., 2017, J. Phys. Conf. Ser., 840, 012021
  * Hinton (2016) Hinton S. R., 2016, The Journal of Open Source Software, 1, 00045
  * Klein et al. (2016) Klein A., et al., 2016, Phys. Rev. D, 93, 024003
  * Korol et al. (2017) Korol V., Rossi E. M., Groot P. J., Nelemans G., Toonen S., Brown A. G., 2017, Mon. Not. Roy. Astron. Soc., 470, 1894
  * Korol et al. (2018) Korol V., Koop O., Rossi E. M., 2018, Astrophys. J. Lett., 866, L20
  * Korol et al. (2020) Korol V., et al., 2020, Astron. Astrophys., 638, A153
  * Kramer & Champion (2013) Kramer M., Champion D. J., 2013, Class. Quant. Grav., 30, 224009
  * Kudoh & Taruya (2005) Kudoh H., Taruya A., 2005, Phys. Rev. D, 71, 024025
  * Lau et al. (2020) Lau M. Y., Mandel I., Vigna-Gómez A., Neijssel C. J., Stevenson S., Sesana A., 2020, Mon. Not. Roy. Astron. Soc., 492, 3061
  * Littenberg et al. (2020) Littenberg T., Cornish N., Lackeos K., Robson T., 2020, Phys. Rev. D, 101, 123021
  * Maggiore (2008) Maggiore M., 2008, Gravitational Waves: Volume 1: Theory and Experiments
  * Manchester et al. (2013) Manchester R. N., et al., 2013, Publ. Astron. Soc. Austral., 30, 17
  * Marsh (2011) Marsh T. R., 2011, Classical and Quantum Gravity, 28, 094019
  * McMillan (2011) McMillan P. J., 2011, Monthly Notices of the Royal Astronomical Society, 414, 2446
  * Meurer et al. (2017) Meurer A., et al., 2017, PeerJ Computer Science, 3, e103
  * Payne et al. (2020) Payne E., Banagiri S., Lasky P., Thrane E., 2020
  * Phinney (2001) Phinney E., 2001
  * Regimbau (2011) Regimbau T., 2011, Res. Astron. Astrophys., 11, 369
  * Renzini & Contaldi (2018) Renzini A., Contaldi C., 2018, Mon. Not. Roy. Astron. Soc., 481, 4650
  * Romano & Cornish (2017) Romano J. D., Cornish N. J., 2017, Living Rev. Rel., 20, 2
  * Schilling (1997) Schilling R., 1997, Classical and Quantum Gravity, 14, 1513
  * Speagle (2020) Speagle J. S., 2020, Monthly Notices of the Royal Astronomical Society, 493, 3132
  * Tamanini & Danielski (2019) Tamanini N., Danielski C., 2019, Nature Astronomy, 3, 858–866
  * Taruya (2006) Taruya A., 2006, Phys. Rev. D, 74, 104022
  * Taruya & Kudoh (2005) Taruya A., Kudoh H., 2005, Phys. Rev. D, 72, 104015
  * Taylor & Gair (2013) Taylor S. R., Gair J. R., 2013, Phys. Rev. D, 88, 084001
  * Taylor et al. (2015) Taylor S. R., et al., 2015, Phys. Rev. Lett., 115, 041101
  * Taylor et al. (2020) Taylor S. R., van Haasteren R., Sesana A., 2020, Phys. Rev. D, 102, 084039
  * Thrane et al. (2009) Thrane E., Ballmer S., Romano J. D., Mitra S., Talukder D., Bose S., Mandic V., 2009, Phys. Rev. D, 80, 122002
  * Tinto & Dhurandhar (2005) Tinto M., Dhurandhar S. V., 2005, Living Rev. Rel., 8, 4
  * Ungarelli & Vecchio (2001) Ungarelli C., Vecchio A., 2001, Phys. Rev. D, 64, 121501