# Spinning Black Binaries in de Sitter space

Óscar J. C. Dias<EMAIL_ADDRESS>STAG research centre & Mathematical
Sciences, University of Southampton, Highfield Campus, UK  Jorge E. Santos
<EMAIL_ADDRESS>DAMTP, Centre for Mathematical Sciences, University of
Cambridge, Wilberforce Road, Cambridge CB3 0WA, UK  Benson Way<EMAIL_ADDRESS>Maryland Center for Fundamental Physics, University of Maryland, College Park,
MD 20742, USA

###### Abstract

We construct stationary, rotating black binaries in general relativity with a
positive cosmological constant. We consider identical black holes with either
aligned or anti-aligned spins. Both cases have less entropy than the
corresponding single Kerr/Schwarzschild de Sitter black hole with the same
total angular momentum and cosmological horizon entropy. Our solutions
establish continuous non-uniqueness in general relativity without matter. They
also provide initial data for the spinning binary merger problem (when orbital
angular momentum is added).

††preprint: APS/123-QED

## Introduction.

A single black hole in isolation is an unchanging, time-independent
gravitational object. But can two or more black holes exist in such a state of
stationary equilibrium? This question has been around since the early days of
general relativity. Early work by Bach and Weyl [1], later generalized by
Israel and Khan [2] found multi-black hole solutions that were static, but
contained unphysical conical singularities. Indeed, several theorems [3, 4, 5]
preclude the existence of static, asymptotically flat multi black holes. With
the addition of electric charge, a well-known multi-black hole solution was
found by Majumdar and Papapetrou [6, 7], but these solutions require that the
black holes be maximally charged, which is unrealistic. Multi black holes were
also found in more exotic situations, such as higher-dimensional spacetimes
with compact directions [8].

A natural consideration that was not included above is the effect of rotation.
While orbital rotation breaks stationarity due to the emission of
gravitational waves, spin-spin interactions [9] do not. Notably, spin-spin
interactions are a key ingredient to molecular stability [10, 11, 12, 13].
Newtonian attraction scales with distance as $1/r^{2}$, while spin-spin
interaction scales as $1/r^{3}$ and therefore might have a dominant role at
smaller distances. Nevertheless, all known asymptotically flat, stationary
black binaries are singular, even those that maximize spin-spin repulsion (see
_e.g._ [14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26]), though the
existence of singularity-free spinning binaries has not been ruled out by any
theorem. We do mention here that with the addition of a massive complex scalar
field, the scalar pressure together with the spin-spin can balance a spinning
binary of hairy black holes [27].

What about the inclusion of a positive cosmological constant ($\Lambda>0$)?
After all, this ingredient explains the accelerated expansion of the Universe
[28, 29, 30, 31, 32], at least according to the $\Lambda$-Cold Dark Matter
($\Lambda$CDM) model. With a positive cosmological constant, it was recently
found that non-rotating black binary solutions do indeed exist [33]. These
solutions are called de Sitter black binaries. These binaries exist because
the expansion from the cosmological constant balances out the mutual
attraction of the black holes, and their existence could be anticipated from
Newton-Hooke analysis [33]. However, this configuration is unstable, as any
perturbation would cause the black holes to merge or fly apart.

In this Letter, we add rotation to the de Sitter binaries found in [33], study
their properties, and discuss the prospect that spin-spin interactions can
stabilize the binaries. For simplicity, we take the individual black holes to
be identical and focus on the cases where the spin axes are aligned or anti-
aligned, which happens to also be the situations with the strongest spin-spin
interaction.

Our results are also of mathematical interest, as they establish continuous
non-uniqueness of de Sitter black holes, in contrast to the asymptotically
flat case [34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44].

## Newton-Hooke.

Some de Sitter binaries can be described by Newton-Hooke analysis [33], which
we review here with the inclusion of spin-spin interactions. We adopt
geometrized units $c=G_{N}=k_{B}=\hbar=1$.

Consider $N$ black holes with masses $m_{a}$, positions $\mathbf{x}_{a}$, and
spin vectors $\mathcal{S}_{a}$, $a=1,\ldots,N$. Including a gravitational
force, cosmological Hooke’s law, and dipolar spin-spin interactions, Newton’s
second law becomess

$\displaystyle m_{a}\frac{\mathrm{d}^{2}\mathbf{x}_{a}}{\mathrm{d}t^{2}}=$
$\displaystyle\,m_{a}\frac{\mathbf{x}_{a}}{\ell^{2}}+\nabla\Big{(}\frac{\,m_{a}\,m_{b}}{|\mathbf{r}_{ab}|}\Big{)}$
(1)
$\displaystyle+\nabla\Bigg{[}\frac{\mathcal{S}_{a}\cdot\mathcal{S}_{b}}{|\mathbf{r}_{ab}|^{3}}-\frac{3\left(\mathcal{S}_{a}\cdot\mathbf{r}_{ab}\right)\left(\mathcal{S}_{b}\cdot\mathbf{r}_{ab}\right)}{|\mathbf{r}_{ab}|^{5}}\Bigg{]}\,,$

for $b\neq a$, de Sitter radius $\ell=\sqrt{\Lambda/3}$, and
$\mathbf{r}_{ab}\equiv\mathbf{x}_{a}-\mathbf{x}_{b}$. Stationary solutions
exist when $\frac{\mathrm{d}^{2}\mathbf{x}_{a}}{\mathrm{d}t^{2}}=0$.

Consider now the case with two equal mass black holes separated by a distance
$d$, aligned along the $z$-axis, with their spins also aligned along the
$z$-axis: $m_{1}=m_{2}\equiv m$,
$\mathbf{x}_{1}=-\mathbf{x_{2}}=\frac{z}{2}\,\mathbf{e}_{z}$ and
$\mathcal{S}_{1,2}=m\,\sigma_{1,2}\,\mathbf{e}_{z}$. Take the symmetric cases
$\sigma_{2}=\gamma\sigma_{1}\equiv\gamma\sigma$ with $\gamma=1$ (repulsive
spin-spin force) or $\gamma=-1$ (attractive spin-spin force). Then (1) admits
stationary spinning binary solutions when
$\frac{d^{3}}{\ell^{3}}=\frac{2m}{\ell}\left(1-6\,\gamma\,\frac{\sigma^{2}}{d^{2}}\right)$,
which implies

$\frac{d}{\ell}\simeq\frac{1}{\left(4\pi
T_{+}\ell\right)^{\frac{1}{3}}}\left\\{1-\left[\frac{1}{3}+\frac{2\gamma}{\left(4\pi
T_{+}\ell\right)^{4/3}}\right]\frac{\Omega_{+}^{2}}{\left(4\pi
T_{+}\right)^{2}}\right\\}\,$ (2)

In deriving the above, we discarded $\mathcal{O}(\Omega_{+}^{4})$ terms. We
also expressed the mass and spin of the test body in terms of the temperature
$T_{+}$ and angular velocity $\Omega_{+}$ of a Kerr black hole, so that

$\displaystyle 2\,m=\frac{1}{2\pi
T_{+}+\sqrt{4\pi^{2}T_{+}^{2}+\Omega_{+}^{2}}}\,,$ (3)
$\displaystyle\sigma=\frac{m\,\Omega_{+}}{\sqrt{4\pi^{2}T_{+}^{2}+\Omega_{+}^{2}}}\,.$

For vanishing spin, $\Omega_{+}=0$, (2) yields $\frac{d^{3}}{\ell^{3}}=(4\pi
T_{+}\ell)^{-\frac{1}{3}}=\frac{r_{S}}{\ell}$, where $r_{S}=2m$ is the
Schwarzschild radius of the particles, in agreement with [33]. The equilibrium
condition (2) can fall within the regime of validity of Newton-Hook theory,
which occurs when $r_{S}\ll d\ll\ell$ (i.e. large $T_{+}\ell$).

With nonzero spin, $\Omega_{+}\neq 0$, an inspection of (2) for a given
(large) $T_{+}\ell$ shows that $d/\ell$ is a monotonically decreasing function
of $\Omega_{+}\ell$. Moreover, the $d(\Omega_{+})$ curve for the aligned
binary ($\gamma=1$) is always below the curve for the anti-aligned case
($\gamma=-1$). This agrees with expectations as spin-spin forces are repulsive
for aligned spins, so the black holes need to be closer apart to remain in
equilibrium (for fixed gravitational and cosmological forces). The opposite is
true for anti-aligned spins. Our numerical (anti-)aligned spinning solutions
of de Sitter general relativity will match this behaviour.

## Numerical construction.

Now we construct spinning de Sitter binaries using numerical relativity. We
use the Einstein-DeTurck method [45] (see [46, 47] for a review), which first
involves the selection of a reference metric $\bar{g}$ which has the same
causal (regularity, horizon, and asymptotic) structure and, usually, the same
symmetries as the solution we seek.

Our reference metric is based on the one used in [33], which is a carefully
chosen combination of the $\Lambda=0$ Bach-Weyl solution with two identical
black holes [1, 2, 48] and the static patch of de Sitter space. Like the
reference metric in [33], ours has a $\mathbb{Z}_{2}$ symmetry, temporal and
azimuthal Killing vectors $k=\partial/\partial t$ and
$m=\partial/\partial\phi$, two black hole horizons, and a cosmological
horizon. Our main modification is the addition of a cross term
$\mathrm{d}t\mathrm{d}\phi$ that introduces the angular velocities
$\Omega_{+}^{(i)}$ ($i=1,2$) to the black holes in a frame where the de Sitter
horizon is not rotating. The $\mathbb{Z}_{2}$ symmetry forces the spin axes of
the binaries to be aligned or anti-aligned. Further details on the reference
metric can be found in the Supplementary Material.

With a reference metric, we then solve the Einstein-DeTurck equation

$R_{ab}-\nabla_{(a}\xi_{b)}=\frac{3}{\ell^{2}}g_{ab}\,,$ (4)

where $g_{ab}$ and $R_{ab}$ are the metric and its Ricci tensor,
$\xi^{a}\equiv
g^{bc}\left[\Gamma^{a}_{bc}(g)-\Gamma^{a}_{bc}(\bar{g})\right]$, and $\Gamma$
is the Christoffel connection. Unlike the Einstein equation (without proper
gauge fixing), the equation (4) is elliptic [45, 46, 49, 47], and yields a
well-posed boundary problem with suitable boundary conditions. In our case,
boundary conditions are set by regularity and symmetry requirements. The
aligned $\Omega_{+}^{(2)}=\Omega_{+}^{(1)}$ and anti-aligned
$\Omega_{+}^{(2)}=-\Omega_{+}^{(1)}$ cases have the same boundary conditions
except at the $\mathbb{Z}_{2}$ symmetry plane, where the $g_{t\phi}$ metric
component has Neumann and Dirichlet boundary conditions, respectively.

A solution to (4) solves the Einstein equation only when $\xi=0$ (the solution
will then be in the gauge $\xi=0$). In certain cases, solutions to the
Einstein-DeTurck equation are guaranteed to satisfy $\xi=0$ [49, 50]. However,
there are no such guarantees when $\Lambda>0$. Fortunately, ellipticity
guarantees local uniqueness. That is, solutions with $\xi=0$ cannot be
arbitrarily close to those with $\xi\neq 0$, and thus we can monitor the norm
$\xi^{a}\xi_{a}$ to verify that our numerical discretization converges in the
continuum to a true Einstein solution (see Supplementary Material).

## Results.

With our symmetry requirements, the rotating de Sitter binaries form a two-
parameter family. Numerically, we take these parameters to be the horizon
temperature $T_{+}\ell$ and angular velocity $\Omega_{+}\ell$. Because
scanning the full parameter space is computationally expensive, we leave a
thorough sweep to future work. Instead we focus our attention on solutions
with a fixed temperature $T_{+}\ell$, and varying $\Omega_{+}\ell$. Note that
the Newton-Hooke regime occurs for $T_{+}\ell\gg 1$ (_i.e._ very small black
holes).

In Fig. 1, we show the Komar angular momentum $J_{+}/\ell^{2}$ of one of the
black holes versus its angular velocity $\Omega_{+}\ell$. The highest
angularity velocity we have reached is $\Omega_{+}\ell=1.72$. We see no
evidence of pathologies or singularities (see Supplementary Material for
curvature invariants), so it is possible that these binaries exist for higher
$\Omega_{+}$ (however, it cannot be ruled out that solutions stop existing at
a critical $\Omega_{+}\ell$ due to a lack of bound states).

We show both aligned and anti-aligned results, but they do not differ much. We
see that the angular momentum increases with angular velocity. If solutions
exist for higher $\Omega_{+}$, we would expect that the angular momentum will
reach a maximum, as the same occurs for constant-temperature Kerr de Sitter
black holes. A polynomial extrapolation gives a peak around $\Omega_{+}\ell=3$
111In a forthcoming publication [DSW-chargedBinary], we construct charged de
Sitter binaries, and we see a maximum charge as the electric potential
increases..

Figure 1: Black hole angular momentum of one of binary’s black hole versus its
angular velocity for binaries with $T_{+}\ell=\frac{3.75}{2\pi}$.

We now discuss the proper distance $\mathcal{P}_{\phi}^{in}$ along the $\phi$
symmetry axis between the black hole horizons as a function of
$\Omega_{+}\ell$. This is shown in Fig. 2 for a temperature
$T_{+}\ell=\frac{9}{\pi}$ that seems sufficiently large to expect qualitative
matching with Newton-Hooke expectations. We see that, in agreement with the
analysis of the Newton-Hooke relation (2), by increasing spin the black holes
move closer together with the aligned binary always having a smaller distance.

Figure 2: Proper distance along the symmetry axis between the two black holes,
$\mathcal{P}_{\phi}^{in}/\ell$, versus angular velocity $\Omega_{+}\ell$ for
$T_{+}\ell=\frac{9}{\pi}$.

We now discuss the thermodynamics of the system. In our case, the first law of
black hole mechanics reads [52]

$2T_{+}\,\mathrm{d}S_{+}+2\Omega_{+}\,\mathrm{d}J_{+}=-T_{c}\,\mathrm{d}S_{c}\,$
(5)

where $T_{c}$ is the temperature of the cosmological horizon, and $S_{c}$ is
its entropy. We checked that our data obeys this law to within $0.01\%$.

Figure 3: Projection of known asymptotically de Sitter (dS) static and
stationary solutions in the plane that displays the total black hole entropy
$S/\ell^{2}$ versus the cosmological horizon entropy $S_{c}/\ell^{2}$. The
inset plots are zoom-ins of regions of particular interest as discussed in the
text. The spinning binaries have $T_{+}\ell=\frac{3.75}{2\pi}$.

Figure 4: Contour plot showing the level sets of the gauge invariant metric
component $-g_{tt}$. The cosmological horizon is the outer solid black
semicircle. The two black hole horizons are located along $r=0$ at
$z\in[-1,-2]$ and $z\in[1,2]$, and the outer and inner axes are the dashed
black lines at $r=0$. The inset plot is a zoom around the two black holes. The
green triangle is where $-g_{tt}$ takes its maximum value and the magenta
inverted-triangle is where $-g_{tt}$ takes its minimum (negative) value.

In Fig. 3, we show the total black hole entropy $S/\ell^{2}$ against the
cosmological horizon entropy $S_{c}/\ell^{2}$ [52]. Again, we are only showing
the spinning de Sitter binaries for the constant temperature
$T\ell=\frac{3.75}{2\pi}$. For reference, we also included the corresponding
values for Schwarzschild de Sitter, Kerr de Sitter, and the static binaries
from [33]. In the Kerr case, we only show black holes with the same total
angular momentum as the aligned de Sitter binaries. The anti-aligned binaries
have zero total angular momentum, so should be compared with Schwarzschild and
the static binaries. We see that in all cases, the spinning binaries have less
entropy $S$ than the corresponding Schwarzschild or Kerr black holes, and so
do not dominate the microcanonical ensemble.

From a dynamical perspective, the cosmological horizon entropy does not
necessarily stay constant during evolution. But, both horizons satisfy a
second law, so it cannot decrease during evolution. We note that our maximal
entropy in Fig. 3 is monotonically decreasing with $S_{c}$. That is, the
maximum entropy $S_{\mathrm{max}}(S_{c})$ is larger than any entropy
$S_{\mathrm{max}}(S_{c}^{\prime})$ for $S_{c}^{\prime}>S_{c}$. This suggests
that the endpoint that maximizes the entropy is indeed the Schwarzschild de
Sitter black hole and also suggests its stability. Similar arguments can be
made in asymptotically flat space, where energy and angular momentum can be
radiated to infinity, and hence are not conserved during evolution.

## Conclusions.

We constructed the first examples of stationary, spinning black binary
solutions of general relativity with a positive cosmological constant where
the gravitational, cosmological and spin-spin interactions combine to produce
an equilibrium configuration. We have found these solutions well away from the
regime where Newton-Hooke theory is valid, but the properties of our binaries
match the Newton-Hooke expectations.

We point out that, like in [33], the available uniqueness theorems [53, 54,
55, 56] do not actually forbid the existence of these de Sitter binaries. This
is largely because such theorems make assumptions about the level sets of the
gauge-invariant metric component $-g_{tt}$. Fig. 4 shows such level sets for a
binary with $\Omega_{+}^{(2)}\ell=\Omega_{+}^{(1)}\ell=1.7$ (other solutions
are qualitatively similar). The location of the maximum of $-g_{tt}$ is marked
by a green triangle. The fact that this maximum is a point in this plot and
not a line demonstrates that the assumptions in these theorems are violated
(we refer the reader to [33] for a more detailed discussion on this matter).

For completeness, we also show the minima of $-g_{tt}$ marked by the two
magenta inverted-triangles. These minima take negative values because of the
ergoregion around each black hole (where energy and angular momentum can be
extracted from the Penrose process or superradiant scattering).

It is clear that black hole uniqueness is strongly violated in de Sitter,
unlike in flat space. Indeed, we have an infinite non-uniqueness. For example,
in the anti-aligned case, we have a whole continuous two-parameter family of
solutions with zero total angular momentum. Therefore, there is no analogous
no-hair conjecture [57] for pure gravity with a $\Lambda>0$.

Although the solutions we have found do not dominate the microcanonical
ensemble, this does not necessarily mean that they are dynamically unstable.
Because spin-spin interactions act on shorter length scales than that of the
Newtonian potential, they may provide a mechanism for stabilizing the binary
in some regions of parameter space, much like the stability of molecules [10,
11, 12, 13]. It would be interesting to see if this intuition also holds
within general relativity, which would require a proper stability analysis of
our numerical spinning binaries (or their extensions in the parameter space).

We have left a more thorough search of parameter space to future work. This
could include higher angular velocities as well as a sweep over the full
2-parameter family. It would also be interesting to drop the $\mathbb{Z}_{2}$
symmetry restrictions and consider unequal black holes with unequal spins. In
the case where one of the black holes is much smaller than the other, the
Mathisson-Papapetrou-Dixon’s formalism [58, 59, 60, 61, 62, 63, 9] can be
used. Work in this direction is underway.

Finally, we note that our solutions have spin but no orbital angular momentum.
Adding orbital angular momentum will also create quadrupole moments that emit
gravitational radiation and hence, the system will not be stationary.
Nevertheless, our solutions can provide initial data that satisfy the elliptic
constraint equations of the time evolution problem, that ultimately leads to
the black hole binary merger and associated gravitational (and
electromagnetic) wave emission.

###### Acknowledgements.

## Acknowledgments.

O.D. acknowledges financial support from the STFC “Particle Physics Grants
Panel (PPGP) 2018” Grant No. ST/T000775/1 and PPGP 2020 grant No.
ST/X000583/1. O.D.’s research was also supported in part by the YITP-ExU long-
term workshop Quantum Information, Quantum Matter and Quantum Gravity
(QIMG2023), Yukawa Institute for Theoretical Physics, Kyoto Univ., where part
of this work was completed. J. E. S. has been partially supported by STFC
consolidated grant ST/T000694/1. The authors acknowledge the use of the IRIDIS
High Performance Computing Facility, and associated support services at the
University of Southampton, for the completion of this work.

## Supplementary Material

##  Appendix A Metric Ansatz for Spinning Black Binaries in de Sitter

In the Einstein-DeTurck formalism, (4) requires a reference metric $\bar{g}$
which has the same casual (asymptotic and horizon) structure and, whenever
possible, also the same symmetries as the solution we seek. We will follow the
same procedure in [33] to obtain the reference metric, making a few
modifications to accommodate spin.

The reference metric in [33] was built from a careful combination of the Bach-
Weyl [1] (i.e., the Israel-Khan solution [2] with two black holes) and de
Sitter space. The Bach-Weyl solution is a closed form, asymptotically flat
solution of the Einstein equation that consists of two non-spinning black
holes held together by a conical singularity. The overall strategy is to
combine this singular binary with the de Sitter horizon in hopes of finding a
de Sitter solution with a conical singularity, and then remove the conical
singularity.

The Bach-Weyl solution is usually written in Weyl coordinates $(r,z)$ because
the resulting Einstein equation in this gauge can be expressed in an
integrable form. These coordinates have a “rod-structure” so that the
coordinate line $r=0$ includes both inner and outer segments of the $\phi$
symmetry axis as well as the black hole horizons. This rod-structure is
inconvenient for numerical purposes, largely because it introduces more
coordinate singularities and makes the application of boundary conditions
difficult. Instead, in [33], a transformation

$z=\frac{x\sqrt{2-x^{2}}\sqrt{(1-y^{2})^{2}+k^{2}y^{2}(2-y^{2})}}{(1-y^{2})^{2}+k^{2}x^{2}(2-x^{2})y^{2}(2-y^{2})}\;,\qquad
r=\frac{(1-x^{2})\sqrt{1-k^{2}x^{2}(2-x^{2})}y\sqrt{2-y^{2}}(1-y^{2})}{(1-y^{2})^{2}+k^{2}x^{2}(2-x^{2})y^{2}(2-y^{2})}\;,$
(6)

was used to map both horizons, as well as the outer and inner segments of the
axis into the four sides of a coordinate rectangle $(x,y)$. The Bach-Weyl
solution [1, 2] then takes the form [33]

$\mathrm{d}s^{2}=\ell^{2}\left\\{-f\mathrm{d}t^{2}+\frac{\lambda^{2}}{m^{2}\Delta_{xy}^{2}}\left[p^{2}\left(\frac{4\mathrm{d}x^{2}}{(2-x^{2})\Delta_{x}}+\frac{4\mathrm{d}y^{2}}{(2-y^{2})\Delta_{y}}\right)+y^{2}(2-y^{2})(1-y^{2})^{2}\mathrm{d}\phi^{2}\right]\right\\}\;,$
(7)

where $\ell$ is an arbitrary dimensionful length scale that we have introduced
for later use in de Sitter, and

$\displaystyle\Delta_{x}=1-k^{2}x^{2}(2-x^{2})\;,\qquad\Delta_{y}=1-(1-k^{2})y^{2}(2-y^{2})\;,\qquad\Delta_{xy}=(1-y^{2})^{2}+k^{2}x^{2}(2-x^{2})y^{2}(2-y^{2})\;,$
(8) $\displaystyle f=(1-x^{2})^{2}\Delta_{x}m^{2}\;,\quad
p=\frac{k\left(1+\sqrt{\Delta_{y}}\right)^{2}}{(1+k)^{2}}\;,\quad
m=\frac{k\Big{[}1-(1-k)y^{2}(2-y^{2})+\sqrt{\Delta_{y}}\Big{]}}{(1-k)\Delta_{x}(1-y^{2})^{2}+(k+\sqrt{\Delta_{y}})\Big{[}\Delta_{x}+(1-k)(\sqrt{\Delta_{xy}}-1)\Big{]}}\;.$

All functions $f$, $\Delta_{x}$, $\Delta_{y}$, $p$, and $m$ are smooth and
positive definite in the domain. $\Delta_{xy}$ vanishes at $(x,y)=(0,1)$
(asymptotic infinity), and is positive and smooth otherwise. The solution is
parametrized by $k\in(0,1)$.

To eventually accommodate a cosmological horizon, we also consider the polar
Weyl coordinates $(\rho,\xi)$ defined by:

$z=\rho\,\xi\sqrt{2-\xi^{2}}\;,\qquad r=\rho(1-\xi^{2})\;,$ (9)

where the Bach-Weyl solution [1, 2] takes the form [33]

$\mathrm{d}s^{2}=\ell^{2}\left\\{-f\mathrm{d}t^{2}+\frac{\lambda^{2}h}{f}\left[\mathrm{d}\rho^{2}+\rho^{2}\left(\frac{4\mathrm{d}\xi^{2}}{2-\xi^{2}}+\frac{(1-\xi^{2})^{2}}{h}\mathrm{d}\phi^{2}\right)\right]\right\\}\;,$
(10)

with

$\displaystyle
R_{\pm}=\sqrt{\rho^{2}+\frac{1}{k^{2}}\pm\frac{2}{k}\rho\,\xi\sqrt{2-\xi^{2}}}\;,\qquad
r_{\pm}=\sqrt{\rho^{2}+1\pm 2\rho\,\xi\sqrt{2-\xi^{2}}}\;,$ $\displaystyle
f=\left(\frac{k(R_{+}+r_{+})-(1-k)}{k(R_{+}+r_{+})+(1-k)}\right)\left(\frac{k(R_{-}+r_{-})-(1-k)}{k(R_{-}+r_{-})+(1-k)}\right)\;,$
$\displaystyle
h=\left(\frac{\rho^{2}+\tfrac{1}{k}[1+(1+k)\rho\,\xi\sqrt{2-\xi^{2}}]+R_{+}r_{+}}{2R_{+}r_{+}}\right)\left(\frac{\rho^{2}+\tfrac{1}{k}[1-(1+k)\rho\,\xi\sqrt{2-\xi^{2}}]+R_{-}r_{-}}{2R_{-}r_{-}}\right)$
(11)
$\displaystyle\qquad\times\left(\frac{\rho^{2}-1+r_{+}r_{-}}{\rho^{2}-\tfrac{1}{k}[1+(1-k)\rho\,\xi\sqrt{2-\xi^{2}}]+r_{+}R_{-}}\right)\left(\frac{\rho^{2}-(1/k^{2})+R_{+}R_{-}}{\rho^{2}-\tfrac{1}{k}[1-(1-k\,\rho\,\xi\sqrt{2-\xi^{2}})]+R_{+}r_{-}}\right)\;.$

Note that $h$ and $f$ approach unity when $\rho\to\infty$, where the spacetime
becomes asymptotically flat.

Now we separately consider de Sitter space, which can be written in isotropic
coordinates [33]

$\mathrm{d}s^{2}=\frac{\ell^{2}}{g_{+}^{2}}\left\\{-g_{-}^{2}\mathrm{d}t^{2}+\lambda^{2}\left[\mathrm{d}\rho^{2}+\rho^{2}\left(\frac{4\mathrm{d}\xi^{2}}{2-\xi^{2}}+(1-\xi^{2})^{2}\mathrm{d}\phi^{2}\right)\right]\right\\}\;,\qquad\hbox{with}\>\>g_{\pm}=1\pm\frac{\lambda^{2}\rho^{2}}{4},$
(12)

where $\ell$ is the de Sitter radius $\ell=\sqrt{\Lambda/3}$. In these
coordinates, the de Sitter horizon has a constant temperature of
$T_{c}=1/(2\pi)$. In (12), as in (10), $\lambda$ is a gauge parameter that
simply scales the radial coordinate $\rho$.

The form (12) for de Sitter space is suggestively similar to the Bach-Weyl
solution (10) in polar Weyl coordinates. Aside from some factors of $f$ and
$h$ (which approach unity at large $\rho$), the only differences are that de
Sitter in isotropic coordinates has an overall conformal factor of
$1/g_{+}^{2}$ and a factor of $g_{-}^{2}$ in the $\mathrm{d}t^{2}$ term whose
zero defines the de Sitter horizon. This effectively constrains the $\rho$
coordinate in our integration domain to $\rho<2/\lambda$. We will make use of
these similarities below.

Now we can form a reference metric as was done in [33]. To do so note that
from (6) and (9), we can find an explicit coordinate transformation between
the $(x,y)$ coordinates and the polar Weyl $(\rho,\xi)$ coordinates:

$\displaystyle\rho$
$\displaystyle=\frac{\sqrt{y^{2}(2-y^{2})+x^{2}(2-x^{2})(1-y^{2})^{2}}}{\sqrt{(1-y^{2})^{2}+k^{2}x^{2}(2-x^{2})y^{2}(2-y^{2})}}\;,$
(13a) $\displaystyle\xi$
$\displaystyle=\sqrt{1-\frac{(1-x^{2})y\sqrt{2-y^{2}}(1-y^{2})\sqrt{1-k^{2}x^{2}(2-x^{2})}}{\sqrt{y^{2}(2-y^{2})+x^{2}(2-x^{2})(1-y^{2})^{2}}\sqrt{(1-y^{2})^{2}+k^{2}x^{2}(2-x^{2})y^{2}(2-y^{2})}}}\;.$
(13b)

which is needed in what follows.

With these ingredients in place, we can finally describe our reference metric:

$\displaystyle\mathrm{d}s^{2}_{\mathrm{ref}}$
$\displaystyle=\frac{\ell^{2}}{g_{+}^{2}}\left\\{-fg_{-}^{2}\,F\,\mathrm{d}t^{2}+\frac{\lambda^{2}}{m^{2}\Delta_{xy}^{2}}\left[p^{2}\left(\frac{4\mathrm{d}x^{2}}{(2-x^{2})\Delta_{x}}+\frac{4\mathrm{d}y^{2}}{(2-y^{2})\Delta_{y}}\right)+y^{2}(2-y^{2})(1-y^{2})^{2}\,s\,\left(\mathrm{d}\phi-\beta\,g_{-}^{2}w\,\mathrm{d}t\right)^{2}\right]\right\\}$
$\displaystyle=\frac{\ell^{2}}{g_{+}^{2}}\left\\{-fg_{-}^{2}\,F\,\mathrm{d}t^{2}+\frac{\lambda^{2}h}{f}\left[\mathrm{d}\rho^{2}+\rho^{2}\left(\frac{4\mathrm{d}\xi^{2}}{2-\xi^{2}}+\frac{(1-\xi^{2})^{2}}{h}\,s\,\big{(}\mathrm{d}\phi-\beta\,g_{-}^{2}w\,\mathrm{d}t\big{)}^{2}\right)\right]\right\\}\;,$
(14)

with several functions that appear in theses metrics defined in (8) or (A). We
will describe the remaining functions further below. The equality between the
first and second lines of (A) (here and in the remainder of this section) is
understood to be through the coordinate transformation (13). We will use
$(x,y)$ coordinates in the region near the black holes and inner segment of
the axis, and the $(\rho,\xi)$ coordinates near the cosmological horizon. The
DeTurck metric (A) is a simple upgrade of the static DeTurck reference metric
(B5) of [33] where we have simply introduced a $g_{t\phi}$ term to accommodate
the spins of the black holes. More concretely, we have simply taken (B5) of
[33] and made the replacement $\mathrm{d}\phi^{2}\to\left(\mathrm{d}\phi-\beta
g_{-}^{2}w\,\mathrm{d}t\right)^{2}$ (which also introduces a new function $w$
defined below). This reference metric depends on four parameters,
$\lambda,k,\alpha$ and $\beta$ whose physical interpretation is given below.

We only made four changes to the Bach-Weyl (Israel-Khan) solution to arrive at
the reference metric (A). The first is the inclusion of a conformal factor
$1/g_{+}^{2}$ to enable the matching to de Sitter spacetime (12) in isotropic
coordinates in the far-region. The second is the inclusion of a function $s$
in the $\mathrm{d}\phi^{2}$ term with a conical parameter $\alpha$ that we
will adjust to remove any conical singularities. This is achieved with the
choices

$s=1-\alpha(1-y^{2})^{2}\;,\qquad\alpha=\frac{(1-k)^{2}\left(k^{2}+6k+1\right)}{(k+1)^{4}}\,.$
(15)

The third change is the inclusion of a factor of $g_{-}^{2}\,F$ in the
$\mathrm{d}t^{2}$ term which introduces a cosmological horizon at
$\rho=2/\lambda$ ($F$ is a complicated function that will be given further
below in (17)). The fourth change is the replacement
$\mathrm{d}\phi^{2}\to\left(\mathrm{d}\phi-\beta
g_{-}^{2}w\,\mathrm{d}t\right)^{2}$ where $w$ must satisfy $g_{-}^{2}w|_{x=\pm
1}=1$ but can otherwise be freely chosen. We choose

$w=\frac{1}{\left(1-x^{2}\right)^{2}+g_{-}^{2}}\,.$ (16)

This fourth change introduces an angular velocity at the horizons $x=\pm 1$
(which is given by the novel parameter $\beta$) but leaves the cosmological
horizon (where $g_{-}$ vanishes) with zero rotation.

We have some freedom to choose the function $F$, but the choice is delicate
[33]. For numerical reasons, we want $F$ to be smooth both in $(x,y)$ and
$(\rho,\xi)$ coordinates. The DeTurck method further requires that $F$
preserves the regularity of both the black hole and cosmological horizons [45,
46, 47]. In addition, to aid in finding a solution using a Newton-Raphson
algorithm, it would be convenient to also choose $F$ to match physical
expectations in certain limits. Specifically, we expect that when the
cosmological horizon is large compared to other length scales (_i.e._
$\lambda\ll 1$), the spacetime near the cosmological horizon should approach
de Sitter and the spacetime closer to the origin should be approximately
described by the Bach-Weyl solution (when $\alpha=0$ and $\beta=0$). As
detailed in [33], these requirements are achieved if we choose $F$ to be [33]

$F=\frac{G}{f+g_{-}^{2}G-fg_{-}^{2}G}\;,\qquad\hbox{with }\quad
G=\frac{\tfrac{\hat{h}}{\hat{f}}(1-x^{2})y^{2}(2-y^{2})+g_{-}^{2}}{(1-x^{2})y^{2}(2-y^{2})+g_{-}^{4}}\;,$
(17)

where $\hat{f}$ and $\hat{h}$ are any smooth, positive definite functions that
agree with $f$ and $h$, respectively at $\rho=2/\lambda$. To choose $\hat{f}$
and $\hat{h}$, we first take the expressions for $f$ and $h$ as written in
(A), and treat them as functions $f(\rho,z)$ and $h(\rho,z)$. We then set
$\hat{f}(\rho,\xi)=f(2/\lambda,\rho\xi\sqrt{2-\xi^{2}})$ and similarly for
$\hat{h}$. Note that we cannot use a choice like
$\hat{f}(\rho,\xi)=f(\tfrac{2}{\lambda},\tfrac{2}{\lambda}\xi\sqrt{2-\xi^{2}})$
as it is not smooth in the $(x,y)$ coordinates at $x=0,y=1$.

The reference metric (A) is parametrized by $\lambda$, $k$, $\alpha$ and
$\beta$. This metric describes a solution that is everywhere regular when the
conical parameter $\alpha$ to be given by (15), with black hole horizons at
$x=\pm 1$. These are Killing horizons generated by the Killing vector field
$K=\partial_{t}+\Omega_{+}\partial_{\phi}$, _i.e._ $|K|_{x=\pm 1}=0$ (where
we use the fact that $f=(1-x^{2})^{2}\Delta_{x}m^{2}$ and $g_{-}^{2}w|_{x=\pm
1}=1$). The temperature $T_{+}$ and angular velocity $\Omega_{+}$ are given by

$T_{+}=\frac{1}{2\pi}\frac{k(1+k)}{4\lambda(1-k)}\,,\qquad\Omega_{+}=\beta\,,$
(18)

where $\lambda\in(0,\infty)$ and $k\in(0,1)$. Note that because of the overall
factor of $\ell$ in the metric, $T_{+}$ and $\Omega_{+}$ as expressed here are
dimensionless. Any $k$ and $\lambda$ that give the same $T_{+}$ are physically
equivalent. Effectively, $\lambda$ is a gauge parameter that only sets an
overall scale (it fixes the cosmological horizon location to be at
$\rho=2/\lambda$) and $k$ then fixes the black hole temperature. Solutions
described by (A) also have a cosmological horizon at $\rho=2/\lambda$ (where
$g_{-}$ vanishes) with temperature $T_{c}=1/(2\pi)$ and angular velocity
$\Omega_{c}=0$. This is generated by the Killing vector field
$\zeta=\partial_{t}+\Omega_{c}\partial_{\phi}=\partial_{t}$, _i.e._
$|\zeta|_{\rho=2/\lambda}=0$. Note that, unlike the event horizon, the
cosmological horizon has vanishing angular velocity. The reference metric
therefore only has two physical parameters, $T_{+}$ and $\Omega_{+}$ (or $k$
and $\beta$).

With a reference metric, we can finally write our metric ansatz for the
spinning binary in de Sitter, which is essentially the most general ansatz
compatible with the symmetries of the problem:

$\displaystyle\mathrm{d}s^{2}$
$\displaystyle=\frac{\ell^{2}}{g_{+}^{2}}\Bigg{\\{}-fg_{-}^{2}\,F\,\mathcal{T}\,\mathrm{d}t^{2}+\frac{\lambda^{2}}{m^{2}\Delta_{xy}^{2}}\Bigg{[}p^{2}\left(\frac{4\mathcal{A}\,\mathrm{d}x^{2}}{(2-x^{2})\Delta_{x}}+\frac{4\mathcal{B}}{(2-y^{2})\Delta_{y}}\left(\mathrm{d}y-x\,(1-x^{2})\,y\,(2-y^{2})(1-y^{2})\mathcal{F}\,\mathrm{d}x\right)^{2}\right)$
$\displaystyle\hskip
150.79968pt+y^{2}(2-y^{2})(1-y^{2})^{2}\,s\,\mathcal{S}\left(\mathrm{d}\phi-\,g_{-}^{2}w\,\mathcal{W}\,\mathrm{d}t\right)^{2}\Bigg{]}\Bigg{\\}}$
$\displaystyle=\frac{\ell^{2}}{g_{+}^{2}}\Bigg{\\{}-fg_{-}^{2}\,F\,\widetilde{\mathcal{T}}\,\mathrm{d}t^{2}+\frac{\lambda^{2}h}{f}\Bigg{[}\widetilde{\mathcal{A}}\,\mathrm{d}\rho^{2}+\rho^{2}\Bigg{(}\frac{4\widetilde{\mathcal{B}}}{2-\xi^{2}}\left(\mathrm{d}\xi-\xi\,(2-\xi^{2})(1-\xi^{2})\,\rho\,\widetilde{\mathcal{F}}\,\mathrm{d}\rho\right)^{2}$
$\displaystyle\hskip
136.5733pt+\frac{(1-\xi^{2})^{2}}{h}\,s\,\widetilde{\mathcal{S}}\big{(}\mathrm{d}\phi-\,g_{-}^{2}w\,\widetilde{\mathcal{W}}\,\mathrm{d}t\big{)}^{2}\Bigg{)}\Bigg{]}\Bigg{\\}}\;.$
(19)

where
$\\{\mathcal{Q}\\}\equiv\\{\mathcal{T},\mathcal{A},\mathcal{B},\mathcal{F},\mathcal{S},\mathcal{W}\\}$
are unknown functions of $\\{x,y\\}$ and
$\\{\widetilde{\mathcal{Q}}\\}\equiv\\{\widetilde{\mathcal{T}},\widetilde{\mathcal{A}},\widetilde{\mathcal{B}},\widetilde{\mathcal{F}},\widetilde{\mathcal{S}},\widetilde{\mathcal{W}}\\}$
are unknown functions of $\\{\rho,\xi\\}$. The known functions, already
present in (A) and defined in (8) or (A), should be treated as scalars,
transforming between coordinate systems as (13). Note that when we set
$\\{\mathcal{T}=1,\mathcal{A}=1,\mathcal{B}=1,\mathcal{F}=0,\mathcal{S}=1,\mathcal{W}=1,\mathcal{W}=\beta\\}$
(and similarly for the associated set of tilde functions
$\\{\widetilde{\mathcal{Q}}\\}$), we recover the DeTurck reference metric (A).
The choice of our treatment of the factor $\beta$ allows us to more easily
recover the static results [33] in the limit $\beta\to 0$.

Most of the boundary conditions are determined by regularity, and we refer
readers to Section V of the review [47] for a detailed discussion of such
boundary conditions. Many of these are set by including factors of $x$, $1-x$,
$y$ and $1-y$ (and similar factors of $\xi$, $1-\xi$, $\rho$ and $1-\rho$ into
the metric ansatz. The $\mathbb{Z}_{2}$ reflection plane at $x=0$ ($\xi=0$)
deserves a few further comments. Our reference metric is even about this
surface, which agrees with the desired solutions with aligned spins, for which
we demand Neumann boundary conditions. For the anti-aligned case, we take the
same boundary conditions except for $\mathcal{W}$ ($\widetilde{\mathcal{W}}$),
for which we choose a Dirichlet condition. The fact that the reference metric
does not itself satisfy this condition does not lead to any contradictions.
This is indeed not a problem since this is happening at an axis of symmetry
(not at a boundary of the spacetime) and the reference metric still has well
defined parity (moreover, one explicitly checks that the de Turck norm indeed
vanishes along this axis even with the Neumann boundary conditions).

We supplement the boundary conditions by a set of patching conditions to join
the two different coordinate systems together. We defer a discussion of
patching to the next section.

For fixed cosmological constant, our solutions are physically parametrised by
$T_{+}$ and $\Omega_{+}$, but our solutions are parametrised by four
parameters: $\alpha$, $\beta$, $\lambda$, $k$. $\alpha$ is set according to
(15) to remove conical singularities, $\beta$ is uniquely determined by a
choice of $\Omega_{+}$, but any combination of our parameters $\lambda$ and
$k$ that give the same black hole temperature $T_{+}$ are physically
equivalent. After some trial and error, we find good numerical results by
fixing $\lambda=1/10$ and using $k$ and $\beta$ to parametrize the solutions.

Our (aligned or anti-aligned) spinning de Sitter binaries are thus a
2-parameter family of solutions parametrized by $T_{+}/T_{c}$ and
$\Omega_{+}$. Any combination of $k$ and $\lambda$ that give the same black
hole temperature $T_{+}$ are physically equivalent. To collect our numerical
data, we typically have fixed $\lambda=1/10$ and used $k$ and $\beta$ to
parametrize our solutions. We have tried different values of $\lambda$, but
after trial and error, this value typically generated the best numerical
results.

##  Appendix B Patching and Numerical Methods

In this section, we explain how we partition the domain of integration using
patching techniques (reviewed _e.g._ in [47]). The solution we seek contains
five boundaries: the inner segment of the axis ($\partial_{\phi}^{\rm in}$),
the black hole horizon ($\mathcal{H}^{+}$), the outer segment of the axis
($\partial_{\phi}^{\rm out}$), the cosmological horizon ($\mathcal{H}^{c}$),
and the plane of $\mathbb{Z}_{2}$ symmetry. Near the black hole event horizon
we use $(x,y)$ coordinates, while near the cosmological horizon we use
$(\rho,\xi)$ coordinates, and these coordinates at patched together.

For the results presented in this Letter, we used a total of five patches $-$
$I$, $II$,$III$, $IV$ and $V$ $-$ with each of these having four boundaries:
see Fig. 5. Patches $I,II,III$ and $IV$ are defined in $(x,y)$ coordinates,
and patch $V$ in $(\rho,\xi)$ coordinates. The patching boundary (dashed line
in Fig. 5) between patches $I$ and $II$ is given by constant $x=\chi_{1}$, and
between patches $II$ and $III$ it is given by constant $x=\chi_{0}$. The
patching boundary between patches $III$ and $IV$ is given by
$x=x_{0}y\sqrt{2-y^{2}}$. Finally, the patching boundary between patch $IV$
and patch $V$ is given by $\rho=\rho_{0}$.

Figure 5: Illustration of the patches used in our numerical construction. This
is for $\lambda=0.1$, $k=0.5$ and $\rho_{0},x_{0},\chi_{0},\chi_{1}$ as in
(20).

We fix the grid parameters $\rho_{0},x_{0},\chi_{0}$, and $\chi_{1}$ by:

$\rho_{0}=\frac{1}{10}\left(\frac{2}{\lambda}-\frac{1}{k}\right)+\frac{1}{k}\,\quad
x_{0}=\frac{7}{10}\left(1-\sqrt{1-\sqrt{1-\frac{1}{k^{2}\rho_{0}^{2}}}}\right)+\sqrt{1-\sqrt{1-\frac{1}{k^{2}\rho_{0}^{2}}}},\quad\chi_{0}=0.6\quad\text{and}\quad\chi_{1}=0.9\,,$
(20)

We apply the numerical methods detailed in [47], and discretize each of our
patches on a $N\times N$ Chebyshev-Gauss-Lobatto grid using transfinite
interpolation and pseudospectral collocation, for a total grid size of
$(N+N+N+N+N)\times N$. For $N=60$, the associated grid points are explicitly
displayed in Fig. 5, using a different colour code for each patch. After
discretization, the Newton-Raphson equation and boundary conditions are solved
with LU decomposition, using the $\beta=0$ (static) solution of [33] as a
first starting seed.

##  Appendix C Further results

Figure 6: Thermodynamic properties of Kerr-de Sitter.

Let us review the Kerr de Sitter black hole in order to make comparisons with
the spinning binaries. We begin by introducing the parameters
$y_{+}=\frac{r_{+}}{r_{c}}$ and $\tilde{a}=\frac{a}{r_{c}}$, where $r_{+}$ and
$r_{c}$ are the horizon radius and cosmological horizon radius, and $a$ is the
rotation parameter. Then in a frame with a non-rotating cosmological horizon,
the temperature, entropy, angular velocity, and angular momenta of the event
($+$) and cosmological ($c$) horizons of the Kerr-dS (KdS) black hole are:

$\displaystyle T_{+}^{\hbox{\tiny
KdS}}\ell=\frac{(1-y_{+})\left(2y_{+}^{3}+\left(1-3\tilde{a}^{2}\right)y_{+}^{2}-2\tilde{a}^{2}y_{+}-\tilde{a}^{2}\left(1+\tilde{a}^{2}\right)\right)}{4\pi\sqrt{y_{+}}\sqrt{y_{+}-\tilde{a}^{2}}\left(y_{+}^{2}+\tilde{a}^{2}\right)\sqrt{y_{+}^{2}+y_{+}+1+\tilde{a}^{2}}}\,,\qquad\frac{S_{+}^{\hbox{\tiny
KdS}}}{\ell^{2}}=\frac{\pi\left(y_{+}-\tilde{a}^{2}\right)\left(y_{+}^{2}+\tilde{a}^{2}\right)}{y_{+}^{3}+y_{+}^{2}+y_{+}+2\tilde{a}^{2}y_{+}-\tilde{a}^{4}}\,,$
$\displaystyle\Omega_{+}^{\hbox{\tiny
KdS}}\ell=\frac{\tilde{a}\left(1-y_{+}^{2}\right)\left(y_{+}^{3}+y_{+}^{2}+y_{+}+2\tilde{a}^{2}y_{+}-\tilde{a}^{4}\right)}{\left(1+\tilde{a}^{2}\right)\sqrt{y_{+}}\sqrt{y_{+}-\tilde{a}^{2}}\left(y_{+}^{2}+\tilde{a}^{2}\right)\sqrt{y_{+}^{2}+y_{+}+1+\tilde{a}^{2}}}\,,\quad\frac{J_{+}^{\hbox{\tiny
KdS}}}{\ell^{2}}=\frac{\tilde{a}\left(1+\tilde{a}^{2}\right)(y_{+}+1)\left(y_{+}-\tilde{a}^{2}\right)\left(y_{+}^{2}+\tilde{a}^{2}\right)}{2\left(y_{+}^{3}+y_{+}^{2}+y_{+}+2\tilde{a}^{2}y_{+}-\tilde{a}^{4}\right)^{2}}\,;$
$\displaystyle T_{c}^{\hbox{\tiny
KdS}}\ell=\frac{(1-y_{+})\left(\left(1-\tilde{a}^{2}\right)y_{+}^{2}+2\left(1-\tilde{a}^{2}\right)y_{+}-\tilde{a}^{2}\left(\tilde{a}^{2}+3\right)\right)}{4\pi\left(1+\tilde{a}^{2}\right)\sqrt{y_{+}}\sqrt{y_{+}-\tilde{a}^{2}}\sqrt{y_{+}^{2}+y_{+}+1+\tilde{a}^{2}}}\,,\qquad\frac{S_{c}^{\hbox{\tiny
KdS}}}{\ell^{2}}=\frac{\pi\left(1+\tilde{a}^{2}\right)\left(y_{+}-\tilde{a}^{2}\right)}{y_{+}^{3}+y_{+}^{2}+y_{+}+2\tilde{a}^{2}y_{+}-\tilde{a}^{4}}\,,$
$\displaystyle\Omega_{c}^{\hbox{\tiny
KdS}}\ell=0\,,\qquad\frac{J_{c}^{\hbox{\tiny
KdS}}}{\ell^{2}}=\frac{J_{+}^{\hbox{\tiny KdS}}}{\ell^{2}}\,.$ (21)

Regular Kerr-dS black holes exist for $0<y_{+}\leq 1$ and
$0\leq\tilde{a}\leq\tilde{a}_{\rm ext}(y_{+})$ where $\tilde{a}_{\rm ext}$ is
the value of $\tilde{a}$ at extremality where $T_{+}^{\hbox{\tiny KdS}}=0$. In
the left panel of Fig. 6, we display the microcanonical phase diagram with the
Kerr-dS black hole family. This plot displays the full Kerr-dS family of
solutions. The black line ($J_{+}^{\hbox{\tiny KdS}}=0$) describes the
Schwarzschild-dS solution, the red line describes the extremal Kerr-dS family
($T_{+}^{\hbox{\tiny KdS}}=0$) and the grey curve describes Kerr-dS black
holes in the Nariai limit $y_{+}=1$ (_i.e._ $r_{+}=r_{c}$).

To produce the microcanonical phase diagram displayed in Fig. 3 of the main
text, we had to find the corresponding Kerr-dS (or Schwarzschild-dS) solution
with the same $S_{c}/\ell^{2}$ and same angular momentum $J/\ell^{2}$ as the
spinning binary. For that, we used (C) to find the matching Kerr-dS parameters
$(y_{+},\tilde{a})$, and then used these to find the event horizon entropy of
this Kerr-dS black hole. This generates the magenta curve in Fig. 3 of the
main text or the cyan and magenta curves in the left panel of Fig. 3 (note
that for the cyan curve one has $a=J_{+}^{\hbox{\tiny KdS}}=0$, _i.e._ it
describes Schwarzschild-dS black holes) that are also displayed in the left
panel of Fig. 6 (the cyan curve on top of the black curve is barely seen since
it has very small length). As it is clear from Fig. 3 in the main text, the
2-parameter families of regular (anti-)aligned spinning binaries are described
by surfaces that are well below the Kerr-dS surface of Fig. 6 (for values of
$J/\ell^{2}$ and $S_{c}/\ell^{2}$ where they co-exist).

In the main text, we pointed out that for Kerr-dS black holes with fixed
temperature, the Komar angular momentum $J_{+}^{\hbox{\tiny KdS}}/\ell^{2}$
first increases with $\Omega_{+}^{\hbox{\tiny KdS}}\ell$ till it reaches a
maximum, and then decreases as $\Omega_{+}\ell$ increases further. This is
explicitly shown in the right panel of Fig. 6 where we plot
$J_{+}^{\hbox{\tiny KdS}}/\ell^{2}$ as a function of $T_{+}^{\hbox{\tiny
KdS}}\ell$ and $\Omega_{+}^{\hbox{\tiny KdS}}\ell$ and display three families
with constant $T_{+}^{\hbox{\tiny KdS}}\ell=0.05,0.1,0.2$ (yellow, green,
brown curves); note that Ker-dS exists for higher values of temperature and
angular velocity than those shown.

Figure 7: Ricci curvature evaluated at the black hole horizon (left panel) and
at the cosmological horizon (right panel)

Figure 8: Left panel: Gauge invariant metric component $-g_{tt}$, evaluated
along the $z-$axis containing the inner, black hole horizon and outer axis.
Right panel: $-g_{tt}$ as a function of the proper distance
$\mathcal{P}_{\phi}^{\rm out}$ between the event horizon and a point outside
it.

We have designed our numerical code such that the binary solutions are
parametrized by the black hole temperature and angular velocity. In
particular, some of the solutions we present have constant temperature
$T_{+}\ell=\frac{3.75}{2\pi}$ and $|\Omega_{+}^{(i)}\ell|=\beta$ with
$\beta\in[0,1.72]$. We stopped at $\beta=1.72$ because it becomes increasingly
more difficult to find numerical solutions with higher $\beta$. This occurs
because some of the functions start developing very large gradients and/or a
hierarchy of scales between different regions of the integration domain.
However, do not find no evidence for the appearance of curvature singularities
beyond $\beta>1.72$. One such curvature invariant that we monitored is the
pullback of the Ricci curvature into the black hole ($x=1$) and cosmological
horizons ($\rho=2/\lambda$), which are shown in the left and right panels of
Fig. 7, for three different values of $\beta$. These quantities remain clearly
finite as $\beta$ grows which suggests that our binaries should exist for even
higher values of $\beta$ without developing curvature singularities.

An important gauge invariant quantity of the binary system is the time-time
component of the metric $g_{tt}$. In the left panel of Fig. 8 we plot
$-g_{tt}$ as a function of the Weyl coordinate $z$ along the axis that
connects the two black holes (_i.e._ with Weyl coordinate $r=0$) and that
contains the inner $\partial_{\phi}$-axis ($z\in[0,1]$), black hole horizon
($z\in[1,1/k]$) and outer $\partial_{\phi}$-axis ($z\in[1/k,2/\lambda]$) for
anti-aligned binaries with $k=0.5$ for three different values of
$\Omega^{(1)}_{+}=\Omega^{(2)}_{+}=\beta$. We first notice the existence of an
ergoregion with $-g_{tt}<0$. As expected, we see that this ergoregion is
absent for $\beta=0$, and the minimum of $-g_{tt}$ becomes more negative as
$\beta$ increases. One also finds that $-g_{tt}$ becomes considerably more
spiked around its local maximum outside the black hole horizon (_i.e._ just
above $z=1/k=2$) and then decays faster with $z$ as $\beta$ increases. This
partially justifies why it becomes increasingly more difficult to find our
numerical solutions as $\beta$ grows. But it is important to note that this
does not imply the development of any pathology. Indeed, if instead of the
(gauge dependent) coordinate $z$, we plot $-g_{tt}$ as a function of the
(gauge invariant) proper length $\mathcal{P}_{\phi}^{\rm
out}=\int_{1/k}^{z}\sqrt{g_{\phi\phi}|_{r=0}}\mathrm{d}\tilde{z}$ along the
outer $\partial_{\phi}$-axis, as we do in the right panel of Fig. 8, one finds
that the qualitative behaviour of the curve (in particular, its width) does
not change significantly as $\beta$ increases. Anti-aligned spinning binaries
with $\Omega^{(1)}_{+}=-\Omega^{(2)}_{+}=\beta$ have similar qualitative
behaviour.

To have a more holistic illustration of how our solutions look like, (besides
Fig. 4 in the main text) in Figs. 9$-$10 we display the two metric functions
$-g_{tt}$ and $-g_{t\phi}$ (that are gauge invariant since $\partial_{t}$ and
$\partial_{\phi}$ are Killing vector fields) as a function of the original
cylindrical Weyl coordinates $(r,z)$ for two representative aligned spinning
binaries with black hole temperature $T_{+}/T_{c}=3.75$ and angular velocity
$\Omega_{+}^{(1)}\ell=\Omega_{+}^{(2)}\ell=1.0$ (Fig. 9) and
$\Omega_{+}^{(1)}\ell=\Omega_{+}^{(2)}\ell=1.7$ (Fig. 10). The system is
$\mathbb{Z}_{2}$ symmetric about $z=0$, _i.e._ $g_{ab}(r,-z)=\pm g_{ab}(r,z)$
(the minus sign holds only for the component $ab=t\phi$ in the anti-aligned
binary) and thus we just display the solution for $z\geq 0$. The five
different colours represent the five distinct patches that we use and the
colour code is the same one displayed in Fig. 5. Note that, as required, one
always has a smooth transition between patches. Anti-aligned spinning binaries
with $\Omega^{(1)}_{+}=-\Omega^{(2)}_{+}=\beta$ have similar qualitative
behaviour.

Figure 9: The gauge invariant metric functions $-g_{tt}$ and $g_{\phi\phi}$
for $T_{+}/T_{c}=3.75$ ($\lambda=0.1,k=0.5$) and
$\Omega_{+}^{(1)}\ell=\Omega_{+}^{(2)}\ell\equiv\beta=1.0$.

Figure 10: The gauge invariant metric functions $-g_{tt}$ and $-g_{t\phi}$ for
$T_{+}/T_{c}=3.75$ ($\lambda=0.1,k=0.5$) and
$\Omega_{+}^{(1)}\ell=\Omega_{+}^{(2)}\ell\equiv\beta=1.7$.

Recall that in this cylindrical Weyl coordinate chart $(r,z)$, we have a “rod
structure” where the rotation axis and the black hole horizons are all located
at $r=0$. More concretely, at $r=0$, the black horizon lies in the region
$z\in(1,1/k)$ and thus, the norm of the horizon generator
$K=\partial_{t}+\Omega_{+}^{(i)}\partial_{\phi}$ vanishes here (not plotted).
The ergoregion (where $-g_{tt}<0$) is also along this segment $z\in(1,1/k)$ at
$r=0$ but extends for $r>0$ in a neighbourhood of this rod. Figs. 9$-$10
indeed show that $-g_{tt}$ is negative in these ergoregions and that the
negative absolute minimum of $-g_{tt}$ (that is at $r=0$ as better pinpointed
by the magenta inverted triangles in Fig. 4) increases as $\beta$ increases.
Moreover, still at $r=0$, the inner segment of the axis between the black
holes is in the region $z\in(-1,1)$ and the outer segments of the axis are in
$z\in(1/k,\infty)$ and $z\in(-\infty,-1/k)$: not shown (see [33]),
$g_{\phi\phi}$ vanishes in these segments. The solution has a cosmological
horizon at $\sqrt{r^{2}+z^{2}}=2/\lambda$ (with temperature $T_{c}=1/(2\pi)$)
where $g_{tt}$ vanishes (since this horizon is generated by
$\partial_{t}+\Omega_{c}\partial_{\phi}$ with $\Omega_{c}=0$) as clearly
identified in the left panels of Figs. 9$-$10 .

In the right panels of Figs. 9$-$10 we plot $-g_{t\phi}$. This function
vanishes everywhere when $\beta=0$. When the components of the binary are
spinning, $-g_{t\phi}$ vanishes at $r=0$ along the inner $z\in(0,1)$ and outer
$z\in(1/k,\infty)$ segments of the axis. It also vanishes at the cosmological
horizon $\sqrt{r^{2}+z^{2}}=2/\lambda$ since $\Omega_{c}=0$, but is otherwise
non-zero. In particular, it is non-vanishing at $r=0$ along $z\in(1,1/k)$
inside of which it attains its maximum value, and this maximum value increases
as $\beta$ increases. Note that $-g_{t\phi}$ is an even function of $z$ for
aligned spinning binaries (the case that is shown in Figs. 9$-$10), but is an
odd function of $z$ in the anti-aligned binary case (not shown).

##  Appendix D Convergence Tests

In this section, we show that the norm $\chi\equiv\xi^{a}\xi_{a}$ of the
DeTurck vector vanishes in the continuum limit, as expected for a solution of
the Einstein-DeTurck equation that is not a Ricci soliton (_i.e._ that is
instead a true solution to the Einstein equation). Additionally, we find
exponential convergence, which is consistent with the use of Chebyshev
pseudospectral collocation methods.

Figure 11: Convergence test showing both the exponential accuracy of our
pseudospectral numerical method and the fact that we are not converging to a
Ricci soliton.

Let $\chi^{(N)}$ be $\chi$ computed on a (5-patched) grid with
$(N+N+N+N+N)\times N$ spectral collocation points. For concreteness, we take
$\alpha$ as given in (15), $\lambda=1/10$, $k=1/2$ and $\beta=1.5$. In Fig. 11
we show $\lVert\chi^{(N)}\rVert_{\infty}$ as a function of $N$ in a
$\log$-plot. The solid black line shows the best $\chi^{2}$-fit to a straight
line in the $\log$-plot, and yields

$f(N)=-14.97022-0.07694\,N\,.$ (22)

The exponential trend is clear and confirms that the Einstein-DeTurck solution
is converging to a true solution of the Einstein equation (and not to a Ricci
soliton with finite $\chi$).

## References

  * [1] R. Bach and H. Weyl, “New solutions to Einstein’s equations of gravitation. B. Explicit determination of static, axially symmetric fields (By Rudolf Bach). With a supplement on the static two-body problem (By Hermann Weyl). [Republication of: R. Bach and H. Weyl, Mathematische Zeitschrift 13 (1922) 134] ,” Gen Relativ Gravit 44 (2012) 817.
  * [2] W. Israel and K. A. Khan, “Collinear particles and bondi dipoles in general relativity,” Il Nuovo Cimento (1955-1965) 33 no. 2, (Jul, 1964) 331–344.
  * [3] W. Israel, “Event horizons in static vacuum space-times,” Phys. Rev. 164 (1967) 1776–1779.
  * [4] D. C. Robinson, “A simple proof of the generalization of Israel’s theorem.,” General Relativity and Gravitation 8 no. 8, (Aug., 1977) 695–698.
  * [5] G. L. Bunting and A. K. M. Masood-Ul-Alam, “Nonexistence of multiple black holes in asymptotically Euclidean static vacuum space-time.,” General Relativity and Gravitation 19 no. 2, (Feb., 1987) 147–154.
  * [6] S. D. Majumdar, “A class of exact solutions of Einstein’s field equations,” Phys. Rev. 72 (1947) 390–398.
  * [7] A. Papapetrou, “A Static solution of the equations of the gravitational field for an arbitrary charge distribution,” Proc. Roy. Irish Acad.(Sect. A) A51 (1947) 191–204.
  * [8] M. Astorino, R. Emparan, and A. Viganò, “Bubbles of nothing in binary black holes and black rings, and viceversa,” JHEP 07 (2022) 007, arXiv:2204.09690 [hep-th].
  * [9] R. M. Wald, “Gravitational spin interaction,” Phys. Rev. D 6 (1972) 406–413.
  * [10] N. Bohr, “On the constitution of atoms and molecules, part i,” Philosophical Magazine 26 (1913) 1.
  * [11] N. Bohr, “On the constitution of atoms and molecules, part ii, systems containing only a single nucleus,” Philosophical Magazine 26 (1913) 476.
  * [12] N. Bohr, “On the constitution of atoms and molecules, part iii, systems containing several nuclei,” Philosophical Magazine 26 (1913) 857.
  * [13] A. Svidzinsky, M. Scully, and D. Herschbach, “Bohr’s molecular model, a century later,” Physics Today 67 no. 1, (01, 2014) 33–39.
  * [14] Z. Perjés, “Solutions of the coupled einstein-maxwell equations representing the fields of spinning sources,” Phys. Rev. Lett. 27 (Dec, 1971) 1668–1670.
  * [15] W. Israel and G. A. Wilson, “A class of stationary electromagnetic vacuum fields,” J. Math. Phys. 13 (1972) 865–871.
  * [16] J. B. Hartle and S. W. Hawking, “Solutions of the Einstein-Maxwell equations with many black holes,” Commun. Math. Phys. 26 (1972) 87–101.
  * [17] A. Tomimatsu and H. Sato, “New Exact Solution for the Gravitational Field of a Spinning Mass,” Phys. Rev. Lett. 29 (1972) 1344–1345.
  * [18] D. Kramer and G. Neugebauer, “The superposition of two kerr solutions,” Physics Letters A 75 no. 4, (1980) 259–261.
  * [19] G. Neugebauer, “Relarivistic gravitational fields of rotating bodies,” Physics Letters A 86 no. 2, (1981) 91–93.
  * [20] H. Stephani, D. Kramer, M. MacCallum, C. Hoenselaers, and E. Herlt, Exact Solutions of Einstein’s Field Equations.  Cambridge Monographs on Mathematical Physics. Cambridge University Press, 2 ed., 2003.
  * [21] V. S. Manko and E. Ruiz, “Exact solution of the double-kerr equilibrium problem,” Classical and Quantum Gravity 18 no. 2, (Jan, 2001) L11.
  * [22] C. A. R. Herdeiro and C. Rebelo, “On the interaction between two Kerr black holes,” JHEP 10 (2008) 017, arXiv:0808.3941 [gr-qc].
  * [23] V. S. Manko, E. D. Rodchenko, E. Ruiz, and B. I. Sadovnikov, “Exact solutions for a system of two counter-rotating black holes,” Phys. Rev. D 78 (2008) 124014, arXiv:0809.2422 [gr-qc].
  * [24] V. S. Manko and E. Ruiz, “Metric for two equal Kerr black holes,” Phys. Rev. D 96 no. 10, (2017) 104016, arXiv:1702.05802 [gr-qc].
  * [25] V. S. Manko and E. Ruiz, “Metric for two arbitrary Kerr sources,” Phys. Lett. B 794 (2019) 36–40, arXiv:1806.10408 [gr-qc].
  * [26] V. S. Manko and E. Ruiz, “Equatorially symmetric configurations of two Kerr-Newman black holes,” Phys. Rev. D 105 no. 2, (2022) 024036, arXiv:2011.04313 [gr-qc].
  * [27] C. A. R. Herdeiro and E. Radu, “Two Spinning Black Holes Balanced by Their Synchronized Scalar Hair,” Phys. Rev. Lett. 131 no. 12, (2023) 121401, arXiv:2305.15467 [gr-qc].
  * [28] Supernova Cosmology Project Collaboration, S. Perlmutter et al., “Discovery of a supernova explosion at half the age of the Universe and its cosmological implications,” Nature 391 (1998) 51–54, arXiv:astro-ph/9712212.
  * [29] Supernova Search Team Collaboration, A. G. Riess et al., “Observational evidence from supernovae for an accelerating universe and a cosmological constant,” Astron. J. 116 (1998) 1009–1038, arXiv:astro-ph/9805201.
  * [30] Supernova Cosmology Project Collaboration, S. Perlmutter et al., “Measurements of $\Omega$ and $\Lambda$ from 42 high redshift supernovae,” Astrophys. J. 517 (1999) 565–586, arXiv:astro-ph/9812133.
  * [31] A. Wright, “It’s a mystery,” Nature Physics 7 no. 11, (Nov, 2011) 833–833.
  * [32] Planck Collaboration, P. A. R. Ade et al., “Planck 2013 results. XVI. Cosmological parameters,” Astron. Astrophys. 571 (2014) A16, arXiv:1303.5076 [astro-ph.CO].
  * [33] O. J. C. Dias, G. W. Gibbons, J. E. Santos, and B. Way, “Static Black Binaries in de Sitter Space,” Phys. Rev. Lett. 131 no. 13, (2023) 131401, arXiv:2303.07361 [gr-qc].
  * [34] J. E. Chase, “Event horizons in static scalar-vacuum space-times,” Commun. Math. Phys. 19 (1970) 276.
  * [35] R. Penney, “Axially Symmetric Zero-Mass Meson Solutions of Einstein Equations,” Phys. Rev. 174 (1968) 1578–1579.
  * [36] J. D. Bekenstein, “Transcendence of the law of baryon-number conservation in black hole physics,” Phys. Rev. Lett. 28 (1972) 452–455.
  * [37] J. D. Bekenstein, “Nonexistence of baryon number for static black holes,” Phys. Rev. D5 (1972) 1239–1246.
  * [38] J. D. Bekenstein, “Nonexistence of baryon number for black holes. ii,” Phys. Rev. D5 (1972) 2403–2412.
  * [39] C. Teitelboim, “Nonmeasurability of the quantum numbers of a black hole,” Phys. Rev. D5 (1972) 2941–2954.
  * [40] J. B. Hartle, Magic Without Magic.  edited by J. Klauder, Freeman, San Francisco, 1972.
  * [41] M. Heusler, “A No hair theorem for selfgravitating nonlinear sigma models,” J. Math. Phys. 33 (1992) 3497–3502.
  * [42] J. D. Bekenstein, “Novel ’no scalar hair’ theorem for black holes,” Phys. Rev. D51 (1995) 6608–6611.
  * [43] J. D. Bekenstein, “Black hole hair: 25 - years after,” in 2nd International Sakharov Conference on Physics, pp. 216–219.  5, 1996.  arXiv:gr-qc/9605059.
  * [44] D. Sudarsky, “A Simple proof of a no hair theorem in Einstein Higgs theory,,” Class. Quant. Grav. 12 (1995) 579–584.
  * [45] M. Headrick, S. Kitchen, and T. Wiseman, “A New approach to static numerical relativity, and its application to Kaluza-Klein black holes,” Class. Quant. Grav. 27 (2010) 035002, arXiv:0905.1822 [gr-qc].
  * [46] T. Wiseman, Numerical construction of static and stationary black holes.  in Black Holes in Higher Dimensions, CUP 2012, Ed. Gary T. Horowitz, 2011\.  arXiv:1107.5513 [gr-qc].  http://inspirehep.net/record/920553/files/arXiv:1107.5513.pdf.
  * [47] O. J. C. Dias, J. E. Santos, and B. Way, “Numerical Methods for Finding Stationary Gravitational Solutions,” Class. Quant. Grav. 33 no. 13, (2016) 133001, arXiv:1510.02804 [hep-th].
  * [48] R. Emparan and E. Teo, “Macroscopic and microscopic description of black diholes,” Nucl. Phys. B 610 (2001) 190–214, arXiv:hep-th/0104206.
  * [49] P. Figueras, J. Lucietti, and T. Wiseman, “Ricci solitons, Ricci flow, and strongly coupled CFT in the Schwarzschild Unruh or Boulware vacua,” Class. Quant. Grav. 28 (2011) 215018, arXiv:1104.4489 [hep-th].
  * [50] P. Figueras and T. Wiseman, “On the existence of stationary Ricci solitons,” Class. Quant. Grav. 34 no. 14, (2017) 145007, arXiv:1610.06178 [gr-qc].
  * [51] In a forthcoming publication [DSW-chargedBinary], we construct charged de Sitter binaries, and we see a maximum charge as the electric potential increases.
  * [52] G. W. Gibbons and S. W. Hawking, “Cosmological Event Horizons, Thermodynamics, and Particle Creation,” Phys. Rev. D 15 (1977) 2738–2751.
  * [53] P. G. LeFloch and L. Rozoy, “Uniqueness of kottler spacetime and the besse conjecture,” Comptes Rendus Mathematique 348 no. 19, (2010) 1129–1132.
  * [54] S. Borghini, P. T. Chruściel, and L. Mazzieri, “On the uniqueness of Schwarzschild-de Sitter spacetime,” arXiv:1909.05941 [math.DG].
  * [55] D. Katona and J. Lucietti, “Uniqueness of the extremal Schwarzschild de Sitter spacetime,” Lett. Math. Phys. 114 no. 1, (2024) 18, arXiv:2309.04238 [gr-qc].
  * [56] D. Katona, “Uniqueness of extremal charged black holes in de Sitter,” arXiv:2403.08467 [gr-qc].
  * [57] R. Ruffini and J. A. Wheeler, “Introducing the black hole,” Phys. Today 24 no. 1, (1971) 30.
  * [58] M. Mathisson, “Republication of: New mechanics of material systems, acta phys. pol. 6 (1937) 163,” General relativity and gravitation 42 no. 4, (2010) 1011–1048.
  * [59] A. Papapetrou, “Spinning test-particles in general relativity. i,” Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 209 no. 1097, (1951) 248–258.
  * [60] E. Corinaldesi and A. Papapetrou, “Spinning test-particles in general relativity. ii,” Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 209 no. 1097, (1951) 259–268.
  * [61] W. G. Dixon, “Dynamics of extended bodies in general relativity. I. Momentum and angular momentum,” Proc. Roy. Soc. Lond. A 314 (1970) 499–527.
  * [62] W. G. Dixon, “The definition of multipole moments for extended bodies,” General Relativity and Gravitation 4 (1973) 199–209.
  * [63] W. G. Dixon, “Dynamics of extended bodies in general relativity III. Equations of motion,” Phil. Trans. Roy. Soc. Lond. A 277 no. 1264, (1974) 59–119.