Thus given dual sliced potentials $(r_{\theta}(\cdot),s_{\theta}(\cdot))$
which solve $\text{{SOT}}(\bar{\alpha},\bar{\beta})$, one can then perform
Frank-Wolfe updates (77) on $(f_{\theta},g_{\theta})$ and thus
$(f_{avg},g_{avg})$.

###### Proof.

Our goal is to compute the first order variation of the USOT functional.
First, we leverage Theorem 5 such that USOT reads

$\displaystyle\text{{USOT}}(\alpha,\beta)=\sup_{f_{\theta}(\cdot)\oplus
g_{\theta}(\cdot)\leq\text{{C}}_{1}}$
$\displaystyle\int\varphi_{1}^{\circ}\Big{(}\int_{\mathbb{S}^{d-1}}f_{\theta}(\theta^{\star}(x))\mathrm{d}\hat{\bm{\sigma}}_{K}(\theta)\Big{)}\mathrm{d}\alpha(x)$
(85)
$\displaystyle+\int\varphi_{2}^{\circ}\Big{(}\int_{\mathbb{S}^{d-1}}g_{\theta}(\theta^{\star}(y))\mathrm{d}\hat{\bm{\sigma}}_{K}(\theta)\Big{)}\mathrm{d}\beta(y)$
(86) $\displaystyle=\sup_{f_{\theta}(\cdot)\oplus
g_{\theta}(\cdot)\leq\text{{C}}_{1}}$
$\displaystyle\int\varphi_{1}^{\circ}\Big{(}f_{avg}(x)\Big{)}\mathrm{d}\alpha(x)+\int\varphi_{2}^{\circ}\Big{(}g_{avg}(y)\Big{)}\mathrm{d}\beta(y),$
(87)

where

$\displaystyle
f_{avg}(x)=\int_{\mathbb{S}^{d-1}}f_{\theta}(\theta^{\star}(x))\mathrm{d}\hat{\bm{\sigma}}_{K}(\theta),$
$\displaystyle
g_{avg}(y)=\int_{\mathbb{S}^{d-1}}g_{\theta}(\theta^{\star}(y))\mathrm{d}\hat{\bm{\sigma}}_{K}(\theta).$

From this, we derive the translation-invariant formulation as follows.

$\displaystyle\text{{USOT}}(\alpha,\beta)=\sup_{f_{\theta}(\cdot)\oplus
g_{\theta}(\cdot)\leq\text{{C}}_{1}}\sup_{\lambda\in\mathbb{R}}$
$\displaystyle\int\varphi_{1}^{\circ}\Big{(}f_{avg}(x)+\lambda\Big{)}\mathrm{d}\alpha(x)$
(88) $\displaystyle+$
$\displaystyle\int\varphi_{2}^{\circ}\Big{(}g_{avg}(y)-\lambda\Big{)}\mathrm{d}\beta(y),$
(89)

For smooth and strictly concave $\varphi^{\circ}$, there exists a unique
$\lambda^{\star}(f_{avg},g_{avg})$ attaining the supremum. Furthermore, one
can apply the enveloppe theorem and differentiate under the integral sign
(since the support is compact). Consider perturbations
$(r_{\theta}(\cdot),s_{\theta}(\cdot))$ of
$(f_{\theta}(\cdot),g_{\theta}(\cdot))$. Write

$\displaystyle
r_{avg}(x)=\int_{\mathbb{S}^{d-1}}r_{\theta}(\theta^{\star}(x))\mathrm{d}\hat{\bm{\sigma}}_{K}(\theta),$
$\displaystyle
s_{avg}(y)=\int_{\mathbb{S}^{d-1}}s_{\theta}(\theta^{\star}(y))\mathrm{d}\hat{\bm{\sigma}}_{K}(\theta).$

Given that
$\varphi_{1}^{\circ}(f_{avg}+r_{avg})=\varphi_{1}^{\circ}(f_{avg})+r_{avg}\nabla\varphi_{1}^{\circ}(f_{avg})+o(||r_{avg}||_{\infty})$,
the first order variation reads

$\displaystyle\int
r_{avg}(x)\nabla\varphi_{1}^{\circ}\Big{(}f_{avg}(x)+\lambda^{\star}(f_{avg},g_{avg})\Big{)}\mathrm{d}\alpha(x)$
(90) $\displaystyle+$ $\displaystyle\int
s_{avg}(y)\nabla\varphi_{2}^{\circ}\Big{(}g_{avg}(y)-\lambda^{\star}(f_{avg},g_{avg})\Big{)}\mathrm{d}\beta(y).$
(91)

Then we define

$\displaystyle\bar{\alpha}=\nabla\varphi_{1}^{\circ}(f_{avg}+\lambda^{\star}(f_{avg},g_{avg}))\alpha,$
$\displaystyle\bar{\beta}=\nabla\varphi_{2}^{\circ}(g_{avg}-\lambda^{\star}(f_{avg},g_{avg}))\beta,$

such that the first order variation reads

$\displaystyle\int r_{avg}(x)\mathrm{d}\bar{\alpha}(x)+\int
s_{avg}(y)\mathrm{d}\bar{\beta}(y).$ (92)

One can then explicit the definition of $(r_{avg},s_{avg})$, such that it
reads

$\displaystyle\int_{\mathbb{S}^{d-1}}\int
r_{\theta}(\theta^{\star}(x))\mathrm{d}\bar{\alpha}(x)+\int_{\mathbb{S}^{d-1}}\int
s_{\theta}(\theta^{\star}(y))\mathrm{d}\bar{\beta}(y)$ (93) $\displaystyle=$
$\displaystyle\int_{\mathbb{S}^{d-1}}\int
r_{\theta}\mathrm{d}\theta^{\star}_{\sharp}\bar{\alpha}(x)+\int_{\mathbb{S}^{d-1}}\int
s_{\theta}\mathrm{d}\theta^{\star}_{\sharp}\bar{\beta}(y).$ (94)

By optimizing the above over the constraint set $\\{r_{\theta}\oplus
s_{\theta}\leq\text{{C}}_{1}\\}$, we identify the computation of
$\text{{SOT}}(\bar{\alpha},\bar{\beta})$, which concludes the proof. ∎

Since Proposition 7 involves potentials averaged over $\bm{\sigma}$, we thus
need to define the AvgPot routine detailed below.

Algorithm 5 – $\texttt{AvgPot}(f_{\theta})$

Input: sliced potentials $(f_{\theta})$ with $(\theta_{k})_{k}^{K}$
Output: Averaged potential $f_{avg}$ as in Proposition 7

Average $f_{avg}=\tfrac{1}{K}\sum_{k=1}^{K}f_{\theta}$

###  B.4 Implementation of Sliced OT to return dual potentials

Recall from Section 4, Algorithms 1 and 2 and more precisely, Propositions 6
and 7, that FW linear oracle is a sliced OT program, i.e. a set of OT problems
computed between univariate distributions of $\mathcal{M}_{+}(\mathbb{R})$.
Therefore, a key building block of our algorithm is to compute the loss and
dual variables of these univariate OT problems. We explain below how one can
compute the sliced OT loss and dual potentials. The computation of the loss
consists in implementing closed formulas of OT between univariate
distributions, as detailed in [57, Proposition 2.17]. More precisely, when
$\text{{C}}_{1}(x,y)=|x-y|^{p}$ and $(\mu,\nu)\in\mathcal{M}_{+}(\mathbb{R})$,
then

$\displaystyle\text{{OT}}(\mu,\nu)=\int_{0}^{1}|F^{[-1]}_{\mu}(t)-F^{[-1]}_{\nu}(t)|^{p}\mathrm{d}t,$
(95)

where $F^{[-1]}_{\mu}$ denotes the inverse cumulative distribution function
(ICDF) of $\mu$.

Algorithm 6 – $\texttt{SlicedOTLoss}(\alpha,\beta,\\{\theta\\},p)$

Input: $\alpha$, $\beta$, projections $\\{\theta\\}$, exponent $p$
Output:
$\text{{OT}}(\theta^{\star}_{\sharp}\alpha,\theta^{\star}_{\sharp}\beta)$ as
in eq. (95)

for $\theta\in\\{\theta\\}$ do

Project support of $\theta^{\star}_{\sharp}\alpha$ and
$\theta^{\star}_{\sharp}\beta$

Sort weights of $(\theta^{\star}_{\sharp}\alpha,\theta^{\star}_{\sharp}\beta)$
and support $(\theta^{\star}(x))$, $(\theta^{\star}(y))$ s.t. support is non-
decreasing

Compute ICDF of $\theta^{\star}_{\sharp}\alpha$ and
$\theta^{\star}_{\sharp}\beta$

Compute
$\text{{OT}}(\theta^{\star}_{\sharp}\alpha,\theta^{\star}_{\sharp}\beta)$ as
in eq. (95) with exponent $p$

end for

To compute dual potentials using backpropagation, one computes the sliced OT
losses (using Algorithm 6) then calls the backpropagation w.r.t to inputs
$(\alpha,\beta)$, because their gradients are optimal dual potentials [57,
Proposition 7.17]. We describe this procedure in Algorithm 7.

Algorithm 7 –
$\texttt{SlicedOTPotentialsBackprop}(\alpha,\beta,\\{\theta\\},p)$

Input: $\alpha$, $\beta$, projections $\\{\theta\\}$, exponent $p$
Output: Dual potentials $(f_{\theta},g_{\theta})$ solving
$\text{{OT}}(\theta^{\star}_{\sharp}\alpha,\theta^{\star}_{\sharp}\beta)$

Enable gradients w.r.t.
$(\theta^{\star}_{\sharp}\alpha,\theta^{\star}_{\sharp}\beta)$

Call $\texttt{SlicedOTLoss}(\alpha,\beta,\\{\theta\\},p)$

Sum (but do not average) losses
$\mathcal{L}=\sum_{\theta}\text{{OT}}(\theta^{\star}_{\sharp}\alpha,\theta^{\star}_{\sharp}\beta)$.

Backpropagate $\mathcal{L}$ w.r.t. $(\alpha,\beta)$

Return $(f_{\theta},g_{\theta})$ as gradients of $\mathcal{L}$ w.r.t.
$(\alpha,\beta)$.

The implementation of the dual potentials using 1D closed forms relies on the
north-west corner rule principle, which can be vectorized in PyTorch in order
to be computed in parallel. The contribution of our implementation thus
consists in making such algorithm GPU-compatible and allowing for a parallel
computation for every slice simultaneously. We stress that this constitutes a
non-trivial piece of code, and we refer the interested reader to the code in
our supplementary material for more details on the implementation.

###  B.5 Output optimal sliced marginals

In all our algorithms, we focus on dual formulations of SUOT and USOT, which
optimize the dual potentials. However, one might want the output variables of
the primal formulation (See Definition 3). In particular, the marginals of
optimal transport plans are interested because they are interpreted as
normalized versions of inputs $(\alpha,\beta)$ where geometric outliers have
been removed. We detail where this interpretation comes from in the setting of
UOT, and then give how it is adapted to SUOT and USOT. In particular, we
justify that the Norm routine suffices to compute them.

##### Case of UOT.

We focus on the $\text{D}_{\varphi_{i}}=\rho_{i}\text{{KL}}$. As per [4,
Equation 4.21], we have at optimality that the optimal transport $\pi^{\star}$
plan solving $\text{{UOT}}(\alpha,\beta)$ as in Equation (2) has marginals
$(\pi_{1}^{\star},\pi_{2}^{\star})$ which read
$\pi_{1}^{\star}=e^{-f^{\star}/\rho_{1}}\alpha$ and
$\pi_{2}^{\star}=e^{-g^{\star}/\rho_{2}}\beta$, where $(f^{\star},g^{\star})$
are the optimal dual potentials solving Equation (3). Since on
$\mathrm{supp}(\pi^{\star})$ one also has
$f^{\star}(x)+g^{\star}(y)=\text{{C}}_{d}(x,y)$, if the transportation cost
$\text{{C}}_{d}(x,y)$ is large (i.e. we are matching a geometric outlier), so
are $f^{\star}(x)$ and $g^{\star}(y)$, and eventually the weights
$\pi_{1}^{\star}(x)$ and $\pi_{2}^{\star}(y)$ are small, hence the
interpretation of the geometric normalization of the measures. Note that in
that case, one obtain $(\pi_{1}^{\star},\pi_{2}^{\star})$ by calling
$\texttt{Norm}(\alpha,\beta,f^{\star},g^{\star},\rho_{1},\rho_{2})$.

##### Case of SUOT.

Since $\text{{SUOT}}(\alpha,\beta)$ consists in integrating
$\text{{UOT}}(\theta^{\star}_{\sharp}\alpha,\theta^{\star}_{\sharp}\beta)$
w.r.t. $\bm{\sigma}$, it shares many similarities with UOT. For any $\theta$,
we consider $\pi_{\theta}$ and $(f_{\theta},g_{\theta})$ solving the primal
and dual formulation of
$\text{{UOT}}(\theta^{\star}_{\sharp}\alpha,\theta^{\star}_{\sharp}\beta)$.
The marginals of $\pi_{\theta}$ are thus given by
$(e^{-f_{\theta}/\rho_{1}}\alpha,e^{-g_{\theta}/\rho_{2}}\beta)$. In
particular, we retrieve the observation made in Figure 1 that the optimal
marginals change for each $\theta$. In that case we call for each $\theta$ the
routine $\texttt{Norm}(\alpha,\beta,f_{\theta},g_{\theta},\rho_{1},\rho_{2})$.

##### Case of USOT.

Recall that the optimal marginals $(\pi_{1},\pi_{2})$ in
$\text{{USOT}}(\alpha,\beta)$ do not depend on $\theta$, contrary to
$\text{{SUOT}}(\alpha,\beta)$. Leveraging the dual formulation of Theorem 5,
and looking at the Lagrangian which is defined in the proof of Theorem 5 (see
Appendix A.7), we have the optimality condition that
$\pi_{1}=e^{-f_{avg}/\rho_{1}}\alpha$ and
$\pi_{2}=e^{-g_{avg}/\rho_{2}}\beta$. Thus in that case, calling
$\texttt{Norm}(\alpha,\beta,f_{avg},g_{avg},\rho_{1},\rho_{2})$ yields the
desired marginals.

###  B.6 Convergence of Frank-Wolfe iterations: Empirical analysis

We display below an experiment on synthetic dataset to illustrate the
convergence of Frank-Wolfe iterations. We also provide insights on the number
of iterations that yields a reasonable approximation: a few iterations
suffices in our practical settings, typically $F=20$.

The results are displayed in Figure 5. We consider the empirical distributions
$(\alpha,\beta)$ computed over respectively, $N=400$ and $M=500$ samples over
the unit hypercube $[0,1]^{d}$, $d=10$. Moreover, $\beta$ is slightly shifted
by a vector of uniform coordinates $0.5\times\mathbf{1}_{d}$. We choose
$\rho=1$ and report the estimation of $\text{{SUOT}}(\alpha,\beta)$ and
$\text{{USOT}}(\alpha,\beta)$ through Frank-Wolfe iterations. We estimate the
true values by running $F=5000$ iterations, and display the difference between
the estimated score and the ’true’ values. Figure 5 shows that numerical
precision is reached in a few tens of iterations. As learning tasks do not
usually require an estimation of losses up to numerical precision, we think
that it is hence reasonable to take $F\approx 20$ in numerical applications.

Figure 5: $|\text{{SUOT}}(\alpha,\beta)-\widehat{\text{{SUOT}}}_{t}|$ and
$|\text{{USOT}}(\alpha,\beta)-\widehat{\text{{USOT}}}_{t}|$ against iteration
$t$, where $\widehat{\text{{SUOT}}}_{t},\widehat{\text{{USOT}}}_{t}$ are the
estimated SUOT, USOT using $t$ FW iterations. Plots are in log-scale. All
figures are issued from the same run, but zoomed on a subset of first
iterations: _(left)_ 1000 iterations of FW, _(middle)_ 200 iterations,
_(right)_ 20 iterations.

##  Appendix C Additional details on Section 5

###  C.1 Document classification: Technical details and additional results

####  C.1.1 Datasets

Table 2: Dataset characteristics. | BBCSport | Movies | Goodreads genre | Goodreads like
---|---|---|---|---
Doc | 737 | 2000 | 1003 | 1003
Train | 517 | 1500 | 752 | 752
Test | 220 | 500 | 251 | 251
Classes | 5 | 2 | 8 | 2
Mean words by doc | $116\pm 54$ | $182\pm 65$ | $1491\pm 538$ | $1491\pm 538$
Median words by doc | 104 | 175 | 1518 | 1518
Max words by doc | 469 | 577 | 3499 | 3499

We sum up the statistics of the different datasets in Table 2.

##### BBCSport.

The BBCSport dataset contains articles between 2004 and 2005, and is composed
of 5 classes. We average over the 5 same train/test split of [44]. The dataset
can be found in https://github.com/mkusner/wmd/tree/master.

##### Movie Reviews.

The movie reviews dataset is composed of 1000 positive and 1000 negative
reviews. We take five different random 75/25 train/test split. The data can be
found in http://www.cs.cornell.edu/people/pabo/movie-review-data/.

##### Goodreads.

This dataset, proposed in [47], and which can be found at
https://ritual.uh.edu/multi_task_book_success_2017/, is composed of 1003 books
from 8 genres. A first possible classification task is to predict the genre. A
second task is to predict the likability, which is a binary task where a book
is said to have success if it has an average rating $\geq 3.5$ on the website
Goodreads (https://www.goodreads.com). The five train/test split are randomly
drawn with 75/25 proportions.

####  C.1.2 Technical Details

All documents are embedded with the Word2Vec model [45] in dimension $d=300$.
The embedding can be found in
https://drive.google.com/file/d/0B7XkCwpI5KDYNlNUTTlSS21pQmM/view?resourcekey=0-wjGZdNAUop6WykTtMip30g.

In this experiment, we report the results averaged over 5 random train/test
split. For discrepancies which are approximated using random projections, we
additionally average the results over 3 different computations, and we report
this standard deviation in Figure 3. Furthermore, we always use 500
projections to approximate the sliced discrepancies. For Frank-Wolfe based
methods, we use 10 iterations, which we found to be enough to have a good
accuracy. We added an ablation of these two hyperparameters in Figure 7. We
report the results obtained with the best $\rho$ for USOT and SUOT computed
among a grid $\rho\in\\{10^{-4},5\cdot 10^{-4},10^{-3},5\cdot
10^{-3},10^{-2},10^{-1},1\\}$. For USOT, the best $\rho$ is consistently
$5\cdot 10^{-3}$ for the Movies and Goodreads datasets, and $5\cdot 10^{-4}$
for the BBCSport dataset. For SUOT, the best $\rho$ obtained was $0.01$ for
the BBCSport dataset, $1.0$ for the movies dataset and $0.5$ for the goodreads
dataset. For UOT, we used $\rho=1.0$ on the BBCSport dataset. For the movies
dataset, the best $\rho$ obtained on a subset was 50, but it took an
unreasonable amount of time to run on the full dataset as the runtime
increases with $\rho$ (see [48, Figure 3]). On the goodreads dataset, it took
too much memory on the GPU. For Sinkhorn UOT, we used $\varepsilon=0.001$ and
$\rho=0.1$ on the BBCSport and Goodreads datasets, and $\varepsilon=0.01$ on
the Movies dataset. For each method, the number of neighbors used for the k-NN
method is obtained via cross-validation.

####  C.1.3 Additional experiments

Figure 6: Runtime on the BBCSport dataset _(left)_ and on the Goodreads
dataset _(right)_.

##### Runtime.

We report in Figure 6 the runtime of computing the different discrepancies
between each pair of documents. On the BBCSport dataset, the documents have in
average 116 words, thus the main bottleneck is the projection step for sliced
OT methods. Hence, we observe that OT runs slightly faster than SOT and the
sliced unbalanced counterparts. Goodreads is a dataset with larger documents,
with on average 1491 words by document. Therefore, as OT scales cubically with
the number of samples, we observe here that all sliced methods run faster than
OT, which confirms that sliced methods scale better w.r.t. the number of
samples. In this setting, we were not able to compute UOT with the POT
implementation in a reasonable time. Computations have been performed with a
NVIDIA A100 GPU.

##### Ablations.

Figure 7: Ablation on BBCSport of the number of projections _(left)_ and of
the number of Frank-Wolfe iterations _(right)_.

We plot in Figure 7 accuracy as a function of the number of projections and
the number of iterations of the Frank-Wolfe algorithm. We averaged the
accuracy obtained with the same setting described in Section C.1.2, with
varying number of projections $K\in\\{4,10,21,46,100,215,464,1000\\}$ and
number of FW iterations $F\in\\{1,2,3,4,5,10,15,20\\}$. Regarding the
hyperparameter $\rho$, we selected the one returning the best accuracy, i.e.
$\rho=5\cdot 10^{-4}$ for USOT and $\rho=10^{-2}$ for SUOT.

###  C.2 Unbalanced sliced Wasserstein barycenters

We define below the formulation of the USOT barycenter which was used in the
experiments of Figure 4 to average predictions of geophysical data. We then
detail how we computed it.

######  Definition 4.

Consider a set of measures
$(\alpha_{1},\ldots,\alpha_{B})\in\mathcal{M}_{+}(\mathbb{R}^{d})^{B}$, and a
set of non-negative coefficients $(\omega_{1},\ldots,\omega_{B})\geq 0$ such
that $\sum_{b=1}^{B}\omega_{b}=1$. We define the barycenter problem (in the KL
setting) as

$\displaystyle\mathcal{B}((\alpha_{b})_{b},(\omega_{b})_{b})$
$\displaystyle\triangleq\inf_{\beta\in\mathcal{P}(\mathbb{R}^{d})}\sum_{b=1}^{B}\omega_{b}\text{{USOT}}(\alpha_{b},\beta),$
(96)
$\displaystyle=\inf_{\beta\in\mathcal{P}(\mathbb{R}^{d})}\sum_{b=1}^{B}\inf_{(\pi_{b,1},\pi_{b,2})}\text{{SOT}}(\pi_{b,1},\pi_{b,2})+\rho_{1}\text{{KL}}(\pi_{b,1}|\alpha_{b})+\rho_{2}\text{{KL}}(\pi_{b,2}|\beta),$
(97)

where $\mathcal{P}(\mathbb{R}^{d})$ denotes the set of probability measures.

To compute the barycenter, we aggregate several building blocks. First, since
we consider that the barycenter $\beta\in\mathcal{P}(\mathbb{R}^{d})$ is a
probability, we perform mirror descent as in [65, 66]. More precisely, we use
a Nesterov accelerated version of mirror descent. We also tried projected
gradient descent, but it did not yield consistent outputs (due to convergence
speed [65]). Second, we use a Stochastic-USOT version (see Section 4), i.e. we
sample new projections at each iteration of the barycenter update (but not a
each iteration of the FW subroutines in Algorithm 2). This procedure is
described in Algorithm 8.

Algorithm 8 –
$\texttt{Barycenter}((\alpha_{b})_{b},(\omega_{b})_{b},\rho_{1},\rho_{2},lr)$

Input: measures $(\alpha_{b})_{b}$, weights $(\omega_{b})_{b}$, $\rho_{1}$,
$\rho_{2}$, learning rate $lr$, FW iter $F$
Output: Optimal barycenter $\beta$ of Equation (96)

$t\leftarrow 1$

Init $(\beta,\tilde{\beta},\hat{\beta})$ as uniform distribution over a grid

while not converged do do

$\gamma\leftarrow\tfrac{2}{(t+1)}$,

$\beta\rightarrow(1-\gamma)\hat{\beta}+\gamma\tilde{\beta}$

Sample projections $(\theta_{k})_{k=1}^{K}$

Compute $\mathcal{B}((\alpha_{b})_{b},(\omega_{b})_{b})$ by calling
$\text{{USOT}}(\alpha_{b},\beta,F,(\theta_{k})_{k=1}^{K},\rho_{1},\rho_{2})$
in Algorithm 2 for each $b$

Compute $g$ as the gradient of
$\mathcal{B}((\alpha_{b})_{b},(\omega_{b})_{b})$ w.r.t. variable $\beta$

$\tilde{\beta}\leftarrow\exp(-lr\times\gamma^{-1}\times g)\beta$

$\tilde{\beta}\leftarrow\tilde{\beta}/m(\tilde{\beta})$

$\hat{\beta}\leftarrow(1-\gamma)\hat{\beta}+\gamma\tilde{\beta}$

$t\leftarrow t+1$

end while

We illustrate this algorithm with several examples of interpolation in Figure
8. We propose to compute an interpolation between two measures located on a
fixed grid of size $200\times 200$ with different values of $\rho_{i}$ in
$\text{D}_{\varphi_{i}}=\rho_{i}\text{{KL}}$. For illustration purposes, we
construct the source distribution as a mixture of two Gaussians with a small
and a larger mode, and the target distribution as a single Gaussian. Those
distributions are normalized over the grid such that both total norms are
equal to one (which is not required by our unbalanced sliced variants but
grants more interpretability and possible comparisons with SOT). Figure 8(a)
shows the result of the interpolation at three timestamps ($t=0.25,\;0.5$ and
$0.75$) of a SOT interpolation (within this setting, $\omega_{1}=1-t$ and
$\omega_{2}=t$). As expected, the two modes of the source distribution are
transported over the target one. We verify in Figure 8(b) that for a large
value of $\rho_{1}=\rho_{2}=100$, the USOT interpolation behaves similarly as
SOT, as expected from the theory. When $\rho_{1}=\rho_{2}=0.01$, the smaller
mode is not moved during the interpolation, whereas the larger one is
stretched toward the target (Figure 8(c)). Finally, in Figure 8(d), an
asymmetric configuration of $\rho_{1}=0.01$ and $\rho_{2}=100$ allows to get
an interpolation when only the big mode of the source distribution is
displaced toward the target. In all those cases, the mirror-descent algorithm
8 is run for $500$ iterations. Even for a large grid of $200\times 200$, those
different results are obtained in a $2-3$ minutes on a commodity GPU, while
the OT or UOT barycenters are untractable with a limited computational budget.

(a) SOT barycenters

(b) USOT barycenters with $\rho_{1}=\rho_{2}=100$

(c) USOT barycenters with $\rho_{1}=\rho_{2}=0.01$

(d) USOT barycenters with $\rho_{1}=0.01$ and $\rho_{2}=100$

Figure 8: Interpolation with USOT as a barycenter computation. We compare
different interpolations using SOT or USOT with different settings for the
$\rho$ values


11institutetext: Institute of Astronomy, KU Leuven, Celestijnenlaan 200D, 3001
Leuven, Belgium 22institutetext: Department of Astrophysics/IMAPP, Radboud
University, PO Box 9010, 6500 GL, The Netherlands 33institutetext: DAMTP,
Centre for Mathematical Sciences, University of Cambridge, Wilberforce Road,
Cambridge CB3 0WA, United Kingdom 44institutetext: SRON, Netherlands Institute
for Space Research, Niels Bohrweg 4, 2333 CA Leiden, The Netherlands

# The Astrophysical Gravitational Wave Background in the mHz band is likely
dominated by White Dwarf binaries

Seppe Staelens , E-mail<EMAIL_ADDRESS>Gijs Nelemans 22 1 1 4 4

(Received ; accepted )

###### Abstract

Context. The Astrophysical Gravitational Wave Background (AGWB) is a
collective signal of astrophysical gravitational wave sources and is dominated
by compact binaries. Its measurement is one of the science goals of current
and future gravitational wave detectors.

Aims. We aim to determine what population of compact binaries dominates the
AGWB in the mHz band.

Methods. We revisit and update earlier work by Farmer & Phinney (2003) to
model the astrophysical gravitational wave background sourced by extragalactic
white dwarf binaries in the mHz frequency band.We calculate the signal using a
single-metallicity model for the white dwarf population in the Universe using
a global star formation history.

Results. We estimate the white dwarf AGWB amplitude to be a factor $\sim 2$
larger than the earlier estimate and find that the overall shape of the white
dwarf AGWB is well fitted by a broken power law.

Conclusions. We compare the results to the present-day best estimates for the
background due to black hole and neutron star binaries, and find that the
white dwarf component likely dominates in the mHz band. We provide an order of
magnitude estimate that explains this hierarchy, and comment on the
implications for future missions that aim to detect the AGWB. The black hole
AGWB may only be detectable at high frequency. We outline several improvements
that can be made to our estimate, but this is unlikely to change our main
conclusion that the white dwarf AGWB dominates in the mHz band.

###### Key Words.:

gravitational waves – binaries: close – white dwarfs – black holes

##  1 Introduction

Apart from individual Gravitational Wave (GW) sources, such as the black hole
(BH) and neutron star (NS) mergers detected by the LIGO/Virgo/Kagra (LVK)
collaboration (Abbott et al., 2021) in the relatively nearby Universe, there
is an isotropic stochastic signal of all GW emitting sources in the Universe
(e.g. Schneider et al., 2010; Christensen, 2019). In particular, all compact
binaries that are evolving due to the emission of GW form a signal with a
specific spectral shape (Phinney, 2001) and this is generally referred to as
the _Astrophysical GW Background_ (AGWB, to distinguish it from potential
stochastic GW signals originating in the Early Universe, e.g. Schneider et
al., 2001; Farmer & Phinney, 2003; Kowalska-Leszczynska et al., 2015; Amaro-
Seoane et al., 2023). It is a broadband signal that extends from the high-
frequency band covered by the LVK detectors to the mHz band covered by future
space GW detectors such as LISA (Amaro-Seoane et al., 2017), Taiji (Luo et
al., 2021) or TianQin (Luo et al., 2016).

Although the AGWB is a collective signal, it provides information about the
global properties of the population of compact binaries that are not
individually detectable because they are too far away for their intrinsic
luminosity. Following the discovery that the binary BH (BBH) population is
both larger and contains more massive BH than thought before (see Abbott et
al., 2016), several investigations calculated the expected AGWB due to BBH
(and binary NS, BNS, e.g. Dvorkin et al., 2016; Bavera et al., 2022; Babak et
al., 2023; Lehoucq et al., 2023). These studies concluded that it would be
detectable by the LISA mission and future ground based detectors such as
Einstein Telescope (ET, Punturo et al., 2010) and Cosmic Explorer (CE, Reitze
et al., 2019).

Detailed understanding of the AGWB is also necessary to fully characterise and
remove it from the data to uncover underlying signals, such as GWs from the
Early Universe that are an increasingly popular research topic (see e.g.
Caprini et al., 2009; Binétruy et al., 2012; Auclair et al., 2023), or GWs
from faint individual sources.

In the mHz band, apart from BH and NS binaries there is also a contribution
from white dwarf (WD) and other binaries. Farmer & Phinney (2003) performed a
detailed calculation of the AGWB excluding BH and NS systems but did not
discuss in detail if it would be detectable by the LISA mission, while
Schneider et al. (2001) concluded that the WD AGWB would dominate. Note that
due to the much larger size of WD and other stars, they do not contribute to
the AGWB above frequencies of about 100 mHz (Farmer & Phinney, 2003). We also
stress that this WD AGWB should not be confused with the stochastic signal of
unresolved WD binaries within our Milky Way, referred to as the _GW
foreground_ , that although stochastic is far from isotropic (e.g. Edlund et
al., 2005). Given the conclusion that the BH AGWB can (easily) be detected by
the LISA mission (e.g. Caprini et al., 2016; Babak et al., 2023) it is time to
reinvestigate the relative importance of the WD AGWB.

In this paper we will compare the BH AGWB as derived from the LVK measurements
with a new estimate of the WD AGWB. In section 2 we describe the methods we
use, in section 3 we show the results and in section 4 we discuss the results,
their limitations and conclude. We use a cosmology with parameters from Planck
Collaboration et al. (2020).

##  2 Methods

The AGWB is described in terms of the dimensionless energy density spectrum

$\Omega_{\text{GW}}(f_{r})=\frac{1}{\rho_{c}c^{2}}\frac{\text{d}\mathcal{E}_{\text{GW}}}{\text{d}\ln
f_{r}}\,,$ (1)

where $\rho_{c}=\frac{3H_{0}^{2}}{8\pi G}$ is the critical density of the
Universe, $\mathcal{E}_{\text{GW}}$ is the present-day energy density in GWs
and $f_{r}$ is the received GW frequency. In an isotropic and homogeneous
Universe, the present-day GW energy density due to a population of compact
binaries in the inspiral phase is given by (Phinney, 2001)

$\displaystyle\mathcal{E}_{\text{GW}}$
$\displaystyle=\int_{0}^{\infty}\frac{\text{d}f_{r}}{f_{r}}f_{r}\int_{0}^{\infty}\text{d}z\,N(z)\frac{\text{d}E_{\text{GW}}}{\text{d}f_{e}}\,.$
(2)

In this expression, $N(z)$ is the number density of sources with energy
spectrum $\frac{\text{d}E_{\text{GW}}}{\text{d}f_{e}}$. The latter can be
approximated by the leading quadrupole order of the emitted GW radiation
(Hawking & Israel, 1987)

$\frac{\text{d}E_{\text{GW}}}{\text{d}f_{e}}=\frac{\pi^{2/3}}{3}G^{2/3}\mathcal{M}^{5/3}f_{e}^{-1/3}\qquad\text{
for }f_{\text{min}}<f_{e}<f_{\text{max}}\,.$ (3)

The frequency of the GW is related to the orbital frequency $\nu$ as
$f_{e}=2\nu$. Note that $f_{e}$ is the GW frequency in the rest frame of the
emitting source, which is redshifted as $f_{e}=f_{r}(1+z)$ by the time the
signal is detected. The limiting frequencies depend on the source properties:
the lower limit $f_{\text{min}}$ is set by the separation of the system after
circularisation, and is such that the final merger is smaller than the Hubble
time. The frequency at which the source emits increases over time as (Farmer &
Phinney, 2003)

$\dot{\nu}\equiv\frac{\text{d}\nu}{\text{d}t}=K\nu^{11/3}\,,$ (4)

until the binary components come into Roche lobe contact, or the evolution
becomes dominated by tidal forces. This sets a minimal orbital separation, and
therefore a maximal frequency $f_{\text{max}}$, for the binary. Equation (4)
can be integrated to find

$\nu(t)^{-8/3}-\nu_{0}^{-8/3}=\frac{8K(t_{0}-t)}{3}\,,$ (5)

where $t_{0}$ and $\nu_{0}$ are the time and frequency of the binary at birth,
respectively. The constant in (4) depends on the chirp mass as

$K=\frac{96}{5}\left(2\pi\right)^{8/3}\left(\frac{G\mathcal{M}}{c^{3}}\right)^{5/3}\approx
3.7\cdot 10^{-6}\left(\frac{\mathcal{M}}{M_{\odot}}\right)^{5/3}\text{
s}^{5/3}\,.$ (6)

If we allow the population to have varying chirp mass, the background is given
by (Farmer & Phinney, 2003; Renzini et al., 2022)

$\Omega_{\text{GW}}^{\text{CB}}(f_{r})=\frac{(\pi
Gf_{r})^{2/3}}{3\rho_{c}c^{2}}\int_{0}^{\infty}\text{d}z\frac{1}{(1+z)^{1/3}}\int_{\mathcal{M}_{\text{min}}}^{\mathcal{M}_{\text{max}}}\text{d}\mathcal{M}\frac{\text{d}N}{\text{d}\mathcal{M}}\mathcal{M}^{5/3}\,.$
(7)

Therefore, if the integrals do not have a frequency dependence, the expected
spectral dependence of the AGWB sourced by compact binaries is

$\Omega\left(f\right)=\Omega\left(f_{\text{ref}}\right)\left(\frac{f}{f_{\text{ref}}}\right)^{2/3}\,,$
(8)

where $f_{\text{ref}}$ is a reference frequency to which the background is
normalized.

###  2.1 Estimating the AGWB

This section aims to provide an estimate of the different contributions to the
compact binary AGWB in the LISA frequency band. We first focus on the
contributions of the binary BH (BBH) and binary NS (BNS) signals, leaving out
the NS+BH population for simplicity (and it is subdominant). These
contributions are obtained by extrapolating the current best estimates in the
LVK band using (8). In reality the signal may drop off somewhat more steeply,
especially at low-frequencies (e.g. Schneider et al., 2001), so this
extrapolation is an upper bound. The current upper limit on the AGWB in the
LVK band is (Abbott et al., 2021)

$\Omega_{\text{GW}}(25\text{ Hz})\leq 3.4\cdot 10^{-9}\,$ (9)

at the $95\%$ confidence interval. The best estimates, based on current
knowledge of the compact binary population are (Abbott et al., 2023)

$\displaystyle\Omega_{\text{BBH}}(25\text{ Hz})$
$\displaystyle=5.0^{+1.4}_{-1.8}\cdot 10^{-10}\,,$ (10)
$\displaystyle\Omega_{\text{BNS}}(25\text{ Hz})$
$\displaystyle=0.6^{+1.7}_{-0.5}\cdot 10^{-10}\,.$ (11)

The WD contribution is obtained by updating the work of Farmer & Phinney
(2003). The method we followed is based on their Sec. 6. It starts from an
alternative formulation of (1):

$\Omega(f_{r})=\frac{1}{\rho_{c}c^{3}}f_{r}F_{f_{r}}\,,$ (12)

where $F_{f_{r}}=\frac{\text{d}F(f_{r})}{\text{d}f_{r}}$ is the specific flux
received from an object with specific luminosity $L_{f_{e}}$ at redshift $z$.
It is given by

$F_{f_{r}}=\frac{L_{f_{e}}}{4\pi
d_{L}(z)^{2}}\left(\frac{\text{d}f_{e}}{\text{d}f_{r}}\right)\,,$ (13)

where $d_{L}=(1+z)d_{M}$ is the luminosity distance to redshift $z$, and
$d_{M}$ is the proper motion distance. In the case of a large number of
systems, all these contributions can be added together by introducing the
specific luminosity density $l_{f_{e}}$, defined as
$\text{d}L_{f_{e}}(z)=l_{f_{e}}(z)\text{d}V(z)$. $\text{d}V(z)$ is the
comoving volume element at redshift $z$, i.e. $\text{d}V(z)=4\pi
d_{M}(z)^{2}\text{d}\chi$, with $\chi(z)$ the proper motion distance. The
specific flux (13) is then rewritten as

$F_{f_{r}}=\int_{z=0}^{\infty}\frac{l_{f_{e}}(z)}{(1+z)^{2}}\left(\frac{\text{d}f_{e}}{\text{d}f_{r}}\right)\text{d}\chi(z)\,.$
(14)

In order to evaluate this numerically, the integral over the redshift is
discretized in 20 bins over the interval $z\in[0,8]$. Furthermore, the
frequency range $[10^{-5},1]$ Hz is divided into 50 bins, and the received
flux per bin is calculated as

$\displaystyle F_{f_{r_{1}}\to f_{r_{2}}}$
$\displaystyle=\sum_{i}\int_{f_{r_{1}}(1+z_{i})}^{f_{r_{2}}(1+z_{i})}\frac{l_{f_{e}}}{(1+z_{i})^{2}}\text{d}f_{e}\Delta\chi(z_{i})\,,$
(15)

where we integrated over $f_{r}$ and performed a change of variables $f_{r}\to
f_{e}$. This expression essentially shows that we integrate over all the
source frequencies that get redshifted to the correct observed frequency bin.
The specific luminosity density is determined as the sum of the contributions
of different systems, labeled by $k$, and is equal to

$l_{f_{e}}=\sum_{k}L_{k}(f_{e})n_{k}(f_{e},z)\,.$ (16)

The luminosity is given by (Peters & Mathews, 1963)

$L_{k}(f_{e})=\frac{32\pi^{10/3}}{5}\frac{G^{7/3}}{c^{5}}\mathcal{M}_{k}^{10/3}f_{e}^{10/3}\,,$
(17)

and $n_{k}(f_{e},z)$ is the specific number density of systems of type $k$ at
redshift $z$, emitting GWs at a frequency $f_{e}$. The specific number density
scales as $n_{k}(f)=A_{k}f^{-11/3}$ (see eq. 23), where the proportionality
constant $A_{k}$ is such that

$\int_{\text{bin}}A_{k}f^{-11/3}\text{d}f=n_{\text{bin}}(k;z)\,,$ (18)

with $n_{\text{bin}}(k;z)$ the number density of systems $k$ at redshift $z$
emitting GWs with frequencies in the correct bin. The latter can be estimated
by multiplying the star formation rate per unit volume $\psi$ with the time it
takes a binary to traverse the frequency bin:

$n_{\text{bin}}(k;z)\approx\mathcal{P}_{k}\cdot\psi(z+\Delta z)\cdot\Delta
t(k;\text{bin})\,.$ (19)

The factor $\mathcal{P}_{k}$ reflects the proportion of the formed stars that
end up in a binary with properties corresponding to label $k$. The $\Delta z$
reflects an additional redshift to reflect the conditions at formation of the
binary, before evolving to the frequency at the lower end of the bin, which
happens at redshift $z$. This time delay, as well as $\Delta t$ can be
calculated using (5).
The population model and star formation rate history that are used are
explained in the next section. Finally, the energy density in the frequency
bin is calculated as

$\Omega(f_{r})\approx\frac{1}{\rho_{c}c^{3}}\frac{f_{r}F_{f_{r_{1}}\to
f_{r_{2}}}}{f_{r_{2}}-f_{r_{1}}}\,,$ (20)

where $f_{r}$ is now an appropriately chosen frequency that represents the
bin.

###  2.2 The WD population model and star formation rate history

In order to sample the population of WD in the different redshift shells, we
calculate a model for the population of double WD based on the SeBa population
synthesis code (Portegies Zwart & Verbunt, 1996; Nelemans et al., 2001b;
Toonen et al., 2012). SeBa also has also been used to make synthetic models
for the population of double WD detectable by LISA in the Galaxy and nearby
galaxies (Nelemans et al., 2001a; Korol et al., 2020; Keim et al., 2023) and
has been used in the LISA Data Challenges (Baghi, 2022). For this first
estimate of the WD AGWB, we use a single (solar) metallicity and simulate
250,000 binaries with initial masses between 0.9 and 11 $M_{\odot}$
distributed according to the IMF of Kroupa (2001), with flat mass-ratio
distribution and initial separations flat in the $\log a$. For each of the
systems the evolution is followed and we select the 14,418 systems that form a
close ($a<500R_{\odot})$ double WD. The properties with which the WD binaries
are formed (chirp mass and GW frequency) are shown in Figure 1, with the
majority of sources having low chirp mass and frequencies at formation below 1
mHz, but a tail of systems forming with frequencies above 1 mHz.

The population we obtained corresponds to $4\cdot 10^{6}M_{\odot}$ of total
star formation, taking the complete IMF and the mass dependent binary fraction
(Moe & Di Stefano, 2017) into account. We then take this population to be
representative for every redshift shell in (15). This means that we sum over
all the binaries of our population in (16), and the factor $\mathcal{P}_{k}$
in (19) can be replaced as

$n_{\text{bin}}(k;z)\approx\frac{\psi(z+\Delta z)}{4\cdot
10^{6}M_{\odot}}\cdot\Delta t(k;\text{bin})\,.$ (21)

The star formation history is taken to be the one in Madau & Dickinson (2014):

$\psi(z)=0.015\frac{(1+z)^{2.7}}{1+[(1+z)/2.9]^{5.6}}\quad M_{\odot}\text{
yr}^{-1}\text{ Mpc}^{-3}\,.$ (22)

This exhibits a peak approximately 3.5 Gyr after the Big Bang, at $z\approx
1.9$.

Figure 1: Density plot of the initial properties of the WD population: chirp
mass $\mathcal{M}$ versus GW frequency at formation.

##  3 Results

###  3.1 Order of magnitude

Before we share the results of the estimate of the AGWB, we can make an order
of magnitude estimate of the relative importance of WD, NS and BH to the AGWB,
by comparing the number density of sources as function of frequency of each of
the classes and their GW luminosity. The number density of sources can easily
be estimated from (4) via

$\frac{dN}{df}=\frac{dN}{dt}\frac{dt}{df}=\mathcal{R}/\dot{f}=\mathcal{R}\mathcal{M}^{-5/3}f^{-11/3}$
(23)

with $\mathcal{R}$ the rate at which systems arrive at a certain frequency.
For sources that are inspiralling, $\mathcal{R}$ can be approximated by the
merger rate of the sources. The contribution to the AGWB of two sets of
sources with such density is simply this number density times the signal per
source (eq. 12)

$\Omega_{\rm GW,i}\propto\frac{dN}{df}fF_{\rm
GW}=\mathcal{R}\mathcal{M}^{5/3}f^{2/3}$ (24)

where we use eqs. (13) and (17) to estimate the flux (assuming they have the
same distance distribution).

Table 1: Order of magnitude estimates of the AGWB using the merger rate per Milky Way equivalent galaxy (MWEG) and a typical chirp mass of WD, BH and NS binaries. source | $\mathcal{R}$ | $\mathcal{M}$ | $\mathcal{R}\mathcal{M}^{5/3}$
---|---|---|---
| (/MWEG/yr) | ($M_{\odot}$) | relative
WD | $1\times 10^{-2}$ | $0.3$ | $1$
BH | $3\times 10^{-6}$ | $20$ | $0.4$
NS | $1\times 10^{-5}$ | $1.2$ | $0.01$

We assume a typical chirp mass for BBH of 20 $M_{\odot}$, for BNS of 1.2
$M_{\odot}$ and for WD binaries of 0.3 $M_{\odot}$ (components of 0.3 and 0.4
$M_{\odot}$). For BBH and BNS we estimatie the merger rates in the Milky Way
by scaling from the volume rates of the LVK collaboration (Abbott et al.,
2023) as in Abadie et al. (2010). The estimate of the WD merger rate in the
Milky Way we take from Nelemans et al. (2001b). This results in the numbers
given in Table 1, showing that the WD background likely dominates, with the BH
background getting to the same order of magnitude. The NS background is likely
significantly less important.

###  3.2 The estimated AGWB

The resulting components of the AGWB in the mHz frequency band calculated
according to the methods discussed above are shown in Figure 2, with the LVK
extrapolations in green (BBH) and blue (BNS) and the calculated WD signals in
orange. It is seen that the calculated WD component is significantly larger
than the BH (and NS) component - about an order of magnitude. Furthermore, it
is also seen to be larger than the extrapolated upper limit in the LVK
frequency band, and well above the LISA sensitivity curve in the range
$10^{-3}-10^{-2}$ Hz.

The WD component displays a turnover around 10 mHz and the BH component
becomes dominant around 100 mHz. This can be explained due to the fact that
most WD binaries have $f_{\text{max}}\lesssim 100$ mHz, causing a steady
decline in number of sources that can contribute to the total background
signal. Below $\sim$0.5 mHz the WD signal also drops off, due to the fact that
most of the systems form in that region and also because the merger time for
the systems becomes larger than the Hubble time (see also Farmer & Phinney,
2003).

The WD signal seems to follow a straight power law. In order to investigate
this further we plot in Figure 3 the signal in the region above 0.5 mHz
divided by a pure $f^{2/3}$ power law normalised at 1 mHz. It is clear that
there is a small deviation from a straight power law and also that the slope
is slightly steeper than the expected $2/3$ (eq. 8). The signal is well fitted
by a broken power law

$\Omega_{\rm GW}=\left((1.5\times 10^{-9}f^{0.79})^{-1}+(9.8\times
10^{-13}f^{-1.15})^{-1}\right)^{-1}$ (25)

with smooth transition between the two by inverted summation. The bottom panel
of Figure 3 shows the residuals of the signal with respect to the fit, which
even for this simple fit remain below a few %. So, a broken power law is a
good enough approximation for most studies.

Finally, we note that the WD AGWB is dominated by the $z\leq 1$ Universe, as
one could naively expect since the signals of WD are relatively weak. This is
shown in Figure 4, and justifies the fact that we only used redshift bins in
the range $z\in[0,8]$.

Figure 2: Resulting WD component (thick salmon) of the AGWB compared to the
LVK results (upper limit to BH/NS AGWB (dashed grey) and estimates for the BBH
and BNS components in green and blue), the LISA Powerlaw Integrated
sensitivity (black Thrane & Romano, 2013; Alonso et al., 2020) and an estimate
of the Galactic foreground (brown) based on Karnesis et al. (2021). A fit to
the WD component is shown in dark red. Figure 3: Comparison of the WD AGWB
(salmon) and the broken power law fit (dark red) to a purely $f^{2/3}$ signal
(green dashed) by dividing the curves by $f^{2/3}$. The WD AGWB deviates from
the expected $2/3$ slope but the residuals (bottom) between the WD AGWB and
the fit are below 5% Figure 4: Contribution of the different redshift shells
in (15) to the WD component of the AGWB for different frequencies. The value
quoted in the legend is the central value of the redshift bin, that all have
width 0.8 here (i.e. this simulation had only 10 redshift bins). The high
redshift contributions are small for the simple model we use and the signal is
dominated by the $z\leq 1$ Universe.

##  4 Discussion and Conclusions

Given the amplitude of the WD component in Figure 2, it is expected that it
can be very well measured by LISA. Furthermore, the relative amplitudes show
that, if LISA detects an AGWB signal in the mHz regime, it is likely dominated
by the WDs. This means that it is likely hard to make statements about the BH
(and NS) population based on a measurement of the AGWB unless there is a way
to disentangle the two, or to detect the high-frequency component of the AGWB
at about 100 mHz.

In comparison to Farmer & Phinney (2003), the computation presented in this
work results in a higher signal. They found $\Omega_{\text{WD}}(1\text{
mHz})=3.57\cdot 10^{-12}$, whereas we find $\Omega_{\text{WD}}(1\text{
mHz})=6.7\cdot 10^{-12}$. There are several factors that contribute to this
difference. First, they use a star formation history that is about a factor 2
lower than what we use, at least at the low redshifts that dominate the
background. Secondly, their population synthesis code is different, which will
lead to different results. As an indication, Farmer & Phinney (2003) find a
local WD space density of $9\cdot 10^{-3}\text{pc}^{-3}$, whereas Nelemans et
al. (2001b) \- also using the SeBa code - find $19\cdot
10^{-3}\text{pc}^{-3}$. However, since 2001 the SeBa code has been
significantly updated, so this comparison should not be taken as a specific
value, but as an indication that a difference in normalisation of a factor
$\sim 2$ between codes is not surprising. Farmer & Phinney (2003) find a range
of reference values between $1\cdot 10^{-12}\lesssim\Omega_{\text{WD}}(1\text{
mHz})\lesssim 6\cdot 10^{-12}$ in their models. We defer an in depth
investigation of the variation of the shape and normalisation of the WD AGWB
to future work.

Our estimate of the WD AGWB can be improved and extended in several ways,
although this will not change our main conclusion that the WD AGWB likely
dominates over the BH AGWB. First, we can do a study that takes the
metallicity of the star formation (e.g. Chruslinska & Nelemans, 2019) into
account as done with SeBa before (e.g. van Oirschot et al., 2014; Korol et
al., 2020), since metallicity can change the binary evolution. Additionally,
there is also more and more evidence that metallicity also changes the initial
binary parameters (Badenes et al., 2018; Moe et al., 2019). In addition, there
are many uncertainties in binary evolution so a study varying these will give
an indication of the astrophysical uncertainties in the model.

A more detailed model could also look into deviations from isotropy, since a
significant fraction of the signal is found to originate in the nearby
Universe, where the assumption of isotropy is not completely correct. This
significantly complicates the modelling, but several studies showed that
anisotropies could potentially be detected and contain information on the
source population (e.g. Cusin et al., 2020). More detailed modelling will also
give more insight in the expected deviations of the signal from simple power
laws and the shape of the turn over. This then provides input to detailed
studies on the detectability of these features e.g. with LISA. Finally, we
note here that we neglected several types of binaries that will also
contribute to the AGWB, albeit at lower level than the double WD, such as
interacting double WD binaries (AM CVn stars) and binaries containing helium
stars or low-mass main sequence stars.

In conclusion, we simulated the WD AGWB using relatively little, but very
likely sufficient detail to show that it will dominate over the AGWB from BH
and NS binaries in the mHz regime. This offers the opportunity to study the WD
binary population to much larger distances, but hampers the detection of the
BH AGWB with missions such as LISA. The WD signal reaches a peak around 10 mHz
and at higher frequencies the BH AGWB will become the dominant signal.
Detectability of this transition by LISA and other mHz missions has to be
studied in detail.

###### Acknowledgements.

G.N. is supported by the Dutch science foundation NWO. S.S. is supported by
the Centre for Doctoral Training (CDT) at the University of Cambridge funded
through STFC.

## References

  * Abadie et al. (2010) Abadie, J., Abbott, B. P., Abbott, R., et al. 2010, Classical and Quantum Gravity, 27, 173001
  * Abbott et al. (2016) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016, ApJ, 818, L22
  * Abbott et al. (2021) Abbott, R., Abbott, T. D., Abraham, S., et al. 2021, Phys. Rev. D, 104, 022004
  * Abbott et al. (2021) Abbott, R., Abbott, T. D., Acernese, F., et al. 2021, arXiv e-prints, arXiv:2111.03606
  * Abbott et al. (2023) Abbott, R., Abbott, T. D., Acernese, F., et al. 2023, Physical Review X, 13, 011048
  * Alonso et al. (2020) Alonso, D., Contaldi, C. R., Cusin, G., Ferreira, P. G., & Renzini, A. I. 2020, Phys. Rev. D, 101, 124048
  * Amaro-Seoane et al. (2023) Amaro-Seoane, P., Andrews, J., Arca Sedda, M., et al. 2023, Living Reviews in Relativity, 26, 2
  * Amaro-Seoane et al. (2017) Amaro-Seoane, P., Audley, H., Babak, S., et al. 2017, arXiv e-prints, arXiv:1702.00786
  * Auclair et al. (2023) Auclair, P., Bacon, D., Baker, T., et al. 2023, Living Reviews in Relativity, 26, 5
  * Babak et al. (2023) Babak, S., Caprini, C., Figueroa, D. G., et al. 2023, J. Cosmology Astropart. Phys., 2023, 034
  * Badenes et al. (2018) Badenes, C., Mazzola, C., Thompson, T. A., et al. 2018, ApJ, 854, 147
  * Baghi (2022) Baghi, Q. 2022, arXiv e-prints, arXiv:2204.12142
  * Bavera et al. (2022) Bavera, S. S., Franciolini, G., Cusin, G., et al. 2022, A&A, 660, A26
  * Binétruy et al. (2012) Binétruy, P., Bohé, A., Caprini, C., & Dufaux, J.-F. 2012, J. Cosmology Astropart. Phys., 2012, 027
  * Caprini et al. (2009) Caprini, C., Durrer, R., & Servant, G. 2009, J. Cosmology Astropart. Phys., 2009, 024
  * Caprini et al. (2016) Caprini, C., Hindmarsh, M., Huber, S., et al. 2016, J. Cosmology Astropart. Phys., 2016, 001
  * Christensen (2019) Christensen, N. 2019, Reports on Progress in Physics, 82, 016903
  * Chruslinska & Nelemans (2019) Chruslinska, M. & Nelemans, G. 2019, MNRAS, 488, 5300
  * Cusin et al. (2020) Cusin, G., Dvorkin, I., Pitrou, C., & Uzan, J.-P. 2020, MNRAS, 493, L1
  * Dvorkin et al. (2016) Dvorkin, I., Vangioni, E., Silk, J., Uzan, J.-P., & Olive, K. A. 2016, MNRAS, 461, 3877
  * Edlund et al. (2005) Edlund, J. A., Tinto, M., Królak, A., & Nelemans, G. 2005, Phys. Rev. D, 71, 122003
  * Farmer & Phinney (2003) Farmer, A. J. & Phinney, E. S. 2003, MNRAS, 346, 1197
  * Hawking & Israel (1987) Hawking, S. W. & Israel, W. 1987, Three hundred years of gravitation (Cambridge University Press)
  * Karnesis et al. (2021) Karnesis, N., Babak, S., Pieroni, M., Cornish, N., & Littenberg, T. 2021, Phys. Rev. D, 104, 043019
  * Keim et al. (2023) Keim, M. A., Korol, V., & Rossi, E. M. 2023, MNRAS, 521, 1088
  * Korol et al. (2020) Korol, V., Toonen, S., Klein, A., et al. 2020, A&A, 638, A153
  * Kowalska-Leszczynska et al. (2015) Kowalska-Leszczynska, I., Regimbau, T., Bulik, T., Dominik, M., & Belczynski, K. 2015, A&A, 574, A58
  * Kroupa (2001) Kroupa, P. 2001, MNRAS, 322, 231
  * Lehoucq et al. (2023) Lehoucq, L., Dvorkin, I., Srinivasan, R., Pellouin, C., & Lamberts, A. 2023, MNRAS[arXiv:2306.09861]
  * Luo et al. (2016) Luo, J., Chen, L.-S., Duan, H.-Z., et al. 2016, Classical and Quantum Gravity, 33, 035010
  * Luo et al. (2021) Luo, Z., Wang, Y., Wu, Y., Hu, W., & Jin, G. 2021, Progress of Theoretical and Experimental Physics, 2021, 05A108
  * Madau & Dickinson (2014) Madau, P. & Dickinson, M. 2014, Annual Review of Astronomy and Astrophysics, 52, 415, _eprint: https://doi.org/10.1146/annurev-astro-081811-125615
  * Moe & Di Stefano (2017) Moe, M. & Di Stefano, R. 2017, ApJS, 230, 15
  * Moe et al. (2019) Moe, M., Kratter, K. M., & Badenes, C. 2019, ApJ, 875, 61
  * Nelemans et al. (2001a) Nelemans, G., Yungelson, L. R., & Portegies Zwart, S. F. 2001a, A&A, 375, 890
  * Nelemans et al. (2001b) Nelemans, G., Yungelson, L. R., Portegies Zwart, S. F., & Verbunt, F. 2001b, A&A, 365, 491
  * Peters & Mathews (1963) Peters, P. C. & Mathews, J. 1963, Physical Review, 131, 435
  * Phinney (2001) Phinney, E. S. 2001, arXiv e-prints, astro
  * Phinney (2001) Phinney, E. S. 2001, A Practical Theorem on Gravitational Wave Backgrounds, arXiv:astro-ph/0108028
  * Planck Collaboration et al. (2020) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, Astronomy & Astrophysics, 641, A6
  * Portegies Zwart & Verbunt (1996) Portegies Zwart, S. F. & Verbunt, F. 1996, A&A, 309, 179
  * Punturo et al. (2010) Punturo, M., Abernathy, M., Acernese, F., et al. 2010, Classical and Quantum Gravity, 27, 194002
  * Reitze et al. (2019) Reitze, D., Adhikari, R. X., Ballmer, S., et al. 2019, in Bulletin of the American Astronomical Society, Vol. 51, 35
  * Renzini et al. (2022) Renzini, A. I., Goncharov, B., Jenkins, A. C., & Meyers, P. M. 2022, Galaxies, 10, 34
  * Schneider et al. (2001) Schneider, R., Ferrari, V., Matarrese, S., & Portegies Zwart, S. F. 2001, MNRAS, 324, 797
  * Schneider et al. (2010) Schneider, R., Marassi, S., & Ferrari, V. 2010, Classical and Quantum Gravity, 27, 194007
  * Thrane & Romano (2013) Thrane, E. & Romano, J. D. 2013, Phys. Rev. D, 88, 124032
  * Toonen et al. (2012) Toonen, S., Nelemans, G., & Portegies Zwart, S. 2012, A&A, 546, A70
  * van Oirschot et al. (2014) van Oirschot, P., Nelemans, G., Toonen, S., et al. 2014, A&A, 569, A42