# A NICER viewing angle on the accretion stream and disk of Vela X-1

Roi Rahin  NASA Goddard Space Flight Center, 8800 Greenbelt Road, Greenbelt,
MD 20771, USA  Physics Department, Technion, Haifa 32000, Israel  Ehud Behar
Physics Department, Technion, Haifa 32000, Israel

###### Abstract

Vela X-1 is the archetypical eclipsing high-mass X-ray binary, composed of a
neutron star (NS) accreting the B-star wind. It was observed by nearly all
X-ray observatories, often multiple times, featuring a rich spectrum of
variable emission lines. Yet, the precise origin of these lines in the binary
system remains uncertain. We perform a systematic, orbital phase-dependent
analysis of the reflected Fe K$\alpha$ fluorescence line at 6.4 keV using over
100 NICER observations. We resolve the line variability into 500 s time bins
and find that it is predominantly due to variation in the ionizing flux, with
a moderate underlying phase dependence over the 9 day orbital period. Our
analysis reveals a significant reflection component that cannot originate from
the companion B-star alone. We also find that an appreciable portion of the
B-star surface is obscured opposite the eclipse, and this obscuration is not
symmetric around the mid-point (phase=0.5). We argue that an accretion stream,
from the B-star to the NS and distorted by the orbital motion, is responsible
for both the additional fluorescence emission component and for obscuring the
B-star.

††facilities: NICER (Gendreau et al., 2016), Chandra/HETG (Canizares et al.,
2005), XMM-Newton/RGS (Den Herder et al., 2001) ††software: Xspec (Arnaud,
1996), pyXspec, HEAsoft (HEASARC, 2014)

##  1 Introduction

Vela X-1 is a high mass, eclipsing, X-ray Binary (XRB). Vela X-1 consists of a
neutron star (NS) of $\sim 2M_{\odot}$ and a massive B-star ($M_{*}\sim
25M_{\odot}$,$R_{*}\sim 30R_{\odot}$) with a short orbital period of $\sim 9$
days (Kreykenbohm et al., 2008). The system is highly compact, with a
separation of only $\sim 1.8R_{*}$, which implies the NS is deeply embedded in
the dense stellar wind. Vela X-1 is a consistent X-ray source of moderate
luminosity ($\sim 10^{36}$ erg s-1), but at a short distance of $\sim 2$ kpc
(Kretschmar et al., 2021), which makes it one of the brightest X-ray point
sources in the sky. The X-ray spectrum of Vela X-1 is commonly modeled with an
absorbed cut-off (above 20 keV) power-law continuum, photo-ionized emission
line components, and a reflection component with both continuum and (neutral)
fluorescence lines.

Vela X-1 is highly variable in X-ray luminosity, line of sight column density,
and emission line spectrum. The luminosity variations are attributed to giant
flares, off-states, and column density variations (Kreykenbohm et al., 2008;
Kretschmar et al., 2021, and references therein). The presence of an accretion
stream from the B-star enshrouding the NS is evidenced by partial absorption
of the 2-20 keV X-ray continuum, but not the harder X-rays between 15-50 keV,
even when the NS is in front of the B-star at phase 0.5 (see Figure 4 in
Kretschmar et al. (2021)). The variability in column density is associated
with the accretion stream and clumps in the stellar wind (Watanabe et al.,
2006; Grinberg et al., 2017). Simulations of high-mass XRBs show this column
density is expected to be asymmetric around phase 0.5 (Blondin et al., 1991;
Manousakis et al., 2012), with observations in good agreement (Kretschmar et
al., 2021). The emission line spectrum is variable as well. Grinberg et al.
(2017) noted extreme variability in ionized line emission around phase 0.25.
The reflection component and its fluorescence lines vary as well and are
assumed to originate in a cold and dense neutral medium. Previous works
ascribed most of the fluorescence to the surface of the stellar companion.
High observed equivalent width (EW) of the Fe K$\alpha$ line at phases 0.5 and
0.25 implied an additional reflection component (Watanabe et al., 2006).

The cause of the line variability is a standing question in Vela X-1 and may
contribute to understanding the geometry of the binary and its accretion
stream. Furthermore, the origin and physical mechanism of both ionized
emission and fluorescence in Vela X-1 is still debated (Mauche et al., 2008;
Kretschmar et al., 2021). A systematic analysis of multiple phases and
accretion states may improve our understanding of this enigmatic source. The
challenge is to disentangle the slowly varying phase dependence from that of
the stochastic accretion process. In this paper, we harness the high
sensitivity and high cadence observations of Vela X-1 by NICER (NS Interior
Composition Explorer) to do just that.

##  2 Vela X-1 observations and spectral models

###  2.1 NICER observations and modeling

NICER (Gendreau et al., 2016) is a soft X-ray (0.2 – 12 keV) telescope mounted
on the ISS (International Space Station). NICER observes Vela X-1 regularly as
a calibration source, but many of these observations are valid for science
analysis. We use all observations of Vela X-1 up to 2022-05-04 and divide them
into segments of up to 500 continuous seconds using
Ftools111http://heasarc.gsfc.nasa.gov/ftools. In practice, some segments
contain less than 500 seconds of data after applying NICER’s Good-Time-
Intervals procedure. We use a
pyXspec222https://heasarc.gsfc.nasa.gov/xanadu/xspec/python/html/index.html
script to analyze all resulting observation segments that have $>10000$ total
photons in the 3-9 keV band. This excludes data around the eclipse, from
phases 0.9 - 0.1. We use a model combining an absorbed power-law, a gaussian
emission line for the 6.4 keV Fe-line, and unabsorbed continuum reflection
(Magdziarz & Zdziarski, 1995). The final Xspec model used by the fitting
script is:

$TBabs*(TBabs*powerlaw+pexrav+gauss)$ (1)

The first TBabs absorption component is the Galactic absorption derived from
HI maps 333https://heasarc.gsfc.nasa.gov/cgi-bin/Tools/w3nh/w3nh.pl and its
column density is set to $0.37\times 10^{22}$ cm-2. The second TBabs
absorption component is intrinsic to Vela X-1 and is fitted for each NICER
segment. The model assumes the absorbing medium is a neutral cold gas with
solar composition. We originally attempted to use TBfeo to fit the Fe
abundance, but the Fe abundance can not be well constrained using this data,
and we thus returned to solar abundances. The power-law slope and norm are
fitted, and both are tied to the pexrav power-law slope and norm. In the
pexrav model, the abundance is set to the solar value, while the angle cosine,
defined between the line of sight and reflector inclination, is set to 0.95
since the inclination of Vela X-1 is known to be close to $90^{\circ}$, and
the dominant reflection is expected to be the B-star surface. The power-law
cutoff energy is set to 20 keV (Makishima et al., 1999). Since the X-ray
source is a compact object, it is easily obscured by wind clumps and accreting
matter. By comparison, the reflector is at least the size of the B-star, with
a radius of $\sim 30R_{\odot}$, and will only be obscured by large-scale
geometrical elements. Therefore, we remove the pexrav direct power-law
component (by assigning the rel_refl component a negative value), to allow the
power-law component to be absorbed separately while keeping the reflection
component unabsorbed. In summary, the free model parameters are the power-law
slope and norm, the column density of the local absorber, the gaussian
emission line parameters, and the reflection scaling.

###  2.2 Chandra and XMM-Newton observations

We utilized Chandra/HETG observations, processed by The Chandra Grating-Data
Archive and Catalog, TGCat (Huenemoerder et al., 2011), and XMM-Newton/RGS
observations to measure emission lines from Vela X-1. The eight Chandra/HETG
observations of Vela X-1 were analyzed multiple times (e.g., Schulz et al.,
2001; Goldstein et al., 2004; Watanabe et al., 2006; Grinberg et al., 2017;
Rahin & Behar, 2020). The three RGS observations were analyzed by Martínez-
Núñez et al. (2014). For a more comprehensive list of observations see
Kretschmar et al. (2021). In this work, We use the Chandra/HETG observations
1927,1926,19953,19952 and 14654 as well as the XMM-newton observation
0841890201 (see Table 1).

Throughout this work, We calculate the orbital phase of each measurement using
the reported eclipse epoch and orbit period from Kreykenbohm et al. (2008).
Note that phases calculated in this way may show a slight deviation from other
works.

##  3 sporadic vs. phase-dependent spectral variability

As a binary system, Vela X-1 is expected to exhibit periodic, phase-dependent
behavior. Such behavior can be seen in the averaged observed X-ray flux over
the 9 d period, although surprisingly the peak flux in the 2-20 keV range is
found around phase $\sim 0.3$ (Kretschmar et al., 2021). On the other hand,
Vela X-1 varies on much shorter time scales, likely due to the sporadic
accretion process. Spectral elements such as column density, line flux, and
power-law slope show substantial variability on the shortest time scales where
variability can be decidedly measured ($\sim$1 h). In this section, we try to
identify and separate these two forms of variability, phase-dependent and
sporadic, by focusing on the column density and the Fe K$\alpha$ line flux,
with a note about the photo-ionized gas.

###  3.1 Column density and complete obstruction

The column density absorbing the X-ray continuum of Vela X-1 is known for high
variability both with phase and between observations (Kretschmar et al.,
2021). Figure 1 shows the column density across 608 NICER observation segments
of Vela X-1. Each observation is at most 500 s long. Blue dots represent
individual measurements and red stars are averages over 0.05 phase bins. The
red error bars represent the standard deviation between different
measurements. When performing averages, we first average all observations that
are temporally adjacent (defined as being at most one hour apart) and then
take the average of all other observations in the same bin. Since Vela X-1 is
known to have periods of extreme absorption (Liao et al., 2020), we exclude
outliers with $N_{\textrm{H}}>50\times 10^{22}$ cm-2 from both the figure and
the averaging process. This excludes 22 measurements both for figure clarity
and to avoid skewing the average with extremely high outlier measurements.
Almost all excluded measurements lie between phases 0.3 and 0.6. Both the low
number of excluded measurements and their much higher column density lends
support to their classification as outliers. The X-ray source is more absorbed
on average at phases 0.5–0.9 than at 0.1–0.5 indicating a breaking of symmetry
in the system, where high column density gas is located (and rotates with the
system) on one side of the binary axis. This effect, suggesting an accretion
stream as the absorber, was previously noted by others (Kretschmar et al.,
2021, and references therein) and fits high-mass XRB simulations Blondin et
al. (1991); Manousakis et al. (2012), but it is shown here with significantly
more data points, owing to NICER’s sensitivity and numerous observations.

Figure 1:  Equivalent $N_{\textrm{H}}$ measured in various NICER observation
slices of up to 500 s. Blue dots are individual measurements and red stars are
averages over 0.05 phase bins. Note that the red error bars represent the
standard deviation between measurements, not the error of the mean. A
significant trend is visible where $N_{\textrm{H}}$ decreases to a minimum at
phase $\sim\\!0.3$ and then increases. This asymmetric behavior with phase
indicates excess gas on one side of the system, which is likely associated
with a (rotating) accretion stream.

###  3.2 The Fe K$\alpha$ line

One of the most conspicuous features of the Vela X-1 X-ray spectrum is the Fe
K$\alpha$ fluorescence line. This line was reported in many observations and
studied extensively Nagase et al. (1994); Goldstein et al. (2004); Watanabe et
al. (2006); Grinberg et al. (2017); Kretschmar et al. (2021), and yet there is
no clear understanding of the line source or its phase-dependence. We use the
500 s sliced NICER observations (Section 2.1) to obtain 630 measurements of
the Fe K$\alpha$ line. Although the actual variability time of the line may be
even shorter than 500 s, to the best of our knowledge this is the largest
sample with such short times segments to be presented for Vela X-1. We use it
to shed light on the physical origin of this line.

Figure 2:  The Fe K$\alpha$ line flux variations over 1 day as recorded with
several NICER observations. Each observation was broken into segments of at
most 500 s. There is a significant variation of the line flux, sometimes on a
timescale of one hour.

Figure 2 shows Fe K$\alpha$ flux measurements from the start of observation
2020200102 over a duration of 1 day. The measurements are broken into 500 s
segments and show the rapid variability of the line. The figure shows
significant difference up to a factor of 2 in line flux between observations
($\sim 5$ ks). Variability between adjacent 500 s segments seems to exist but
is not significant to the level that we can measure here. Figure 3 shows the
relation between the Fe K$\alpha$ flux and the model unabsorbed flux in the 7
keV to 20 keV range derived from the spectral model described in Section 2.1.
The Pearson correlation coefficient is $0.78$ with $p<0.001$. Since we expect
the line to be primarily driven by the Vela X-1 continuum, this relation is as
expected.

Figure 3:  The Fe K$\alpha$ line flux vs. the unabsorbed continuum flux in the
7-20 keV range. The correlation is readily apparent and expected.

The above time scales are much shorter than the orbital period. We wish to
also understand the effect of the orbital phase on the Fe K$\alpha$ line.
Figure 4 shows the Fe K$\alpha$ line flux across multiple measurements in
various phases. In an attempt to detect trends, We calculated the average flux
over bins of size 0.05. Unlike the correlation with continuum flux, no
discernible pattern emerges, and the mean line flux remains between
0.0015–0.004 photons cm-2 s-1.

Figure 4:  Fe K$\alpha$ line flux with phase. Red stars are average values
over bins of size 0.05. No clear dependence on phase emerges.

Since we aim to discern the phase dependency of the Fe K$\alpha$ line, we must
disentangle the strong dependence on the continuum. We thus divide the line
flux by the ionizing unabsorbed model flux in the $7-20$ keV range. Figure 5
shows this fraction of the iron line relative to the continuum with phase,
revealing a trend that was not visible in Figure 2. Once the continuum effect
is folded out, this trend can shed light on the geometry of the binary. In
Figure 5, we see a rise in relative flux with phase up to $\sim$0.35, which
interestingly is also where the $2-20$ keV continuum peaks Kretschmar et al.
(2021).

Figure 5:  The Fe K$\alpha$ flux relative to the continuum flux in the 7-20
keV range with phase. The red stars are average values over bins of size 0.05.
There is a clear increase in relative flux up to phase $\sim 0.35$ followed by
a slow decrease. The shape is clearly asymmetrical.

Furthermore, in and out of the eclipse the line flux relative to the continuum
shows an increase towards a virtual maximum at the center (Figure 5). This is
expected if the reflector is the half-sphere surface of the B-star, rotating
into and out of our line of sight, analogous to the lunar phases. However, two
strong deviations from this behavior are apparent. First, in the middle of the
period between phases $0.3-0.7$ the relative line flux is much lower than
expected from this trend, implying that part of the reflecting star is
obstructed from our line of sight. The second deviation is that the line flux
relative to the continuum does not tend to zero towards the eclipse (at phases
0.9 and 0.1), as would be expected from a half-lit rotating sphere. On the
contrary, there is significant Fe K$\alpha$ emission, of the order of $0.5\%$
relative to the continuum, throughout the orbit. Note that although some line
emission is present around the eclipse in phases $0.9-0.1$, there is no
continuum to divide by, and no data is available for Figure 5.

###  3.3 Simple geometrical reflection model

The Fe K$\alpha$ line is associated with the reflection of the X-rays from the
direction of the NS to our line of sight, which is highly dependent on
geometry. Since the primary suspect as a reflector is the B-star, we first
attempt a model of reflection from an optically thick sphere, similar to the
initial assumption of Watanabe et al. (2006). The present calculation
generalizes the 1D slab approximation described in Rahin & Behar (2020).
Unlike in a slab, the entry (continuum excitation) and exit (of the
fluorescence photon) paths are different. The exit path length to our line of
sight is calculated geometrically for all incident angles entering the sphere,
and further depends on the entry path length, and the orbital phase. The model
includes the quantum probability to excite a Fe K$\alpha$ line inside the
star, and its probability to escape the medium (see Rahin & Behar, 2020). It
ignores Compton scattering, whose effect is smaller than the measurement
uncertainties (George & Fabian, 1991). See Figure 6 for a 2D drawing of the
geometry. A detailed derivation of the reflected line intensity can be found
in Appendix A.

Figure 6:  The 2D model for reflection of a point source (the NS at the
leftmost point) from an optically thick sphere, into our line of sight (the
observer is at the bottom reflecting an 0.25 or 0.75 phase), drawn with Desmos
(https://www.desmos.com/geometry) . The stellar surface is denoted by the
green circle. The ionizing continuum enters and exits the sphere through a
path length of $x$ and $y$, respectively.

Figure 7 shows the Fe K$\alpha$ flux and the theoretical value from the
spherical reflection model. As can be seen, the model over-predicts the
reflected line flux around phase 0.5 and under-predicts it close to the
eclipse. Specifically, the average flux measured around the 0.5 phase is about
5 times lower than expected and the average flux measured near eclipse (phase
0.1-0.15) is about 5 times higher than expected. The over-prediction could be
a result of one or more factors. The star could be partially obscured, either
from us or from the source itself. In the latter case, the ionizing radiation
exciting the material on the face of the B-star would be lower than assumed by
the model, lowering the expected emission. The ionizing radiation could also
have a different continuum spectrum than that used in the model. For example,
if the power-law photon index of the radiation reflected off the star is
higher than the measured power-law index, the expected fluorescence will be
considerably lower. Compared to the model over-prediction, the under-
prediction is more difficult to explain within the confines of the model. The
model predicts little emission close to the eclipse, but the observed emission
is still considerable. Any attempt to reconcile this discrepancy within the
sphere-reflection model will necessarily only increase the model over-
prediction. It must therefore be deduced that the model describing reflection
from a sphere cannot describe the observed Fe K$\alpha$ flux or its orbital
dependence. This implies both the presence of additional reflecting material
apparent most clearly near the eclipse and the presence of material obscuring
a large portion of the surface of the star, most conspicuously between phases
$0.3-0.7$.

Figure 7:  Relative Fe K$\alpha$ line flux with phase (blue, and 0.05 phase
averaged in red), compared to the expected reflection from the surface of a
sphere (green) relative to the continuum. The model over-predicts the
measurements between phase $0.2-0.7$ and under-predicts it elsewhere. Clearly,
reflection from a sphere alone cannot explain the observed line flux.

###  3.4 Photoionized-gas emission lines

The Vela X-1 Chandra/HETG and XMM-Newton/RGS spectra are rich in emission
lines from a variety of ion species. Here, we rely on the He-like ions to
provide insights into the properties of the emitting gas, in an attempt to
complete the geometrical picture of Vela X-1. Specifically, the $R=f/i$ ratio
of the forbidden ($f$) to intercombination ($i$) lines can probe the electron
density and distance from the B-star (Porquet & Dubau, 2000; Leutenegger et
al., 2006). Both high electron number density and high B-star UV flux incident
on the X-ray gas will reduce the $R$ ratio from its nominal low-density value.
Figure 8 shows the $R$ ratio of He-like O, Ne, and Mg at several orbital
phases. Some phase dependency is suggested with $R$ having higher values
during and around the eclipse.

Figure 8:  The $R=f/i$ ratio of O, Ne, and Mg across several orbital phases of
Vela X-1. Phase 0.43 measurements are from the RGS, and include O, while all
other measurements are from HETG (see Appendix B). Measured $R$ values are
much lower than the nominal ones of 3.7, 3.1, and 2.7 for O, Ne, and Mg
respectively, implying a significant effect of high UV flux or electron
density.

Differentiating UV flux and electron density effects in photo-ionized plasma
can be challenging. They decrease together with the distance from the B-star
(both as $r^{-2}$), implying the effects of both on the R ratio increase and
decrease in tandem. This was previously noted by Peretz et al. (2019) in AGN
observations. We will now argue that in the case of Vela X-1, the dominant
effect is UV pumping. Since the donor in the Vela X-1 binary is a B-star, it
is natural to first consider the impact of the UV radiation field. Figure 9
shows the implied distances from the B-star in stellar radii $r/R_{*}$, given
the measured $R$ ratios, as a function of phase. We rely on the method
introduced by Leutenegger et al. (2006, Eqs. 1-4 and Table 1 therein), and
adopt the stellar parameters used by Grinberg et al. (2017) ($T_{eff}=25,000$
K, $\log g=3.0$), who performed a similar calculation for the Chandra
observation 1928. We find that the He-like lines are emitted from
$2-20\,R_{*}$ from the UV source, presumably the B-star. The distance of the
NS from the B-star is $\sim 1.8R_{*}$. This implies the emitting gas is
farther than the NS, sometimes considerably. For Mg, the difference is small
and may be a result of model assumptions, but for Ne and O, the distance of
the emitting gas is $r>5R_{*}$. The emitting Oxygen gas is particularly far at
a distance $r>12R_{*}$ from the B-star.

Figure 9:  Distance of He-like O, Ne, and Mg ions from the OB donor star
across several phases of Vela X-1. The distance is calculated from the UV
effect on the R-ratio at the low-density limit. Notice the consistently
shorter distance with atomic number, with Mg the closest to the B-star, and O
the farthest at $\sim 20R_{*}$.

In theory, lower than nominal $R$ values can also be caused by a high number-
density of the emitting gas. The densities for which the $R$ ratio of O, Ne,
and Mg deviates from its nominal value are greater than $10^{9}$ cm-3,
$10^{10}$ cm-3, and $10^{11}$ cm-3, respectively (Porquet & Dubau, 2000). We
can obtain a rough approximation for the stellar wind density profile by
assuming for simplicity a density profile $n(r)=n_{0}/(r/R_{*})^{2}$. The
total (radial) column density of such a wind is $N_{\textrm{H}}=n_{0}R_{*}$.
For $N_{\textrm{H}}=3-4\times 10^{22}$ cm-2 (Rahin & Behar, 2020) and
$R_{*}=2\times 10^{12}$ cm, $n_{0}\sim 10^{10}$ cm-3. This estimated density
is consistent with those obtained in simulations (Mauche et al., 2008;
Manousakis, 2011). Since the wind density decreases rapidly from this value
with $r$, the only element that might be affected by it is O. However, O was
already shown to originate far from the B-star at $r/R_{*}>10$, where
$n/n_{0}<100$ (Figs. 8, 9). The effect of gas density on the $R$ ratio in a
smooth wind is therefore negligible. It is still possible that the He-like
lines originate from dense wind clumps, but these would have to be even
farther away from the star than indicated in Fig. 9.

We did not include the $R$ ratio of He-like Si in this discussion, despite it
being well measured with values between $3.5-8$, since this is substantially
higher than the theoretical nominal limit of 2.3 (Porquet & Dubau, 2000). This
anomaly was noted by Grinberg et al. (2017) and by Amato et al. (2021), with
no explanation. Unfortunately, this work cannot offer a definitive explanation
either. Examinations of other photo-ionized X-ray sources, such as the
Circinus Galaxy with $R\textrm{(Si)}\sim 4$, hint that the high value is
prevalent. The peculiarly high R ratio can have several possible explanations.
First, there could be an instrumental artifact in the Chandra/HETG which is
not corrected for by the response functions. Given the close proximity of the
intercombination and forbidden lines, this is unlikely. Second, if a 1s
electron is ionized from the ground state of Li-Like Si, the resulting He-like
ion is in the 1s2s excited state. This will contribute to the emission in the
forbidden line and increase R. However, the fraction of Li-like ions in
steady-state is small. Third, and most likely, the theoretical values taken
from Porquet & Dubau (2000) may require updating. Future observations of Vela
X-1, and other X-ray sources with XRISM/Resolve (Tashiro et al., 2020) may aid
in determining the exact cause of these high R ratios for Si, which may call
for re-examining the theoretical values. Since the theoretical R-value is
expected to increase, the derived distances will only increase. The
conclusions described above will thus remain valid in this case.

The upshot of this analysis of the He-like line ratios is that the ionized gas
is located several stellar radii away from the B-star and is likely unrelated
to the reflecting medium discussed above.

##  4 Discussion

The column density observed throughout the orbit of Vela X-1 shows a clear
asymmetry around phase 0.5, with low $N_{\textrm{H}}$ in phases 0.1–0.5 and
high $N_{\textrm{H}}$ in phases 0.5–0.9. These results agree with other
measurements (Kretschmar et al., 2021, and references therein) and with
simulations of eclipsing high-mass XRBs (Blondin et al., 1991; Manousakis et
al., 2012). As explained in these references, this is a result of the Coriolis
force pushing away the accretion stream to one side of the system. Comparing
our results to the simulations of Manousakis et al. (2012) (Figures 2 and 3
therein), we note the lack of a column density ”bump” around phases 0.4–0.6.
Considering the estimated mass of the NS in Vela X-1 of $1.7-2.2\,M_{\odot}$,
the best matching parameters in their simulations would be $v_{\infty}\sim
1200$ km s-1 and a wind mass-loss rate of $\dot{M}_{W}\sim 10^{-6}M_{\odot}$
yr-1. These rough estimates are in agreement with previous observations
(Nagase et al., 1986; Prinja et al., 1990), but here derived from the phase
dependence of the column density.

The fluorescence spectrum of Vela X-1, considering its high intensity, must
originate as a reflection from a dense, cold region with a high column
density. A simple model for reflection from the spherical surface of the
B-star overestimates the observed flux between phases 0.3 and 0.7 by up to a
factor of 5 but underestimates the observed flux around the eclipse by a
similar factor (Figure 7). These drastic discrepancies are unlikely to be
fixed by small modifications to the model. We propose that an accretion stream
in Vela X-1 is responsible for both of these deviations. The fact that the
reflected line relative to the continuum is not symmetric with phase (Figure
5) indicates the stream is located on one side of the symmetry axis of the
binary. It obscures part of the B-star between phases 0.35–0.75 and adds
steady fluorescence emission, $\sim 0.5\%$ of the continuum, throughout the
orbit. This interpretation requires the accretion stream to be dense enough to
have a low ionization parameter and to emit Fe K$\alpha$. On the other hand,
it can not be Compton thick, since in phases 0.5–0.9 we observe the reflection
through the stream.

The majority of the reflection component was postulated to originate on the
surface of the companion star (Watanabe et al., 2006). However, with the
stellar wind, this surface may not be well defined and other components such
as wind clumps, and accretion streams may contribute to the reflection
components. Indeed, Watanabe et al. (2006) observed with Chandra excess EW at
two discrete phases 0.25 and 0.5 compared to the reflection contribution of
the B-star and its wind, in their model. They thus invoked a reflecting cloud
to explain this excess. The present analysis of the line flux relative to the
unabsorbed continuum above 7 keV is more physical since it is this continuum
that excites the line, rather than the local continuum at 6.4 keV that does
not. In fact, we find a deficiency of relative flux at phase 0.5 (Figure 5).
NICER’s full orbital coverage and shorter observation times presented in this
work enable a more complete picture of the reflection with phase, and thus
consistently point to the Vela X-1 accretion stream. Note that using the EW
for this kind of analysis could be misleading since variable continuum
absorption, for instance, will change the EW without any change in line flux.
It appears that the emission from the accretion stream is visible from both
sides of the binary (phases 0.25 and 0.75), as the deviation from a reflecting
sphere appears on both ends of the phase plot.

Further verifying this hypothesis can be done using high-resolution
spectroscopy of Vela X-1 immediately before and after the eclipse. The
fluorescence line ratios of different elements can then be compared to
estimate the fluorescing medium column density (Rahin & Behar, 2020).
Immediately after the eclipse, the sole reflecting medium that can be observed
is the accretion stream, and this measurement will yield its column density.
Similarly, immediately before the eclipse, the accretion stream is observed
from the other side and a measurement of the fluorescing column density viewed
through the stream should yield similar values. XRISM, expected to launch in
2023, is a suitable instrument for such an observation.

As shown in Section 3.4, the emission lines of the ionized gas in Vela X-1
provide information about its location. Our analysis of the He-like $R$ ratios
and their dependence on UV flux indicates that the minimum distance of the He-
like Ne and Mg ions from the B-star is greater than the distance to the NS.
This result was reported by Grinberg et al. (2017) for phase 0.25, but we show
it to be consistent across all phases. Hence, the ionized gas is more remote
than the putative accretion stream. In hydrodynamical simulations by Mauche et
al. (2008), He-like Si is found between the NS and B-star with the majority of
emission in a region close to the surface of the B-star $(r\approx 1R_{*})$.
Since He-like Mg forms in similar $\xi$ values as He-like Si, our findings of
$r=2R_{*}$ do not agree with these results. The emission of ionized Ne and Mg
during the eclipse (phase 0) appears to originate farther from the star than
between phases 0.1–0.9. This is consistent with the overall weakness of the
ionized lines in eclipse. Indeed, Watanabe et al. (2006) concluded that the
majority of ionized lines are produced in a region obscured by the host B-star
during the eclipse.

Interestingly, in our results the ions appear ordered in distance by
decreasing atomic numbers, with Mg originating closest to the star, Ne
following, and O originating some distance away, at almost $20R_{*}$ (Figure
9). This could imply a gradient in the ionization parameter of the stellar
wind. The emission of He-like O specifically, probably originates in the
sparse wind of the B-star illuminated by the X-ray source. Taken together,
these results imply a co-rotating ionization cone of gradually decreasing
ionization diverging away from the binary system (increasing $r$). This
implies the density decreases less steeply than $n\propto r^{-2}$, suggesting
collimation, possibly by magnetic fields as found by Fukumura et al. (2017)
for the black hole XRB GRO J1655+40. One could test this hypothesis by
systematically measuring line fluxes at various phases, especially near
eclipse. If the line emission truly originates farther away from the binary
system, the fluxes and $R$ ratios of the He-like triplets should be (near)
constant throughout the orbit, including during the eclipse. Testing this
hypothesis will be possible with XRISM/Resolve observations of Vela X-1 during
the eclipse and in phases close to eclipse. We caution that the varying nature
of the X-ray continuum of Vela X-1 is likely to affect ionized line fluxes as
well. Thus, a number of observations at each phase is required to obtain a
mean distance before meaningful conclusions can be drawn, as done here for Fe
K$\alpha$.

##  5 Conclusions

From the results presented in this work, the following conclusions can be
drawn:

  * •

NICER’s complete coverage in phase space of Vela X-1 allowed us to separately
monitor the continuum flux, the fluorescence lines, and the column density as
a function of phase.

  * •

The observed Fe K$\alpha$ line flux relative to the hard continuum between
$7-20$ keV, when compared with a simple geometrical model deviates from a
smooth orbital reflection model from the B-star, and points to a substantial
role of an accretion stream in both obscuring the B-star and reflecting X-rays
from the NS source into our line of sight.

  * •

The UV-sensitive line ratios of He-like ions from Chandra/HETG and XMM-
Newton/RGS indicate distance estimates of several stellar radii from the
B-star.

  * •

A decreasing ionization with distance from the binary suggests collimation of
the ionization cone, perhaps by magnetic fields.

R.R. is supported by an appointment to the NASA Postdoctoral Program at the
NASA Goddard Space Flight Center, through a contract with ORAU. R.R. is
supported in part by the Zuckerman STEM Leadership Program. R.R. was supported
by a Ramon scholarship from the Israeli Ministry of Science and Technology.
E.B. acknowledges support from a Center of Excellence of THE ISRAEL SCIENCE
FOUNDATION (grant No. 2752/19). This work is based in part on observations
obtained with NICER. The scientific results reported in this article are based
in part on observational data obtained from the Chandra Gratings Catalog and
Archive. This work is based in part on observations obtained with XMM-Newton,
an ESA science mission with instruments and contributions directly funded by
ESA Member States and NASA. We thank the referee

## References

  * Amato et al. (2021) Amato, R., Grinberg, V., Hell, N., et al. 2021, Astronomy & Astrophysics, 648, A105
  * Arnaud (1996) Arnaud, K. 1996, in ASP Conf., Vol. 17
  * Blondin et al. (1991) Blondin, J. M., Stevens, I. R., & Kallman, T. R. 1991, The Astrophysical Journal, 371, 684
  * Canizares et al. (2005) Canizares, C. R., Davis, J. E., Dewey, D., et al. 2005, Publications of the Astronomical Society of the Pacific, 117, 1144
  * Den Herder et al. (2001) Den Herder, J., Brinkman, A., Kahn, S., et al. 2001, Astronomy & Astrophysics, 365, L7
  * Fukumura et al. (2017) Fukumura, K., Kazanas, D., Shrader, C., et al. 2017, Nature Astronomy, 1, 1
  * Gendreau et al. (2016) Gendreau, K. C., Arzoumanian, Z., Adkins, P. W., et al. 2016, in Space telescopes and instrumentation 2016: Ultraviolet to gamma ray, Vol. 9905, SPIE, 420–435
  * George & Fabian (1991) George, I., & Fabian, A. 1991, Monthly Notices of the Royal Astronomical Society, 249, 352
  * Goldstein et al. (2004) Goldstein, G., Huenemoerder, D. P., & Blank, D. 2004, The Astronomical Journal, 127, 2310
  * Grinberg et al. (2017) Grinberg, V., Hell, N., El Mellah, I., et al. 2017, Astronomy & Astrophysics, 608, A143
  * HEASARC (2014) HEASARC, N. 2014, Astrophysics Source Code Library, ascl, 1408
  * Huenemoerder et al. (2011) Huenemoerder, D. P., Mitschang, A., Dewey, D., et al. 2011, The Astronomical Journal, 141, 129
  * Kretschmar et al. (2021) Kretschmar, P., El Mellah, I., Martínez-Núñez, S., et al. 2021, Astronomy & Astrophysics, 652, A95
  * Kreykenbohm et al. (2008) Kreykenbohm, I., Wilms, J., Kretschmar, P., et al. 2008, Astronomy & Astrophysics, 492, 511
  * Leutenegger et al. (2006) Leutenegger, M. A., Paerels, F. B., Kahn, S. M., & Cohen, D. H. 2006, The Astrophysical Journal, 650, 1096
  * Liao et al. (2020) Liao, Z., Liu, J., Zheng, X., & Gou, L. 2020, Monthly Notices of the Royal Astronomical Society, 492, 5922
  * Magdziarz & Zdziarski (1995) Magdziarz, P., & Zdziarski, A. A. 1995, Monthly Notices of the Royal Astronomical Society, 273, 837
  * Makishima et al. (1999) Makishima, K., Mihara, T., Nagase, F., & Tanaka, Y. 1999, The Astrophysical Journal, 525, 978
  * Manousakis (2011) Manousakis, A. 2011, PhD thesis, University of Geneva
  * Manousakis et al. (2012) Manousakis, A., Walter, R., & Blondin, J. 2012, Astronomy & Astrophysics, 547, A20
  * Martínez-Núñez et al. (2014) Martínez-Núñez, S., Torrejón, J. M., Kühnel, M., et al. 2014, Astronomy & Astrophysics, 563, A70
  * Mauche et al. (2008) Mauche, C. W., Liedahl, D. A., Akiyama, S., & Plewa, T. 2008, in American Institute of Physics Conference Series, Vol. 1054, Cool Discs, Hot Flows: The Varying Faces of Accreting Compact Objects, ed. M. Axelsson, 3–11, doi: 10.1063/1.3002508
  * Nagase et al. (1986) Nagase, F., Hayakawa, S., Sato, N., Masai, K., & Inoue, H. 1986, Publications of the Astronomical Society of Japan, 38, 547
  * Nagase et al. (1994) Nagase, F., Zylstra, G., Sonobe, T., et al. 1994, The Astrophysical Journal, 436, L1
  * Peretz et al. (2019) Peretz, U., Miller, J. M., & Behar, E. 2019, The Astrophysical Journal, 879, 102
  * Porquet & Dubau (2000) Porquet, D., & Dubau, J. 2000, Astronomy and Astrophysics Supplement Series, 143, 495
  * Prinja et al. (1990) Prinja, R. K., Barlow, M., & Howarth, I. D. 1990, The Astrophysical Journal, 361, 607
  * Rahin & Behar (2020) Rahin, R., & Behar, E. 2020, The Astrophysical Journal, 904, 40
  * Schulz et al. (2001) Schulz, N. S., Canizares, C. R., Lee, J. C., & Sako, M. 2001, The Astrophysical Journal Letters, 564, L21
  * Tashiro et al. (2020) Tashiro, M., Maejima, H., Toda, K., et al. 2020, in Space Telescopes and Instrumentation 2020: Ultraviolet to Gamma Ray, Vol. 11444, SPIE, 293–303
  * Watanabe et al. (2006) Watanabe, S., Sako, M., Ishida, M., et al. 2006, The Astrophysical Journal, 651, 421

##  Appendix A Spherical reflection model

Assuming a binary inclination of $90^{\circ}$, the exit path length from a
circle of radius $R$ when assuming the X-ray source (origin) is at a distance
of $1.8R$ from the center is:

$\displaystyle
y=-1.8\cos(\alpha)+\cos(\alpha+\theta^{\prime})\left(x+1.8-\sqrt{1-3.24\sin^{2}(\theta^{\prime})}\right)+$
(A1)
$\displaystyle+\bigg{(}3.6x\cos(\theta^{\prime})-x^{2}-2x\left(1.8\,-\sqrt{1-3.24\sin^{2}(\theta^{\prime})}\right)$
$\displaystyle\left(\cos(\alpha+\theta^{\prime})\left(-\sqrt{1-3.24\sin^{2}(\theta^{\prime})}+x+1.8\right)-1.8\cos(\alpha)\right)^{2}\bigg{)}^{0.5}$

Where $x$ is the length of the path into the sphere, $\alpha$ is the angle
between the line of sight and the line connecting the center of the sphere and
X-ray source, and $\theta^{\prime}$ is the angle of incidence of the entry
path. All lengths are in units of $R$. Figure 6 is the drawing this 2D model
describes.

Since the B-star itself has an average number density in excess of $10^{16}$
cm-3, we can assume $x\ll R$ and calculate $y(x)$ to first order,
$y=Cx+O\left(x^{2}\right)$:

$\displaystyle
C_{2D}=\frac{\left(1\,-0.556\sqrt{1-3.24\sin^{2}(\theta^{\prime})}\right)\cos^{2}(\alpha+\theta^{\prime})-\cos(\alpha)\cos(\alpha+\theta^{\prime})}{\left(0.556\sqrt{1-3.24\sin^{2}(\theta^{\prime})}-1\right)\cos(\alpha+\theta^{\prime})+\cos(\alpha)}+$
(A2)
$\displaystyle+\frac{0.556\sqrt{1-3.24\sin^{2}(\theta^{\prime})}+\cos(\theta^{\prime})-1}{\left(0.556\sqrt{1-3.24\sin^{2}(\theta^{\prime})}-1\right)\cos(\alpha+\theta^{\prime})+\cos(\alpha)}+\cos(\alpha+\theta^{\prime})$

Where $C_{2D}$ refers to the solution for the 2D problem.

We can now consider the same problem in 3 dimensions. If we assume the
inclination is $90^{\circ}$, the exiting photons always travel on a plane
parallel to the one in Figure 6. Therefore, all we need to do is transform
$\theta^{\prime}$ to $\theta$ and $\phi$. This adds the following 2
corrections:

$\theta^{\prime}=\arctan(\frac{\sin(\theta)\cos(\phi)}{\cos(\theta)})$ (A3)
$C(\theta,\phi)=C_{2D}\cdot\sqrt{\cos(\theta)^{2}+\sin(\theta)^{2}\cos(\phi)^{2}}$
(A4)

where $\theta$ and $\phi$ are the standard spherical coordinates with the
z-axis being the line connecting the NS and the center of the B-star.

The resulting expected flux is therefore:

$F=A_{\textrm{Z}}\omega_{K\alpha}\int\displaylimits_{E_{K,Z}}^{\infty}\left(\int\frac{I_{0}(\frac{E}{E_{0}})^{-\Gamma}\sigma_{Z,K}(E)}{\sigma(E)+C(\theta,\phi)\sigma(E_{K\alpha,Z})}d\Omega\right)dE$
(A5)

Where $\Omega$ covers the visible sphere cap from the point-of-view of the
X-ray source. The other parameters are as described in Rahin & Behar (2020).

##  Appendix B Observations and measurements

Table 1: List of high-resolution observations used Instrument | Observation ID | Date observed | Exposure (s) | Phases
---|---|---|---|---
Chandra HETG | 1926 | 2001-02-11 | 83150 | 0.98 - 0.09
Chandra HETG | 1927 | 2001-02-07 | 29430 | 0.46-0.5
Chandra HETG | 14654 | 2013-07-30 | 45880 | 0.72-0.78
Chandra HETG | 19953 | 2017-02-22 | 70470 | 0.11-0.2
Chandra HETG | 19952 | 2017-06-19 | 54330 | 0.16-0.23
XMM-Newton RGS | 0841890201 | 2019-05-03 | 113500 | 0.36-0.51
Table 2: R measurements for various Vela X-1 phases Observation ID | Phase | Ne R) | mg R | O R | Si R
---|---|---|---|---|---
1926 | 0.98 - 0.09 | $2.07_{-0.69}^{+1.52}$ | $2.48_{-0.83}^{+1.52}$ | - | $4.57_{-1.46}^{+2.67}$
1927 | 0.46-0.5 | $0.67_{-0.15}^{+0.34}$ | $1.23_{-0.28}^{+0.30}$ | - | $3.66_{-0.73}^{+1.0}$
14654 | 0.72-0.78 | $1.28_{-0.48}^{+0.99}$ | $1.36_{-0.48}^{+0.84}$ | - | $6.98_{-2.4}^{+5.93}$
19953 | 0.11-0.2 | $0.91_{-0.72}^{+2.16}$ | $0.91_{-0.25}^{+0.33}$ | - | $12.54_{-6.24}^{+49.87}$
19952 | 0.16-0.23 | $0.77_{-0.43}^{+0.71}$ | $1.81_{-0.56}^{+0.85}$ | - | $5.8_{-1.4}^{+2.26}$
0841890201 | 0.36-0.51 | $0.48_{-0.1}^{+0.12}$ | $1.42_{-0.42}^{+0.76}$ | $0.74_{-0.4}^{+0.07}$ | -


# Static Black Binaries in de Sitter

Óscar J. C. Dias<EMAIL_ADDRESS>STAG research centre & Mathematical
Sciences, University of Southampton, Highfield Campus, UK  Gary W. Gibbons
<EMAIL_ADDRESS>DAMTP, Centre for Mathematical Sciences, University of
Cambridge, Wilberforce Road, Cambridge CB3 0WA, UK  Jorge E. Santos
<EMAIL_ADDRESS>DAMTP, Centre for Mathematical Sciences, University of
Cambridge, Wilberforce Road, Cambridge CB3 0WA, UK  Benson Way
<EMAIL_ADDRESS>Departament de Física Quàntica i Astrofísica, Institut de
Ciències del Cosmos Universitat de Barcelona, Martí i Franquès, 1, E-08028
Barcelona, Spain

###### Abstract

We construct the first four-dimensional multi-black hole solution of general
relativity with a positive cosmological constant. The solution consists of two
static black holes whose gravitational attraction is balanced by the cosmic
expansion. These static binaries provide the first four-dimensional example of
non-uniqueness in general relativity without matter.

††preprint: APS/123-QED

## Introduction.

Black holes are famously featureless. This idea is embodied by the no-hair
theorems, which state in essence that stationary black holes are uniquely
characterized by their mass, angular momentum, and charge [1, 2, 3, 4, 5, 6,
7, 8, 9, 10, 11].

It should be noted that there are many situations where black hole uniqueness,
as we have expressed it, is known to be violated. A well-known example
involves multi-horizon configurations of charged, extremal black holes [12,
13]. Other examples include higher dimensions [14], anti-de Sitter asymptotics
[15, 16, 17], or exotic matter like classical Yang-Mills fields, complex
scalars and Proca fields [18, 19, 20].

Additionally, there are some mathematical gaps in fully establishing black
hole uniqueness, even in the more limited case of four-dimensional pure
gravity in flat space. Indeed, asymptotically flat multi-Kerr black holes,
where their gravitational attraction might be balanced by spin-spin
interactions, have not been ruled out (see _e.g._ [21, 22, 23, 24, 25, 26, 27,
28, 29, 30, 31, 32, 33] for attempted constructions that yield singular
configurations). Though for static solutions, a classic theorem due to [34,
35, 36] precludes the existence of regular asymptotically flat multi black
holes.

Despite these (and potentially more) counterexamples, there is currently no
experimental or observational evidence that black hole non-uniqueness can be
realized in astrophysical or cosmological contexts. Indeed, the no-hair
theorems are fully consistent with observational results from the LIGO
consortium [37].

However, the no-hair theorems assume that spacetime is asymptotically flat, a
feature which is violated in our universe at the longest scales by the
presence of a cosmological constant [38, 39, 40, 41, 42]. The resulting cosmic
expansion might balance out the gravitational attraction of two or more black
holes, allowing multiple black holes to exist in static equilibrium. Such a
configuration would share the same mass and angular momentum as some single-
horizon black hole and therefore serve as a more realistic counterexample to
black hole uniqueness.

The aim of this Letter is to demonstrate that such a multi-horizon
configuration indeed occurs within general relativity. We will focus on the
simplest case with two equal-mass black holes that do not rotate nor contain
charge, but our results and methods can be straightforwardly generalized. We
will first show that these black binaries can be anticipated using intuition
from Newton-Hooke theory, and then construct these solutions by solving the
Einstein equation numerically. Finally, we study the properties of these
binaries in detail.

Our results, along with physical intuition, suggest that the static de Sitter
binaries are dynamically unstable. Nevertheless, there remains a possibility
that they can be stabilized with the introduction of charge or angular
momentum. We will comment on this and other matters in the conclusions.

Before we continue, we mention some closely related work. Dynamical (_i.e._
out of equilibrium) multi-black holes in Einstein-Maxwell theory with a
positive cosmological constant were found in [43]. The “rod-structure”
corresponding to our static binaries were anticipated and examined in detail
in [44]. In [45], a novel mechanism for balancing multi-black holes was
proposed. These constructions provide Ricci-flat, closed-form solutions for
static binaries supported by expanding bubbles of nothing. Mechanically, these
solutions behave similarly to the static binaries we find.

Finally, we mention the mathematical papers [46, 47, 48], which might seem to
rule out the existence of static black binaries in de Sitter. We will show
that the assumptions made in [46, 47] do not apply, and that (for technical
reasons) this conclusion from [48] is not correct.

## Newton-Hooke.

Let us first set out to see if the aforementioned multi-black hole
configurations are allowed within Newtonian gravity. We adopt geometrized
units in which $c=G=k_{B}=\hbar=1$.

Consider a configuration of $N$ black holes with masses $m_{a}$, with
$a=1,\ldots,N$. For the Newtonian approximation to be valid, we assume that
the distances between the black holes are much larger than their masses. We
now include the effects of the cosmological constant $\Lambda\equiv
3/\ell^{2}>0$, where $\ell$ is the de Sitter length scale. Accordingly, we
assume that the entire configuration of black holes lies within a distance
much smaller than $\ell$ and consider the Newton-Hooke equations of motion
[49, 50]

$m_{a}\frac{\mathrm{d}^{2}\mathbf{x}_{a}}{\mathrm{d}t^{2}}-m_{a}\frac{\mathbf{x}_{a}}{\ell^{2}}=-\sum_{b\neq
a}^{b=N}\frac{\,m_{a}\,m_{b}(\mathbf{x}_{a}-\mathbf{x}_{b})}{|\mathbf{x}_{a}-\mathbf{x}_{b}|^{3}}\,,$
(1)

where $\mathbf{x}_{a}$ are the positions of the black holes.

Static solutions exist when

$\frac{\mathbf{x}_{a}}{\ell^{2}}=\sum_{b\neq
a}^{b=N}\frac{\,m_{b}(\mathbf{x}_{a}-\mathbf{x}_{b})}{|\mathbf{x}_{a}-\mathbf{x}_{b}|^{3}}\,.$
(2)

Such solutions are known as central configurations and provide homothetic
solutions of the Newtonian $N$-body problem which has applications to
Newtonian cosmology. The equation (2) can also be obtained from the Dmitriev-
Zeldovich equations [51] by using the scale factor $S(t)=e^{\frac{t}{\ell}}$,
corresponding to a de Sitter background in “stead-state” coordinates [49, 50,
52, 53].

Consider a central configuration with two equal mass black holes aligned along
the $z$ axis and separated by a distance $d$. That is, $N=2$,
$\mathbf{x}_{1}=-\mathbf{x}_{2}=\frac{d}{2}\,\mathbf{e}_{z}$, and
$m_{a}=m_{b}=M$. Then (2) imposes

$\frac{d^{3}}{\ell^{3}}=\frac{r_{+}}{\ell}$ (3)

where $r_{+}\equiv 2M$ is the Schwarzschild radius.

The requirements that the Newton-Hooke approximation should be valid and that
the black holes are inside a single cosmological event horizon amount to

$r_{+}\ll d\,,\quad d\ll\ell\quad\text{and}\quad r_{+}\ll\ell\,.$ (4)

If the distance between the black holes is given as in (3), then we see that
the first two conditions above are satisfied if we assume the third, _i.e._ if
the black holes are small enough. We therefore conclude that static de Sitter
binaries with small black holes are consistent with Newton-Hooke theory.

For later use, we introduce the event horizon Hawking temperature $T_{+}=(4\pi
r_{+})^{-1}$ and rewrite (3) as

$\frac{d}{\ell}=\frac{1}{(4\pi\ell\,T_{+})^{1/3}}\,.$ (5)

We will confirm that our numerical solutions to the Einstein equation satisfy
this scaling in the appropriate limit.

## Numerical construction.

We now construct static binaries in general relativity by numerically solving
the Einstein equation with a positive cosmological constant:

$R_{ab}=\frac{3}{\ell^{2}}\,g_{ab}\,,$ (6)

where $R_{ab}$ is the Ricci tensor and $g_{ab}$ is the metric tensor.

We use the DeTurck method, first introduced for general relativity in [54] and
reviewed in [55, 56]. This method provides a convenient way of addressing the
issue of gauge invariance, which ultimately causes the Einstein equation (6)
to yield a set of ill-posed, non-elliptic PDEs.

The DeTurck method involves choosing any reference metric $\bar{g}$ with the
same symmetries and causal structure as the solution we seek. In this case,
our reference metric is static, contains two identical black holes, a
cosmological horizon, and is axisymmetric. There is therefore a discrete
$\mathbb{Z}_{2}$ symmetry, as well as two Killing vector fields
$k=\partial/\partial t$ and $m=\partial/\partial\phi$. We further assume that
the black holes and cosmological horizon are Killing horizons generated by
$k$. Our specific choice of reference metric involves a combination of the
Israel-Khan solution [57] and the static patch of de Sitter space. Its design
is detailed in the Supplementary Material.

We then write down the most general metric ansatz $g$ that respects the
desired symmetries and causal structure. In this case, the metric ansatz
depends non-trivially on two coordinates (_i.e._ it is cohomogeneity-two, and
will yield two-dimensional PDEs).

We then solve the Einstein-DeTurck equation

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

where $\xi^{a}\equiv
g^{bc}\left[\Gamma^{a}_{bc}(g)-\Gamma^{a}_{bc}(\bar{g})\right]$, and
$\Gamma(\mathfrak{g})$ is the metric-preserving Christoffel connection
associated to a metric $\mathfrak{g}$. Unlike the Einstein equation, the
Einstein-DeTurck equation (7) yields a set of elliptic PDEs [58], which gives
a well-posed boundary-value problem with appropriate physical boundary
conditions.

The Einstein-DeTurck equation (7) is solved numerically. One complication is
that the integration domain contains five boundaries: the $\mathbb{Z}_{2}$
reflection surface, two disconnected symmetry axes, the black hole horizons,
and the cosmological horizon. We handle this domain using patching techniques.
This and other numerical methods we use are described in [56] and detailed in
the Supplementary Material.

After solving (7), we must verify that the solution actually solves the
Einstein equation, _i.e._ that $\xi=0$, and is therefore not a Ricci soliton
(for which $\xi\neq 0$). Under many circumstances [58, 59], it can be proved
that these unwanted Ricci solitons do not exist. Unfortunately, the present
case is not one of these circumstances. Indeed, with a positive cosmological
constant, Ricci solitons are known to exist (see _e.g._ [60]). Nevertheless,
ellipticity guarantees local uniqueness. That is, solutions with $\xi=0$
cannot be arbitrarily close to solutions with $\xi\neq 0$, and thus the norm
$\xi^{a}\xi_{a}$ can be monitored to identify whether our numerical
discretization converges in the continuum to a Ricci soliton or to a true
solution of the Einstein equation. In the Supplementary Material, we provide
ample evidence that the numerical solutions we construct are _not_ Ricci
solitons.

## Results.

Having numerical solutions corresponding to static black binaries in de
Sitter, we can now describe their properties and compare the numerical results
to Newton-Hooke theory when the latter is valid.

We expect to find agreement with Newton-Hooke theory when the black holes
become sufficiently small, or alternatively, when the black hole temperature
becomes sufficiently large $T_{+}\ell\gg 1$. In FIG. 1, we provide a log-log
plot of the proper distance between the horizons of the two black holes along
the symmetry axis $\mathcal{P}_{\phi}/\ell$, as a function of temperature
$4\pi T_{+}\ell$. The solid black line is the scaling (5), and the blue dots
are the numerical data. The agreement at large values of $T_{+}\ell$ shows the
validity of the Newton-Hooke analysis and corroborates our numerical
construction.

Figure 1: Proper distance between the black hole horizons versus the black
hole temperature. The solid black line shows the scaling (5) according to
Newton-Hooke analysis and the blue dots show numerical data according to
general relativity.

We have not managed to find solutions with large black holes (small $4\pi
T_{+}\ell$). Because our solutions do not have regions of large curvature,
there might be a “turning point” to a new branch of solutions. A similar
phenomenon occurs for localized Kaluza-Klein black holes when the black holes
are large relative to the Kaluza-Klein circle [61, 62, 63, 64, 65, 66, 67, 68,
69, 70, 71, 72, 73, 74, 75, 54, 55, 76, 77, 78, 79, 80]. We leave the
exploration of this region of parameter space for future work.

Let us now discuss black hole thermodynamics. For a central configuration
containing $N$ black holes inside the static patch of de Sitter, the covariant
phase space formalism [81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92] shows
that the following form of the first law of black hole mechanics holds

$\sum_{i=1}^{N}T^{(i)}_{+}\,\mathrm{d}S^{(i)}_{+}=-T_{c}\,\mathrm{d}S_{c}\,,$
(8)

where $T_{c}$ is the temperature of the cosmological horizon, and $S_{c}$ is
its entropy (_i.e_ horizon area). $T^{(i)}_{+}$ and $S_{+}^{(i)}$ are the same
quantities, respectively, for the $i-$th black hole. With $N=2$ and equal-mass
black holes, we find

$2T_{+}\,\mathrm{d}S_{+}=-T_{c}\,\mathrm{d}S_{c}\,.$ (9)

We have checked that our data satisfies this form of the first law to within
$0.01\%$.

Following [93], we now consider the entropy, which must increase during time
evolution. The blue dots in FIG. 2 show the entropy of the static binary as a
function of the entropy $S_{c}$ of the cosmological horizon. The black curve
shows the entropy for the single Schwarzschild-de Sitter black hole (also
known as the Kottler black hole). We see that for any given $S_{c}$, the
single Schwarzschild black hole has higher entropy than the binary. This
indicates that the binary can evolve towards the single black hole but not the
other way around. The static black binary is therefore thermodynamically
unstable.

The fact that (at least) two solutions exist for a given cosmological horizon
entropy implies that the Schwarzschild-de Sitter black hole is not unique.
This is the first counterexample to the no-hair conjecture [94] for pure
gravity with a positive cosmological constant.

Figure 2: Total black hole entropy versus the cosmological horizon entropy.
The blue dots are numerical data for static binaries ($S=2S_{+}$) and the
solid black line is for the single Schwarzschild-de Sitter black hole. Figure
3: Contour plot showing the level sets of the lapse function $N$. The
cosmological horizon is the outer solid black semicircle. The horizon axes has
the two black hole horizons as solid magenta lines, and the outer and inner
axes in dashed black lines. The green square is where $N$ takes its maximum
value.

We now comment on the uniqueness theorems for de Sitter black holes [46, 47,
48] which would, under certain assumptions, rule out the existence of static
de Sitter binaries. In [46], the level sets of the lapse function
$N\equiv\sqrt{-g_{tt}}$ are assumed to be surface forming. In particular, this
means that the level sets must consist only of 2D surfaces. In [47], the set
$\mathrm{MAX}(N)=\\{x\in\mathcal{M}\mathrel{\mathop{\mathchar
58\relax}}N(x)=N_{\max}\\}$, where $N_{\max}$ is the maximum value of $N$ in
the manifold $\mathcal{M}$, is assumed to disconnect $\mathcal{M}$ into an
inner region $\mathcal{M}_{-}$ and an outer region $\mathcal{M}_{+}$ with the
same virtual mass. Our static binaries do not satisfy either of these
assumptions. Indeed, in FIG. 3 we show the level sets $N$ in our domain of
integration for a typical solution (all of our solutions show the same
qualitative behavior). The coordinates $(r,z)$ are defined in the Supplemental
Material. The cosmological horizon is represented by the outer solid black
semicircle, the two black hole horizons are marked by solid magenta lines
along the horizontal axis, and the inner and outer axes are given by the
dashed horizontal line. Finally, the green square marks the location of the
maximum of $N$ in $\mathcal{M}$. This maximum represents an $S^{1}$ on the
manifold, which is not a 2D surface, and it also does not partition the
manifold into two regions. Therefore, our static binaries fail to satisfy the
assumptions in [46, 47].

Finally, we comment on [48]. We believe that this work is not correct for a
rather technical reason. Beginning with the Schwarzschild-de Sitter black
hole, the authors in [48] argue that they can construct an asymptotically flat
metric that is conformal to the original one, is topologically $S^{1}\times
S^{2}$ deprived of one point, and has zero ADM mass. If that were true, the
rigidity statement in the positive mass theorem [95, 96, 97, 98, 99] would not
only imply that the original metric is conformally flat, but also that
$S^{1}\times S^{2}$ with one point removed is diffeomorphic to
$\mathbb{R}^{3}$, which is impossible.

## Conclusions.

We constructed the first example of a multi-black hole solution within general
relativity with a positive cosmological constant and established that the
leading behavior of these solutions agrees (for small black holes) with
estimates from Newton-Hooke theory. Based on thermodynamic considerations, we
argued that these solutions are thermodynamically unstable. Because the
configuration requires a delicate balance between gravitational attraction and
cosmic expansion, we expect these solutions to also be dynamically unstable.

We have focused on the static configuration of two identical black holes, but
our results and methods can be generalized. First, consider the case where the
black holes have different masses. When one of the black holes is much smaller
than the other, one can use the geodesic approximation to predict the
existence of such a configuration. Indeed, one can easily confirm the
existence of static orbits for timelike particles on a Schwarzschild-de Sitter
black hole background, thus providing further evidence for the existence of
this more general central configuration. Note that if [48] were correct, this
asymmetric binary would also not exist.

We can also include rotation, which will introduce spin-spin interaction of
the black holes. This opens the possibility of continuous non-uniqueness.
Consider, for example, the case with two identical black holes rotating in
opposite directions along the axis of symmetry. This configuration will have
vanishing total angular momentum, and will thus be in the same class as the
Schwarzschild-de Sitter black hole. Work in this direction is underway.

Perhaps more interestingly, because spin-spin interactions act on shorter
length scales, they could provide a mechanism for stabilizing the binary. This
possibility resembles the mechanism that provides stability for molecules.
Work in this direction is underway.

We could also consider central configurations containing $N>2$ static black
holes in the static patch of de Sitter. These configurations can show
interesting properties within the Newton-Hooke approximation. For instance,
when $N\geq 13$, minimal energy central configurations do not lie on a regular
polyhedron [49]. We thus expect the equivalent property within general
relativity. The study of these configurations is within the reach of the
numerical methods employed in this Letter.

###### Acknowledgements.

The authors would like to thank Stefano Borghini, Lorenzo Mazzieri and Ryan
Unger for valuable discussions and correspondence regarding [48]. J. E. S.
would like to thank Jay Armas and Troels Harmark for mentioning this
interesting problem to him in 2012 at the 28th Pacific Coast Gravity Meeting,
held at University of California, Santa Barbara and Don Marolf for discussions
on the microcanonical properties of black holes contained inside the static
patch of de Sitter. O. C. D. acknowledges financial support from the STFC
“Particle Physics Grants Panel (PPGP) 2018” Grant No. ST/T000775/1. J. E. S.
has been partially supported by STFC consolidated grant ST/T000694/1. B. W.
would like to thank Roberto Emparan for fruitful discussions and acknowledges
support from ERC Advanced Grant GravBHs-692951 and MEC grant
FPA2016-76005-C2-2-P. 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.

## References

  * [1] J. E. Chase, “Event horizons in static scalar-vacuum space-times,” Commun. Math. Phys. 19 (1970) 276.
  * [2] R. Penney, “Axially Symmetric Zero-Mass Meson Solutions of Einstein Equations,” Phys. Rev. 174 (1968) 1578–1579.
  * [3] J. D. Bekenstein, “Transcendence of the law of baryon-number conservation in black hole physics,” Phys. Rev. Lett. 28 (1972) 452–455.
  * [4] J. D. Bekenstein, “Nonexistence of baryon number for static black holes,” Phys. Rev. D5 (1972) 1239–1246.
  * [5] J. D. Bekenstein, “Nonexistence of baryon number for black holes. ii,” Phys. Rev. D5 (1972) 2403–2412.
  * [6] C. Teitelboim, “Nonmeasurability of the quantum numbers of a black hole,” Phys. Rev. D5 (1972) 2941–2954.
  * [7] J. B. Hartle, Magic Without Magic.  edited by J. Klauder, Freeman, San Francisco, 1972.
  * [8] M. Heusler, “A No hair theorem for selfgravitating nonlinear sigma models,” J. Math. Phys. 33 (1992) 3497–3502.
  * [9] J. D. Bekenstein, “Novel ’no scalar hair’ theorem for black holes,” Phys. Rev. D51 (1995) 6608–6611.
  * [10] 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.
  * [11] D. Sudarsky, “A Simple proof of a no hair theorem in Einstein Higgs theory,,” Class. Quant. Grav. 12 (1995) 579–584.
  * [12] S. D. Majumdar, “A class of exact solutions of Einstein’s field equations,” Phys. Rev. 72 (1947) 390–398.
  * [13] 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.
  * [14] R. Emparan and H. S. Reall, “A Rotating black ring solution in five-dimensions,” Phys. Rev. Lett. 88 (2002) 101101, arXiv:hep-th/0110260.
  * [15] S. S. Gubser, “Breaking an Abelian gauge symmetry near a black hole horizon,” Phys. Rev. D 78 (2008) 065034, arXiv:0801.2977 [hep-th].
  * [16] S. A. Hartnoll, C. P. Herzog, and G. T. Horowitz, “Holographic Superconductors,” JHEP 12 (2008) 015, arXiv:0810.1563 [hep-th].
  * [17] O. J. C. Dias, G. T. Horowitz, and J. E. Santos, “Black holes with only one Killing field,” JHEP 07 (2011) 115, arXiv:1105.4167 [hep-th].
  * [18] M. S. Volkov and D. V. Gal’tsov, “Gravitating nonAbelian solitons and black holes with Yang-Mills fields,” Phys. Rept. 319 (1999) 1–83, arXiv:hep-th/9810070 [hep-th].
  * [19] C. A. R. Herdeiro and E. Radu, “Kerr black holes with scalar hair,” Phys. Rev. Lett. 112 (2014) 221101, arXiv:1403.2757 [gr-qc].
  * [20] C. Herdeiro, E. Radu, and H. Rúnarsson, “Kerr black holes with Proca hair,” Class. Quant. Grav. 33 no. 15, (2016) 154001, arXiv:1603.02687 [gr-qc].
  * [21] Z. Perjés, “Solutions of the coupled einstein-maxwell equations representing the fields of spinning sources,” Phys. Rev. Lett. 27 (Dec, 1971) 1668–1670.
  * [22] W. Israel and G. A. Wilson, “A class of stationary electromagnetic vacuum fields,” J. Math. Phys. 13 (1972) 865–871.
  * [23] J. B. Hartle and S. W. Hawking, “Solutions of the Einstein-Maxwell equations with many black holes,” Commun. Math. Phys. 26 (1972) 87–101.
  * [24] A. Tomimatsu and H. Sato, “New Exact Solution for the Gravitational Field of a Spinning Mass,” Phys. Rev. Lett. 29 (1972) 1344–1345.
  * [25] D. Kramer and G. Neugebauer, “The superposition of two kerr solutions,” Physics Letters A 75 no. 4, (1980) 259–261.
  * [26] G. Neugebauer, “Relarivistic gravitational fields of rotating bodies,” Physics Letters A 86 no. 2, (1981) 91–93.
  * [27] 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.
  * [28] V. S. Manko and E. Ruiz, “Exact solution of the double-kerr equilibrium problem,” Classical and Quantum Gravity 18 no. 2, (Jan, 2001) L11.
  * [29] C. A. R. Herdeiro and C. Rebelo, “On the interaction between two Kerr black holes,” JHEP 10 (2008) 017, arXiv:0808.3941 [gr-qc].
  * [30] 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].
  * [31] 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].
  * [32] V. S. Manko and E. Ruiz, “Metric for two arbitrary Kerr sources,” Phys. Lett. B 794 (2019) 36–40, arXiv:1806.10408 [gr-qc].
  * [33] 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].
  * [34] W. Israel, “Event horizons in static vacuum space-times,” Phys. Rev. 164 (1967) 1776–1779.
  * [35] D. C. Robinson, “A simple proof of the generalization of Israel’s theorem.,” General Relativity and Gravitation 8 no. 8, (Aug., 1977) 695–698.
  * [36] 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.
  * [37] KAGRA, LIGO Scientific, Virgo, VIRGO Collaboration, B. P. Abbott et al., “Prospects for observing and localizing gravitational-wave transients with Advanced LIGO, Advanced Virgo and KAGRA,” Living Rev. Rel. 21 no. 1, (2018) 3, arXiv:1304.0670 [gr-qc].
  * [38] 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.
  * [39] 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.
  * [40] 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.
  * [41] A. Wright, “It’s a mystery,” Nature Physics 7 no. 11, (Nov, 2011) 833–833.
  * [42] 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].
  * [43] D. Kastor and J. H. Traschen, “Cosmological multi - black hole solutions,” Phys. Rev. D 47 (1993) 5370–5375, arXiv:hep-th/9212035.
  * [44] J. Armas, P. Caputa, and T. Harmark, “Domain Structure of Black Hole Space-Times with a Cosmological Constant,” Phys. Rev. D 85 (2012) 084019, arXiv:1111.1163 [hep-th].
  * [45] 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].
  * [46] P. G. LeFloch and L. Rozoy, “Uniqueness of kottler spacetime and the besse conjecture,” Comptes Rendus Mathematique 348 no. 19, (2010) 1129–1132.
  * [47] S. Borghini, P. T. Chruściel, and L. Mazzieri, “On the uniqueness of Schwarzschild-de Sitter spacetime,” arXiv:1909.05941 [math.DG].
  * [48] A. M. ul Alam and W. Yu, “Uniqueness of de sitter and schwarzschild-de sitter spacetimes,” Communications in Analysis and Geometry 23 no. 2, (2014) 377–387.
  * [49] R. A. Battye, G. W. Gibbons, and P. M. Sutcliffe, “Central configurations in three-dimensions,” Proc. Roy. Soc. Lond. A 459 (2003) 911–943, arXiv:hep-th/0201101.
  * [50] G. W. Gibbons and C. E. Patricot, “Newton-Hooke space-times, Hpp waves and the cosmological constant,” Class. Quant. Grav. 20 (2003) 5225, arXiv:hep-th/0308200.
  * [51] I. Dmitriev and V. Nikolaev, “Semi-empirical method for the calculation of the equilibrium distribution of charges in a fast-ion beam,” Sov. Phys. JETP 20 (1965) 409–415.
  * [52] G. F. R. Ellis and G. W. Gibbons, “Discrete Newtonian Cosmology,” Class. Quant. Grav. 31 (2014) 025003, arXiv:1308.1852 [astro-ph.CO].
  * [53] G. F. R. Ellis and G. W. Gibbons, “Discrete Newtonian Cosmology: Perturbations,” Class. Quant. Grav. 32 no. 5, (2015) 055001, arXiv:1409.0395 [gr-qc].
  * [54] 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].
  * [55] 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.
  * [56] 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].
  * [57] 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.
  * [58] 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].
  * [59] 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].
  * [60] M. Feldman, T. Ilmanen, and D. Knopf, “Rotationally Symmetric Shrinking and Expanding Gradient Kähler-Ricci Solitons,” Journal of Differential Geometry 65 no. 2, (2003) 169 – 209.
  * [61] B. Kol, “Topology change in general relativity, and the black hole black string transition,” JHEP 10 (2005) 049, arXiv:hep-th/0206220 [hep-th].
  * [62] T. Wiseman, “Static axisymmetric vacuum solutions and nonuniform black strings,” Class. Quant. Grav. 20 (2003) 1137–1176, arXiv:hep-th/0209051 [hep-th].
  * [63] B. Kol and T. Wiseman, “Evidence that highly nonuniform black strings have a conical waist,” Class. Quant. Grav. 20 (2003) 3493–3504, arXiv:hep-th/0304070 [hep-th].
  * [64] T. Harmark, “Small black holes on cylinders,” Phys. Rev. D69 (2004) 104015, arXiv:hep-th/0310259 [hep-th].
  * [65] D. Gorbonos and B. Kol, “A Dialogue of multipoles: Matched asymptotic expansion for caged black holes,” JHEP 06 (2004) 053, arXiv:hep-th/0406002 [hep-th].
  * [66] T. Harmark and N. A. Obers, “New phases of near-extremal branes on a circle,” JHEP 09 (2004) 022, arXiv:hep-th/0407094 [hep-th].
  * [67] V. Asnin, B. Kol, and M. Smolkin, “Analytic evidence for continuous self similarity of the critical merger solution,” Class. Quant. Grav. 23 (2006) 6805–6827, arXiv:hep-th/0607129.
  * [68] T. Harmark and N. A. Obers, “Black holes on cylinders,” JHEP 05 (2002) 032, arXiv:hep-th/0204047 [hep-th].
  * [69] T. Wiseman, “From black strings to black holes,” Class. Quant. Grav. 20 (2003) 1177–1186, arXiv:hep-th/0211028.
  * [70] H. Kudoh and T. Wiseman, “Properties of Kaluza-Klein black holes,” Prog. Theor. Phys. 111 (2004) 475–507, arXiv:hep-th/0310104 [hep-th].
  * [71] H. Kudoh and T. Wiseman, “Connecting black holes and black strings,” Phys. Rev. Lett. 94 (2005) 161102, arXiv:hep-th/0409111 [hep-th].
  * [72] E. Sorkin, “Nonuniform black strings in various dimensions,” Phys. Rev. D74 (2006) 104027, arXiv:gr-qc/0608115.
  * [73] B. Kleihaus, J. Kunz, and E. Radu, “New nonuniform black string solutions,” JHEP 06 (2006) 016, arXiv:hep-th/0603119 [hep-th].
  * [74] T. Harmark, V. Niarchos, and N. A. Obers, “Instabilities of black strings and branes,” Class. Quant. Grav. 24 (2007) R1–R90, arXiv:hep-th/0701022 [hep-th].
  * [75] O. J. C. Dias, T. Harmark, R. C. Myers, and N. A. Obers, “Multi-black hole configurations on the cylinder,” Phys. Rev. D76 (2007) 104025, arXiv:0706.3645 [hep-th].
  * [76] P. Figueras, K. Murata, and H. S. Reall, “Stable non-uniform black strings below the critical dimension,” JHEP 11 (2012) 071, arXiv:1209.1981 [gr-qc].
  * [77] G.T. Horowitz, et. al, Black Holes in Higher Dimensions.  Cambridge University Press, 2012.
  * [78] M. Kalisch and M. Ansorg, “Pseudo-spectral construction of non-uniform black string solutions in five and six spacetime dimensions,” Class. Quant. Grav. 33 no. 21, (2016) 215005, arXiv:1607.03099 [gr-qc].
  * [79] O. J. C. Dias, J. E. Santos, and B. Way, “Localised $AdS_{5}\times S^{5}$ Black Holes,” Phys. Rev. Lett. 117 no. 15, (2016) 151101, arXiv:1605.04911 [hep-th].
  * [80] O. J. C. Dias, J. E. Santos, and B. Way, “Localised and nonuniform thermal states of super-Yang-Mills on a circle,” JHEP 06 (2017) 029, arXiv:1702.07718 [hep-th].
  * [81] R. M. Wald, “Black hole entropy is the noether charge,” Physical Review D 48 no. 8, (1993) R3427.
  * [82] V. Iyer and R. M. Wald, “Some properties of the noether charge and a proposal for dynamical black hole entropy,” Physical review D 50 no. 2, (1994) 846.
  * [83] V. Iyer and R. M. Wald, “Comparison of the noether charge and euclidean methods for computing the entropy of stationary black holes,” Physical Review D 52 no. 8, (1995) 4430.
  * [84] R. M. Wald and A. Zoupas, “General definition of “conserved quantities” in general relativity and other theories of gravity,” Physical Review D 61 no. 8, (2000) 084027.
  * [85] I. Papadimitriou and K. Skenderis, “Thermodynamics of asymptotically locally ads spacetimes,” Journal of High Energy Physics 2005 no. 08, (2005) 004.
  * [86] I. M. Anderson and C. G. Torre, “Asymptotic conservation laws in field theory,” Phys. Rev. Lett. 77 (1996) 4109–4113.
  * [87] G. Barnich and F. Brandt, “Covariant theory of asymptotic symmetries, conservation laws and central charges,” Nucl. Phys. B633 (2002) 3–82.
  * [88] G. Barnich, “Boundary charges in gauge theories: Using Stokes theorem in the bulk,” Class. Quant. Grav. 20 (2003) 3685–3698.
  * [89] G. Barnich and G. Compere, “Surface charge algebra in gauge theories and thermodynamic integrability,” J. Math. Phys. 49 (2008) 042901.
  * [90] G. Compere, “Note on the First Law with p-form potentials,” Phys. Rev. D75 (2007) 124020.
  * [91] D. D. K. Chow and G. Compere, “Dyonic AdS black holes in maximal gauged supergravity,” Phys. Rev. D89 no. 6, (2014) 065003.
  * [92] G. Compere, Symmetries and conservation laws in Lagrangian gauge theories with applications to the mechanics of black holes and to gravity in three dimensions.  PhD thesis, Brussels U., 2007.  arXiv:0708.3153 [hep-th].
  * [93] G. W. Gibbons and S. W. Hawking, “Cosmological Event Horizons, Thermodynamics, and Particle Creation,” Phys. Rev. D 15 (1977) 2738–2751.
  * [94] R. Ruffini and J. A. Wheeler, “Introducing the black hole,” Phys. Today 24 no. 1, (1971) 30.
  * [95] R. Schon and S.-T. Yau, “On the Proof of the positive mass conjecture in general relativity,” Commun. Math. Phys. 65 (1979) 45–76.
  * [96] R. Schon and S.-T. Yau, “Proof of the positive mass theorem. 2.,” Commun. Math. Phys. 79 (1981) 231–260.
  * [97] E. Witten, “A Simple Proof of the Positive Energy Theorem,” Commun. Math. Phys. 80 (1981) 381.
  * [98] G. W. Gibbons, S. W. Hawking, G. T. Horowitz, and M. J. Perry, “Positive Mass Theorems for Black Holes,” Commun. Math. Phys. 88 (1983) 295.
  * [99] D. A. Lee, M. Lesourd, and R. Unger, “Density and Positive Mass Theorems for Initial Data Sets with Boundary,” Commun. Math. Phys. 395 no. 2, (2022) 643–677, arXiv:2112.12017 [math.DG].

## Supplementary Material

##  Appendix A The Israel-Khan Solution

The Israel-Khan solution [57] is an exact solution that describes two
asymptotically flat black holes separated by a conical strut. Because it is an
exact multi-horizon solution, it will be useful for our numerical construction
of static de Sitter binaries. Here, we review this solution and describe some
of its coordinate representations that we use.

Aside from the conical strut, the Israel-Khan spacetime is completely regular
outside the horizons. Here, we are only concerned with the case where both
black holes are equal. The solution is often presented in Weyl coordinates:

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

where

$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)\;,$
(11) $\displaystyle h$
$\displaystyle=\left(\frac{r^{2}+(z+1/k)(z+1)+R_{+}r_{+}}{2R_{+}r_{+}}\right)\left(\frac{r^{2}+(z-1/k)(z-1)+R_{-}r_{-}}{2R_{-}r_{-}}\right)$
$\displaystyle\qquad\times\left(\frac{r^{2}+z^{2}-1+r_{+}r_{-}}{r^{2}+(z-1/k)(z+1)+r_{+}R_{-}}\right)\left(\frac{r^{2}+z^{2}-(1/k^{2})+R_{+}R_{-}}{r^{2}+(z+1/k)(z-1)+R_{+}r_{-}}\right)\;,$
(12)

with

$R_{\pm}=\sqrt{r^{2}+\left(z\pm\frac{1}{k}\right)^{2}}\;,\qquad
r_{\pm}=\sqrt{r^{2}+(z\pm 1)^{2}}\;,$ (13)

and where $\ell$ is an arbitrary length scale that we have introduced for
later use in de Sitter. The solution is parametrized by $\lambda\in(0,\infty)$
and $k\in(0,1)$. The temperature of the black holes is given by

$T=\frac{1}{2\pi}\frac{4\lambda(1-k)}{k(1+k)}\,.$ (14)

A peculiarity of the Weyl form is that the axes and horizons of the solution
are all located at $r=0$, and so is described as a “rod structure.” At $r=0$,
the horizons lie in the regions $z\in(1,1/k)$ and $z\in(-1/k,-1)$, with the
inner axis between the black holes in the region $z\in(-1,1)$, and the outer
axes in $z\in(1/k,\infty)$ and $z\in(-\infty,-1/k)$. The inner axes contains a
conical singularity which holds the two black holes apart.

The Weyl form is useful for obtaining this solution as it provides a means of
simplifying the Einstein equation into an integrable form. But in order to
accommodate the rod structure, the coordinates cannot be smooth along the line
$r=0$ (there are coordinate singularities at the “joints” between rods),
making the Weyl form ill-suited for our numerical purposes.

We therefore seek a coordinate transformation that maps the outer axes, inner
axis, and horizons into a coordinate rectangle. This can be accomplished by a
conformal Schwarz-Christoffel transformation. The standard formulas for this
type of transformation will give mappings that use Jacobi Elliptic functions,
but here we convert these functions to a more algebraic form. The mappings we
use are defined by

$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})}\;.$
(15)

Lines of constant $x$ and $y$, along with the rod structure of Israel-Khan are
shown in FIG. 4.

Figure 4: Lines of constant $x$ and $y$, shown in Weyl coordinates. There is a
reflection symmetry about $z=0$, so only positive $z$ is shown. The rod
structure of the Israel-Khan solution can be seen at $r=0$, with one of the
black hole horizons lying between the two red dots.

With this coordinate transformation, the Israel-Khan solution takes the form

$\mathrm{d}s^{2}=\ell^{2}\left\\{-f\mathrm{d}t^{2}+\frac{\lambda^{2}}{m^{2}\Delta_{xy}^{2}}\left[w^{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\\}\;,$
(16)

where in these coordinates

$f=(1-x^{2})^{2}\Delta_{x}m^{2}$ (17)

and

$\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})\;,$
(18) $w=\frac{k}{(1+k)^{2}}\left(1+\sqrt{\Delta_{y}}\right)^{2}\;,\qquad
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{]}}\;.$
(19)

The horizons are at $x=\pm 1$, the inner axis is at $y=0$, and the outer axes
are at $y=1$. There is also a $\mathbb{Z}_{2}$ symmetry about $x=0$. All
functions $\Delta_{x}$, $\Delta_{y}$, $w$, 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.

Eventually, we wish to join the Israel-Khan solution with a de Sitter horizon.
In anticipation of doing so, we present the Israel-Khan solution in polar-Weyl
coordinates defined by

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

where the Israel-Khan solution takes the form

$\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\\}\;,$
(21)

with $f$ and $h$ also transformed accordingly. For later use, we find it
convenient to express $f$ and $h$ as functions of $\rho$ and
$z=\rho\xi\sqrt{2-\xi^{2}}$:

$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)\;,$
(22) $\displaystyle h$
$\displaystyle=\left(\frac{\rho^{2}+\tfrac{1}{k}[1+(1+k)z]+R_{+}r_{+}}{2R_{+}r_{+}}\right)\left(\frac{\rho^{2}+\tfrac{1}{k}[1-(1+k)z]+R_{-}r_{-}}{2R_{-}r_{-}}\right)$
$\displaystyle\qquad\times\left(\frac{\rho^{2}-1+r_{+}r_{-}}{\rho^{2}-\tfrac{1}{k}[1+(1-k)z]+r_{+}R_{-}}\right)\left(\frac{\rho^{2}-(1/k^{2})+R_{+}R_{-}}{\rho^{2}-\tfrac{1}{k}[1-(1-kz)]+R_{+}r_{-}}\right)\;,$
(23)

with

$R_{\pm}=\sqrt{\rho^{2}+\frac{1}{k^{2}}\pm\frac{2z}{k}}\;,\qquad
r_{\pm}=\sqrt{\rho^{2}+1\pm 2z}\;,$ (24)

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

From (15) and (20), we can derive an explicit coordinate transformation
between the Schwarz-Christoffel $(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})}}\;,$
(25a) $\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})}}}\;.$
(25b)

Like the original Weyl coordinates, these polar-Weyl coordinates are not
smooth along $\xi=\pm 1$. This will not be an issue for us as we will only use
these coordinates in our numerical construction for sufficiently large $\rho$,
where the coordinates are smooth.

##  Appendix B Designing the Reference Metric for Einstein-DeTurck

Our strategy for designing a reference metric for the Einstein-DeTurck problem
(described in the main text) is to attach a de Sitter horizon to the Israel-
Khan solution. De Sitter space in four dimensions is most commonly written in
the form

$\mathrm{d}s^{2}=-\left(1-\frac{R^{2}}{\ell^{2}}\right)\mathrm{d}\tau^{2}+\frac{\mathrm{d}R^{2}}{1-\frac{R^{2}}{\ell^{2}}}+r^{2}(\mathrm{d}\theta^{2}+\sin^{2}\theta\mathrm{d}\phi^{2})\;,$
(26)

where $\ell$ is the de Sitter length scale. De Sitter can also be written in
isotropic coordinates with the transformation

$\frac{R}{\ell}=\frac{\lambda\rho}{1+\frac{\lambda^{2}\rho^{2}}{4}}\;,\qquad\sin\theta=1-\xi^{2}\;,\qquad\tau=\ell\,t\;,$
(27)

which yields

$\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\\}\;,$
(28)

where

$g_{\pm}=1\pm\frac{\lambda^{2}\rho^{2}}{4}\;.$ (29)

In these coordinates, the de Sitter horizon has a constant temperature of
$1/2\pi$. $\lambda$ is a gauge parameter that merely scales the radial
coordinate $\rho$. There is an origin at $\rho=0$, the de Sitter horizon is
located at $\rho=2/\lambda$, there is an axis of symmetry at $\xi=\pm 1$ and a
$\mathbb{Z}_{2}$ symmetry at $\xi=0$.

This form of de Sitter is suggestively close to the Israel-Khan solution in
polar-Weyl form (21). 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. We will make use of these similarities in our construction.

Now to begin engineering a reference metric, we make some slight modifications
to the Israel-Khan solution (16) and (21):

$\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[w^{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\mathrm{d}\phi^{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\mathrm{d}\phi^{2}\right)\right]\right\\}\;.$
(30)

Here,

$S=1-\alpha(1-y^{2})^{2}\;,$ (31)

where $\alpha$ is a new parameter, $F$ is a complicated function that we will
describe later in (35), and equality between the first and second lines of (B)
(here and in the remainder of this section) is understood to be through the
coordinate transformation (25). We will use $(x,y)$ coordinates in the region
near the black holes and inner axis, and the $(\rho,\xi)$ coordinates near the
cosmological horizon.

We only made three changes to the Israel-Khan solution to arrive at the
reference metric (B). The first is the inclusion of a conformal factor
$1/g_{+}^{2}$ to facilitate the matching to de Sitter. The second is a factor
of $S$ in the $\mathrm{d}\phi^{2}$ term, which we will use to remove the
conical singularity in the inner axis by adjusting the parameter $\alpha$. The
conical singularity is removed when $\alpha$ takes the value

$\alpha=\frac{(1-k)^{2}\left(k^{2}+6k+1\right)}{(k+1)^{4}}\,.$ (32)

The last change is a factor of $g_{-}^{2}\,F$ in the $\mathrm{d}t^{2}$ term
which introduces a cosmological horizon.

We have freedom to choose the function $F$, but the choice is a delicate
matter. For numerical purposes, we wish for $F$ to be smooth in $(x,y)$ or
$(\rho,\xi)$ coordinates, depending on where the coordinates are being used.
The DeTurck method also requires that $F$ be chosen to preserve the regularity
of both the cosmological horizon and of the black hole horizons [54, 55, 56].
That is, $F$ must be positive definite and satisfy

$F|_{x=\pm 1}=\frac{1}{g_{-}^{2}}|_{x=\pm 1}\;,\qquad
F|_{\rho=2/\lambda}=\frac{h}{f^{2}}|_{\rho=2/\lambda}\;,$ (33)

where we have chosen equality in the above instead of proportionality in order
to preserve the de Sitter and Israel-Khan temperatures.

In order to make it easier to find a solution in a Newton-Raphson algorithm
(see _e.g._ [56]), we should 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 Israel-Khan
solution (when $\alpha=0$). The cosmological horizon is already accommodated
by the fact that $f$ and $h$ approach unity for large $\rho$. That is, by
requiring (33), we already have $F\approx 1$ near the cosmological horizon
$\rho\approx 2/\lambda$ when $\lambda$ is small.

As for near the origin, we add the requirement that

$F|_{y=0}=\frac{1}{g_{-}^{2}}|_{y=0}\;,$ (34)

which is consistent with (33). If we set $\alpha=0$, then when $\lambda$ and
$\rho$ are small, $g_{\pm}\approx 1$, and the metric approaches that of the
Israel-Khan solution as desired.

All of these requirements can be satisfied by choosing $F$ to take the form

$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}}\;,$
(35)

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
(22) and (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{g}$. 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 entire Einstein-DeTurck reference metric is parametrized by $k$,
$\lambda$, and $\alpha$. On the desired solution that is free of conical
singularities, $\alpha$ is fixed according to (32). Because the de Sitter
horizon temperature is held fixed, the only remaining physical parameter is
the black hole horizon temperature, which is given by (14). Any $k$ and
$\lambda$ that give the same black hole temperature are physically equivalent.
To collect all our numerical data, we have fixed $\lambda=1/10$ and used $k$
to parametrize our solutions. We have tried different values of $\lambda$, but
after trial and error, this value generated the best numerical results.

##  Appendix C Patching and Numerical Methods

In this section we explain how we partition the domain of integration using
patching techniques (see _e.g._ [56]). As mentioned in the main text, the
solution we seek contains five boundaries: the inner axis
($\partial_{\phi}^{\rm in}$), the black hole horizon ($\mathcal{H}^{+}$), the
outer axis ($\partial_{\phi}^{\rm out}$), the cosmological horizon
($\mathcal{H}^{c}$) and the plane of $\mathbb{Z}_{2}$ symmetry.

We numerically solve the Einstein-DeTurck equations with the boundary
conditions of regularity at each of our five boundaries. To do so, we use both
the coordinates $(\rho,\xi)$ and $(x,y)$ of the preceding section. Near the
black hole event horizon we use $(x,y)$ coordinates, while near the
cosmological horizon we use $(\rho,\xi)$ coordinates. We used a total of three
patches $-$ $I$, $II$ and $III$ $-$ with each of these having four boundaries:
see FIG. 5. Patches $I$ and $II$ are defined in $(x,y)$ coordinates, and patch
$III$ in $(\rho,\xi)$ coordinates. The patching boundary between the patch $I$
and $II$ (dashed line in FIG. 5) is given by $x=x_{0}y\sqrt{2-y^{2}}$, with
$y\in(0,1)$. The patching boundary between patch $II$ and patch $III$ is
simply given by $\rho=\rho_{0}$.

Figure 5: Illustration of the patches used in our numerical construction.

Furthermore, we fix $\rho_{0}$ and $x_{0}$ through

$\rho_{0}=\frac{2}{3}\left(\frac{2}{\lambda}-\frac{1}{k}\right)+\frac{1}{k}\quad\text{and}\quad
x_{0}=\frac{1}{2}\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}}}}\,,$
(36)

so that the only free parameters are $k$,$\lambda$, and $\alpha$. Note that
for non-singular solutions we require $\alpha$ to be given by (32) and, as
described in the end of the previous section, we fix $\lambda=1/10$ for
numerical convenience. We now apply the numerical methods detailed in [56],
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)\times N$. The discretized Newton-Raphson equation
reduces to an iteration of linear matrix problems, which we solve by LU
decomposition. To find the first solution (which is always the trickiest step
of the Newton-Raphson method), we have set $\lambda=1/10$, $k=1/2$ and used a
Newton-Raphson algorithm with a judiciously chosen damping factor.
Additionally, we used what we call the “$\delta$-trick” as explained in
section VII.A of [56].

##  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 pseudospectral
collocation methods.

Let $\chi^{(N)}$ be $\chi$ computed on a (3-patched) grid with $(N+N+N)\times
N$ spectral collocation points. For concreteness we take $k=1/2$,
$\lambda=1/10$ and $\alpha$ as given in (32). In FIG. 6 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)=27.89409-0.79076\,N\,.$ (37)

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).

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