# Wide post-common envelope binaries containing ultramassive white dwarfs:
evidence for efficient envelope ejection in massive AGB stars

Natsuko Yamaguchi,1 Kareem El-Badry1, Jim Fuller1, David W. Latham2, Phillip
A. Cargile2, Tsevi Mazeh3, Sahar Shahaf4, Allyson Bieryla2, Lars A. Buchhave6,
Melissa Hobson5
1 Department of Astronomy, California Institute of Technology, 1200 E.
California Blvd, Pasadena, CA, 91125, USA
2 Center for Astrophysics $|$ Harvard & Smithsonian, 60 Garden Street,
Cambridge, MA 02138, USA
3 School of Physics and Astronomy, Tel Aviv University, Tel Aviv, 6997801,
Israel
4 Department of Particle Physics and Astrophysics, Weizmann Institute of
Science, Rehovot 7610001, Israel
6 DTU Space, National Space Institute, Technical University of Denmark,
Elektrovej 328, DK-2800 Kgs. Lyngby, Denmark
5 Max-Planck-Institute for Astronomy, Königstuhl 17, 69117 Heidelberg, Germany
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Post-common-envelope binaries (PCEBs) containing a white dwarf (WD) and a
main-sequence (MS) star can constrain the physics of common envelope evolution
and calibrate binary evolution models. Most PCEBs studied to date have short
orbital periods ($P_{\rm orb}\lesssim 1$ d), implying relatively inefficient
harnessing of binaries’ orbital energy for envelope expulsion. Here, we
present follow-up observations of five binaries from Gaia DR3 containing
solar-type MS stars and probable ultramassive WDs ($M\gtrsim 1.2\,M_{\odot}$)
with significantly wider orbits than previously known PCEBs, $P_{\rm
orb}=18-49$ d. The WD masses are much higher than expected for systems formed
via stable mass transfer at these periods, and their near-circular orbits
suggest partial tidal circularization when the WD progenitors were giants.
These properties strongly suggest that the binaries are PCEBs. Forming PCEBs
at such wide separations requires highly efficient envelope ejection, and we
find that the observed periods can only be explained if a significant fraction
of the energy released when the envelope recombines goes into ejecting it.
Using 1D stellar evolution calculations, we show that the binding energy of
massive AGB star envelopes is formally positive if recombination energy is
included in the energy budget. This suggests that the star’s envelope can be
efficiently ejected if binary interaction causes it to recombine, and that a
wide range of PCEB orbital periods can potentially result from Roche lobe
overflow of an AGB star. This evolutionary scenario may also explain the
formation of several wide WD+MS binaries discovered via self-lensing.

###### keywords:

binaries: spectroscopic – white dwarfs – stars: AGB and post-AGB – stars:
evolution

††pubyear: 2015††pagerange: Wide post-common envelope binaries containing
ultramassive white dwarfs: evidence for efficient envelope ejection in massive
AGB stars–B

##  1 Introduction

Common envelope evolution (CEE) is a major unsolved problem in binary
evolution. CEE is the outcome of dynamically unstable mass transfer (MT),
generally from a more massive donor to a less massive accretor. During CEE,
both stars orbit inside a shared envelope, spiraling inward on a dynamical or
thermal timescale. In some cases, the orbital energy liberated during this
inspiral is sufficient to eject the shared envelope, leaving behind a close
binary in which at least one component has lost most of its envelope. If
envelope ejection is not successful, the final outcome of CEE is a stellar
merger. Modeling of CEE is a key uncertainty in our understanding of the
formation of a wide variety of binary systems, including cataclysmic variables
(e.g. Paczynski, 1976; Meyer & Meyer-Hofmeister, 1979; Willems & Kolb, 2004),
X-ray binaries (e.g. Kalogera & Webbink, 1998), Type Ia supernovae (e.g.
Webbink, 1984; Meng & Podsiadlowski, 2017), binary neutron stars (e.g.
Bhattacharya & van den Heuvel, 1991), and binary black holes (Belczynski et
al., 2016; Marchant et al., 2021).

CEE is a dynamical process often involving an enormous range of physical and
temporal scales. Because detailed, end-to-end calculations of CEE are
currently infeasible (e.g. Ivanova et al., 2013) – and because it is often
necessary to model the evolution of large numbers of binaries to understand
the possible formation pathways of a single observed system – binary
population synthesis (BPS; e.g. Hurley et al., 2002) codes are often used to
model the evolution of millions of binaries, making it possible to explore a
broad parameter space at the expense of physical realism. These codes make use
of simplified models of CEE based on energy or angular momentum conservation.
In the most widely used formalism, it is assumed that a fixed fraction
(“$\alpha$”) of the liberated orbital energy goes in ejecting the envelope.
This fraction, and the binding energy of the envelope, then sets the post-CEE
orbital separation (e.g. Livio & Soker, 1988; Tout et al., 1997; De Marco et
al., 2011). Other energy sources, such as the photons released during
recombination when the envelope expands, are often modeled as reducing the
envelope’s binding energy.

CEE models have historically been calibrated by comparing the binary
populations they predict to observed post-common envelope binaries (PCEBs).
The most abundantly observed PCEBs contain white dwarfs or hot subdwarfs in
tight orbits ($P_{\rm orb}\lesssim 1$ d) with a main-sequence (MS) star,
having been produced when the MS star spiraled through the envelope of a red
giant and ultimately ejected it. BPS models have succeeded in explaining some
broad population properties of these binaries when using the
$\alpha$-formalism (Han et al., 2002, 2003; Camacho et al., 2014). Such
modeling makes it possible to empirically constrain $\alpha$, and several
calculations have found that the observations can best be reproduced by models
assuming $\alpha\approx 0.3$, meaning that $\approx 30\%$ of the orbital
energy liberated during inspiral goes into ejecting the envelope (Zorotovic et
al., 2010; Davis et al., 2012; Toonen & Nelemans, 2013; Camacho et al., 2014;
Zorotovic & Schreiber, 2022; Scherbak & Fuller, 2023).

At least one PCEB is known with a (relatively) wide orbit. That system, IK
Peg, has a period of 22 days and hosts an unusually massive WD, with mass
$M\approx 1.2\,M_{\odot}$ (Wonnacott et al., 1993). The system’s wide orbit
($a\approx 0.2$ au) means that less orbital energy was liberated during the MS
star’s inspiral than in typical PCEBs with $a\approx 0.01$ au. Zorotovic et
al. (2010) found that in the $\alpha-$formalism, IK Peg’s orbit can only be
explained when additional sources of energy besides orbital inspiral are taken
into account. Besides IK Peg, several wide WD+MS binaries have been discovered
via self-lensing, with orbital periods ranging from a few months to a few
years (Kruse & Agol, 2014; Kawahara et al., 2018). While it is not clear
whether these systems formed via CEE (see Section 5.1), they are also
candidates for being wide PCEBs and would require additional energy sources
(and/or high $\alpha$ values) to explain (Zorotovic et al., 2014). Energy
released by H and He recombination in the expanding envelope of the WD
progenitor (e.g. Webbink, 2008; Ivanova et al., 2015; Ivanova, 2018) is a
prime suspect for supplying the additional energy.

Most PCEBs studied to date were identified via their composite spectra and RV
variability detectable with low-resolution spectra (Rebassa-Mansergas et al.,
2007; Rebassa-Mansergas et al., 2017; Lagos et al., 2022). This leads to
strong selection effects in favor of PCEBs containing low-mass MS stars (which
are less likely to outshine the WD) in tight orbits (where RV shifts are
larger). The recent 3rd data release of the Gaia mission (DR3; Gaia
Collaboration et al., 2023a) contains orbital solutions for more than $10^{5}$
astrometric binaries, and for more than $10^{5}$ single-lined spectroscopic
binaries identified from medium-resolution spectra (Gaia Collaboration et al.,
2023b). This dataset provides a new opportunity to search for PCEBs with wider
orbits and more massive MS companions.

In this paper, we present five binaries in relatively wide orbits containing
solar-type main sequence stars and probable ultramassive WD candiates. Section
2 describes our identification of wide PCEB candidates from the Gaia DR3
catalog. Section 3 describes follow-up spectroscopic observations to obtain
radial velocities (RVs), spectral analysis to calculate metallicities, and
fitting to the broadband spectral energy distributions to constrain stellar
parameters of the MS stars. In Section 4, we fit the RVs to infer orbital
solutions. In particular, we measure a mass function which, when combined with
the luminous star mass, yields a minimum mass for the compact object. We also
discuss alternative possibilities for the nature of the unseen companions. In
Section 5, we compare our systems to other known PCEBs. Section 6 describes
models of the massive WD progenitors and constraints on CEE. In Section 7, we
briefly describe an alternative CE formalism, the occurrence rate of close and
wide PCEBs, and selection biases in past surveys. Finally, in Section 8, we
summarize our main results and conclude.

##  2 Discovery

The five objects studied in this paper (Table 1) were discovered in the course
of a broader search for compact objects with single-lined spectroscopic
(“SB1”) or astrometric + spectroscopic (“AstroSpectroSB1”) solutions in the
Gaia DR3 non-single star (NSS) catalog (Gaia Collaboration et al., 2023b). We
selected promising candidates for further follow-up based on their mass
functions, CMD positions, and Gaia quality flags. In brief, we targeted
sources whose CMD positions suggested a single luminous source and whose Gaia
mass functions implied a companion mass near the Chandrasekhar limit. For
objects with SB1 solutions, we prioritized those for which Bashi et al. (2022)
reported a robustness “score” above 0.5, corresponding roughly to an expected
20% contamination rate with spurious solutions.

Our spectroscopic follow-up revealed some sources to have spurious Gaia
orbital solutions and others to be double- or triple-lined binaries. Here, we
focus on 5 promising sources that are single-lined and whose Gaia-reported
orbits were validated by our follow-up. All 5 of these sources turned out to
have near-circular orbits, but eccentricity did not enter our initial
selection, and we did not find any similar (single-lined, high mass function)
targets with comparable periods and higher eccentricities. Our full search
will be described in future work.

Four targets have spectroscopic SB1 solutions, but no astrometric binary
solution. We suspect this is a result of the stringent cuts on astrometric
signal-to-noise ratio applied to astrometric solutions with short periods (see
Gaia Collaboration et al., 2023b). For these objects, the inclination is
unknown, and only a minimum companion mass can be inferred. One object,
J1314+3818, has a joint astrometric and spectroscopic (“AstroSpectroSB1”)
solution, meaning that its inclination is constrained. This object was
identified as a likely MS + compact object binary by Shahaf et al. (2023) on
the basis of its large astrometric mass ratio function. Another object in our
sample, J2034-5037, was previously identified by Jayasinghe et al. (2023) as a
candidate neutron star + MS binary.

Name | Gaia DR3 ID | RA [deg] | Dec [deg] |  $G$ [mag] | RUWE | $\varpi[\rm mas]$
---|---|---|---|---|---|---
J2117+0332 | 2692960678029100800 | 319.34490 | 3.54044 | 12.47 | 1.27 | 1.96 $\pm$ 0.02
J1111+5515 | 843829411442724864 | 167.80947 | 55.26410 | 10.61 | 1.47 | 3.24 $\pm$ 0.02
J1314+3818 | 1522897482203494784 | 198.51734 | 38.30119 | 11.05 | - | 12.45 $\pm$ 0.02
J2034-5037 | 6475655404885617920 | 308.60840 | -50.62557 | 12.37 | 2.94 | 3.23 $\pm$ 0.04
J0107-2827 | 5033197892724532736 | 16.98021 | -28.46128 | 12.27 | 1.74 | 2.14 $\pm$ 0.02
Table 1: Basic information from Gaia DR3 of the five objects found in this
work. The format for the name of each object is ‘J’ for J2000 followed by the
coordinates of the right ascension (RA) in hours and minutes, and declination
(Dec) in degrees and minutes. $G$ is the G-band mean magnitude, RUWE is the
Renormalised Unit Weight Error, and $\varpi$ is the parallax.

##  3 Follow-up

Here, we describe the follow-up spectra that we obtained, the process of
measuring metallicities from these spectra, and our constraints on the MS
stars’ parameters from their spectral energy distributions. A log of our
observations and measured RVs can be found in Appendix A.

###  3.1 FEROS

We obtained 59 spectra with the Fiberfed Extended Range Optical Spectrograph
(FEROS; Kaufer et al., 1999) on the 2.2 m ESO/MPG telescope at La Silla
Observatory (programs P109.A-9001, P110.A-9014, and P111.A-9003). Some
observations used $2\times 2$ binning to reduce readout noise at the expense
of spectral resolution; the rest used $1\times 1$ binning. The resulting
spectra have resolution $R\approx 40,000$ ($2\times 2$ binning) and $R\approx
50,000$ ($1\times 1$ binning). Exposure times ranged from 1200 to 1800
seconds.

We reduced the data using the CERES pipeline (Brahm et al., 2017), which
performs bias-subtraction, flat fielding, wavelength calibration, and optimal
extraction. The pipeline measures and corrects for small shifts in the
wavelength solution during the course a night via simultaneous observations of
a ThAr lamp obtained with a second fiber. We first calculate RVs by cross-
correlating a synthetic template spectrum with each order individually and
then report the mean RV across 15 orders with wavelengths between 4500 and
6700 Å. We calculate the uncertainty on this mean RV from the dispersion
between orders; i.e., $\sigma_{{\rm RV}}={\rm std}\left({\rm
RVs}\right)/\sqrt{15}$. We used a Kurucz spectral template from the BOSZ grid
(Bohlin et al., 2017) matched to the effective temperature of each star.

###  3.2 TRES

We obtained 34 spectra using the Tillinghast Reflector Echelle Spectrograph
(TRES; Fűrész, 2008) mounted on the 1.5 m Tillinghast Reflector telescope at
the Fred Lawrence Whipple Observatory (FLWO) atop Mount Hopkins, Arizona. TRES
is a fibrefed echelle spectrograph with a wavelength range of 390–910 nm and
spectral resolution $R\approx 44,000$. Exposure times ranged from 1800 to 3600
seconds. We extracted the spectra as described in Buchhave et al. (2010).

As with the FEROS data, we measured RVs by cross-correlating the normalized
spectra from each of 31 orders with a Kurucz spectrum template, and we
estimate RV uncertainties from the dispersion between RVs measured from
different orders; i.e., $\sigma_{{\rm RV}}={\rm std}\left({\rm
RVs}\right)/\sqrt{31}$.

###  3.3 MIKE

We observed J2117+0332 and J2034-5037 with the Magellan Inamori Kyocera
Echelle (MIKE) spectrograph on the Magellan 2 telescope at Las Campanas
Observatory (Bernstein et al., 2003). We used the 0.7” slit with an exposure
of 600s. This yielded a spectral resolution $R\sim 35,000$ and typical SNR of
$\sim$ 35 and 16 per pixel on the red and blue side respectively. The total
wavelength coverage was $\sim 3330-9680$ Å (though we only used spectra below
6850 Å to avoid telluric line contamination when computing the metallicity).
The spectra were reduced with the MIKE Pipeline using CarPy (Kelson et al.,
2000; Kelson, 2003). We flux-calibrated the spectra using a standard star and
merged the orders into a single spectrum, weighting by inverse variance in the
overlap regions.

###  3.4 Metallicities

Measuring metallicities of the MS stars is important for constraining their
masses and ages.

####  3.4.1 SPC

We fit the TRES spectra using the Stellar Parameter Classification (SPC) tool
(Buchhave et al., 2012). This code cross-correlates a grid of synthetic
spectra with each observed spectrum in the wavelength range of 5050 to 5360 Å,
centered on the Mg I b triplet. It then fits the peaks of the CCF with a three
dimensional third order polynomial to return best-fit values of effective
temperature $T_{\rm eff}$, surface gravity $\log g$, and metallicity [M/H]
that may lie in between the spacings of the grid.

As described in the supplementary material of Buchhave et al. (2012), given
systematic uncertainties in the synthetic stellar spectra, error floors on the
derived [M/H] and $T_{\rm eff}$ values are $\sim$ 0.08 dex and $\sim$ 50 K,
respectively (See also Furlan et al., 2018). We report these floors in in
Table 2.

####  3.4.2 BACCHUS

For the MIKE and FEROS spectra, we used the Brussels Automatic Code for
Characterizing High accUracy Spectra (BACCHUS; Masseron et al., 2016; Hayes et
al., 2022). This code performs 1D LTE spectral synthesis to determine stellar
parameters. It carries out normalization by linearly fitting the continuum 30Å
around a line. It then uses several methods to compare each line of the
observed spectrum to that of synthetic spectra to calculate an abundance. The
effective temperature, surface gravity, and microturbulence are estimated by
determining values that result in null trends between the inferred abundances
of a given element against the excitation potential, ionization potential, and
equivalent widths, respectively. The metallicity [Fe/H] is the mean Fe
abundance calculated over lines in the VALD atomic linelist (Piskunov et al.,
1995; Ryabchikova et al., 2015) with a wavelength coverage of 4200 to 9200Å.
We assume the detailed abundance pattern traces solar values. The errors
reported by BACCHUS represent the scatter in the implied abundances between
the different lines and methods of abundance calculations but do not take into
account other systematic uncertainties (Hayes et al., 2022).

####  3.4.3 Gaia XP

We also compare the values measured with SPC and BACCHUS to those calculated
by Andrae et al. (2023) using the Gaia XP very low-resolution spectra. These
authors derive $T_{\rm eff}$, log(g), and [M/H] for 175 million stars with XP
spectra published in DR3. Although the spectra from which these parameters are
derived have low resolution, Andrae et al. (2023) demonstrated that their
reported metallicities are accurate to within better than 0.1 dex for bright
and nearby stars like our targets with temperatures within the range of our
sample.

####  3.4.4 Results

The metallicities and effective temperatures obtained from spectral fitting
are summarized in Table 2. The metallicities range from $-0.15$ to 0.20 dex.
For J2117+0332, we see that the metallicities from SPC and BACCHUS are in
agreement. The Gaia XP metallicites are not used in our analysis in the
following sections but provide a useful comparison point. Most of our [M/H]
measurements are consistent with the Gaia XP measurements from Andrae et al.
(2023) within $1\sigma$. The good agreement between the three metallicities
shows that XP metallicities are likely sufficiently accurate for analysis of
larger samples in cases where high resolution follow-up would be prohibitively
expensive. We also add a column for the best-fit [Fe/H] values from our SED
fitting (Section 3.6), which uses the SPC and BACCHUS metallicities as a
prior.

[Fe/H] $T_{\rm{eff}}$ [K] name SPC BACCHUS Gaia XP SED SPC BACCHUS Gaia XP SED
J2117+0332 -0.24 $\pm$ 0.08 -0.284 $\pm$ 0.18 -0.380 -0.22 $\pm$ 0.06 6029
$\pm$ 50 6152 +/- 79 6111.0 6226 $\pm$ 19 J1111+5515 -0.15 $\pm$ 0.08 - -0.172
-0.17 $\pm$ 0.06 5987 $\pm$ 50 - 6006.3 6190 $\pm$ 22 J1314+3818 -0.39 $\pm$
0.08 - -0.291 -0.34 $\pm$ 0.05 4707 $\pm$ 50 - 4700.2 4684 $\pm$ 13 J2034-5037
- -0.346 $\pm$ 0.078 -0.352 -0.19 $\pm$ 0.06 - 5789 $\pm$ 17 5758.8 5856 $\pm$
20 J0107-2827 - 0.198 $\pm$ 0.127 0.244 0.04 $\pm$ 0.07 - 5524 $\pm$ 51 5330.4
5387 $\pm$ 21

Table 2: Comparison of the metallicities and $T_{\rm{eff}}$ obtained from
various methods.

###  3.5 Light curves

We retrieved observed light curves for our objects from the All-Sky Automated
Survey for Supernovae (ASAS-SN; Shappee et al., 2014; Kochanek et al., 2017).
We used the $V$ band data, for which the number of photometric points ranged
from 1609 to 3194 across the five objects. The typical uncertainty in
normalized flux is $\sim 0.01$. To search for periodic variability, we
computed Lomb-Scargle periodograms of these light curves (Lomb, 1976; Scargle,
1982; Astropy Collaboration et al., 2022). We did not find any significant
periodicities beyond the lunar cycle and sidereal day. This allows us to rule
out periodic variability with amplitude greater than the strongest noise
peaks, which have amplitude $\sim 0.002-0.003$ for all objects except
J1314+3813, where they have amplitude $\sim 0.006$.

###  3.6 SED fitting

We constructed broadband spectral energy distributions (SEDs) of our targets
using synthetic ugriz SDSS photometry calculated from Gaia XP spectra (Gaia
Collaboration et al., 2022) (with the exception of J2117+0332 where actual
SDSS photometry was available and used instead; Padmanabhan et al. 2008),
2MASS JHK photometry (Skrutskie et al., 2006), and WISE $W_{1}W_{2}W_{3}$
photometry (Wright et al., 2010). We obtained $E(B-V)$ for each object using
the Lallement et al. (2022) 3D dust map for declinations below -30∘ and the
Bayestar2019 3D dust map (Green et al., 2019) for declinations above -30∘.
These are given in Table 3. We assume a Cardelli et al. (1989) extinction law
with $R_{V}=3.1$. The Bayestar2019 map provides $E(g-r)$ which is
approximately equal to $E(B-V)$ (Schlafly & Finkbeiner, 2011), while the
Lallement et al. (2022) map provides the extinction $A_{0}$ at 550 nm which we
take to be $A_{V}$. As all our objects are relatively nearby with
$E(B-V)<0.05$, the uncertainties in these extinction values do not dominate
the uncertainties in the final fitted parameters.

We do not attempt to account for flux contributions from the WD companions,
which must be very small (given their high masses) and faint in the optical.
We justify this assumption in Appendix B, where we show that even very hot WDs
with $T_{\rm eff}=60,000$ K would not significantly contribute to the
photometry of all but one of our targets. For the one exception, we find that
a WD with $T_{\rm eff}\gtrsim$ 30,000 K could contribute to the $u$-band
photometry, so we conservatively excluded the $u-$band measurement from our
fit.

We fit the SEDs using MINEsweeper (Cargile et al., 2020), a code designed for
joint modeling of stellar photometry and spectra. We only use the code’s
photometric modeling capabilities but place a prior on the present-day surface
metallicity from spectroscopy. The free parameters to be fit are each star’s
parallax, mass $M_{\star}$, initial metallicity [Fe/H]init, and Equivalent
Evolutionary Phase (EEP, a monotonic function of age; see Dotter 2016). From
each set of parameters, MINEsweeper generates a predicted SED and photometry
in specified filters using neural network interpolation. We use emcee, a
Python Markov chain Monte Carlo (MCMC) sampler (Foreman-Mackey et al., 2013),
to sample from posterior. Constraints from fitting each source’s SED are
listed in Table 3.

We note that MINEsweeper constrains the initial metallicity, which is not
identical to the present-day surface value measured from spectroscopy. For our
targets, the difference between initial and present-day surface metallicity is
a result of atomic diffusion, where heavier elements settle out of the
atmosphere over time (Dotter et al., 2017). The present-day surface
metallicity [Fe/H] is predicted by the isochrones given a set of $M_{\star}$,
[Fe/H]init, and EEP, so the spectroscopic metallicities found in Section 3.4
are used to add a Gaussian constraint on [Fe/H] to the likelihood.

Putting everything together, the final likelihood function is:

$\begin{split}\ln L=&-\frac{1}{2}\sum_{i}\frac{\left(\rm
mag_{pred,i}-mag_{obs,i}\right)^{2}}{\sigma^{2}_{\rm
mag,i}}-\frac{1}{2}\frac{\left(\rm[Fe/H]_{pred}-[Fe/H]_{obs}\right)^{2}}{\sigma^{2}_{\rm[Fe/H]}}\\\
&-\frac{1}{2}\frac{\left(\rm\varpi_{pred}-\varpi_{obs}\right)^{2}}{\sigma^{2}_{\varpi}},\end{split}$
(1)

where “mag" stands for apparent magnitudes and the summation is over the
appropriate photometric filters for each object, $\sigma_{X}$ is the error on
the observed value of some quantity $X$, and [Fe/H] is the present-day surface
metallicity. We set a floor on $\sigma_{\rm mag}$ of 0.02 dex (given possible
calibration issues) to avoid underestimating the errors.

We report the medians of the marginalized posterior distributions for each
parameter in Table 3. $M_{\star}$, EEP, and [Fe/H]init are the parameters
directly fitted for by MINEsweeper, while [Fe/H], $T_{\rm eff}$, and $R$ are
calculated from the isochrones corresponding to the fitted parameters. The
reported errors are described in Section 3.6.2. We also list constraints on
[Fe/H] and $T_{\rm eff}$ for comparison with the values measured from
spectroscopy (Table 2).

Figure 1 shows MIST isochrones corresponding to the stellar parameters of 100
random posterior samples (gray). The cyan lines show the best-fit parameters
and the red point marks the present inferred parameters of the MS stars. Two
systems host stars that have evolved slightly evolved off the MS. This is
likely to be the result of selection bias, as evolved stars are brighter and
thus over-represented in magnitude-limited samples. In addition, we assumed
that stars were on the MS when estimating their masses in our initial
selection of targets for follow-up. These initial estimates were moderately
overestimated for evolved stars, leading to overestimated companion masses.
Since we targeted massive companions – and massive companions are
intrinsically rare – we expect evolved MS stars to be preferentially selected.

The observed and predicted SEDs are shown in Figure 2. The model SEDs plotted
were generated using pytstelllibs111https://mfouesneau.github.io/pystellibs/
with the best-fit parameters as inputs. We have checked that these models give
roughly consistent photometry to that predicted with MINEsweeper which does
not itself return a continuous SED. The residuals of the photometry typically
lie within 0.1 mag. Where available, GALEX NUV points are shown in Figure 2,
but these were not used in the fitting as they could have significant
contamination from the WD companion. For J1314+3818, we also found that WDs
with $T_{\rm eff}\gtrsim 30000$ K would significantly contribute to the SDSS
$u$-band photometry (See Appendix B), so this point was also excluded in the
fitting.

####  3.6.1 A wide tertiary

One system, J0107-2827, has a resolved tertiary separated by a distance of
2.21 arcseconds, corresponding to a projected physical separation of 1033 AU.
The consistency in the parallaxes and proper motions of the two sources make
it highly likely that they are in fact physically bound, as opposed to a
chance alignment (e.g. El-Badry et al., 2021). While the source is resolved by
Gaia in the $G$ band, the XP, 2MASS, and WISE photometry of the two sources
are likely all unresolved, so we model its SED as a sum of two luminous stars.
Assuming the tertiary is on the main sequence, its $G-$band absolute magnitude
of $\sim 7.3$ (calculated using the reported apparent magnitude of 15.7)
corresponds to a mass of approximately $0.67\,M_{\odot}$. We assume solar
metallicity (consistent with the initial metallicity we infer for the primary)
and an age of $\sim 6$ Gyr. Using these parameters, we generated photometry
for this third star (gray line in Figure 2) which we added to the model
primary (black). This sum (magenta) was fit to the observations.

####  3.6.2 Effect of potentially underestimated parallax errors

Our fitting also leaves the parallax, $\varpi$, free, allowing us to propagate
parallaxes uncertainties through to the stellar parameters. From Table 1, we
see that the Renormalised Unit Weight Errors (RUWEs) from Gaia DR3 for several
of our objects are above 1.4, which may indicate that the reported parallax
uncertainties are underestimated (Lindegren, 2018) as a result of orbital
motion, which is not accounted for in the Gaia single-star astrometric model.

To estimate more realistic parallax uncertainties, we carry out the following
analysis. We select sources from the Gaia NSS catalog with Orbital or
AstroSpectroSB1 solutions. In addition to the “single-star model” parallaxes
reported for these solutions in the gaia_source table, these sources have
improved parallaxes from astrometric solutions that account for wobble induced
by their binarity (Gaia Collaboration et al., 2023b; Halbwachs et al., 2023).
From these, we select those with phot_g_mean_mag < 13 and RUWE values
comparable to our targets. We then calculate the standard deviation of the
difference between the the parallaxes reported from single-star solutions
(gaia_source) and binary solution (NSS catalog). We found a standard deviation
of 0.104 mas for 13,132 sources with RUWE between 1.4 and 2, and a standard
deviation of 0.181 mas for 18,351 sources with RUWE between 2 and 3. The
maximum RUWE for our objects is 2.94 (Table 1).

Based on these values, we re-run the SED fitting with increased parallax
uncertainties of 0.2 mas for J2034-5037 (with a RUWE > 2) and 0.1 mas for the
remaining four objects. We find no significant changes to the best-fit values
of the parameters but a general increase in the uncertainties (i.e. the
standard derivations of the parameters from the posterior). We report these
inflated uncertainties in Table 3. Note that in Section 4.2, we also set an
uncertainity floor of $\pm 0.05\,M_{\odot}$ on $M_{\star}$ to obtain
conservative errors on the inferred masses of the WDs.

Name | $E(B-V)$ |  $M_{\star}$ [$M_{\odot}$] | EEP | [Fe/H]init |  $\varpi$ [mas] | [Fe/H] |  $T_{\rm eff}$ [K] |  $R$ [$R_{\odot}$]
---|---|---|---|---|---|---|---|---
J2117+0332 | 0.045 | 1.11 $\pm$ 0.03 | 383.26 $\pm$ 16.40 | -0.09 $\pm$ 0.05 | 1.96 $\pm$ 0.10 | -0.23 $\pm$ 0.05 | 6297 $\pm$ 20 | 1.25 $\pm$ 0.06
J1111+5515 | 0.009 | 1.15 $\pm$ 0.02 | 444.90 $\pm$ 2.92 | -0.09 $\pm$ 0.06 | 3.24 $\pm$ 0.10 | -0.17 $\pm$ 0.06 | 6169 $\pm$ 22 | 1.78 $\pm$ 0.06
J1314+3818 | 0.000 | 0.71 $\pm$ 0.01 | 276.03 $\pm$ 35.93 | -0.33 $\pm$ 0.04 | 12.45 $\pm$ 0.10 | -0.34 $\pm$ 0.05 | 4670 $\pm$ 13 | 0.71 $\pm$ 0.01
J2034-5037 | 0.024 | 0.96 $\pm$ 0.02 | 321.80 $\pm$ 52.27 | -0.16 $\pm$ 0.06 | 3.24 $\pm$ 0.17 | -0.18 $\pm$ 0.06 | 5857 $\pm$ 19 | 0.89 $\pm$ 0.05
J0107-2827 | 0.027 | 0.97 $\pm$ 0.03 | 459.26 $\pm$ 1.39 | 0.08 $\pm$ 0.07 | 2.14 $\pm$ 0.08 | 0.01 $\pm$ 0.06 | 5325 $\pm$ 19 | 1.71 $\pm$ 0.07
Table 3: Best-fit parameters from SED fitting. We have also added the
extinction $E(B-V)$ for all objects – these were taken from 3D dust maps and
have uncertainties of $\sim 0.02$ mag. $M_{\star}$ is the mass of the luminous
star, EEP is the Equivalent Evolutionary Phase (related to its age),
[Fe/H]init is the initial metallicity of the star, and $\varpi$ is the
parallax. For $M_{\star}$, we set an uncertainity floor of $\pm
0.05\,M_{\odot}$ when calculating WD masses in Section 4.2. The other
parameters are inferred from the isochrone corresponding to the fitted
parameters, where [Fe/H] is the present-day metallicity, $T_{\rm eff}$ is the
effective temperature, and $R$ is the radius. Figure 1: MIST isochrones for
all of our objects. In each plot, gray lines are the isochrones given the
masses and metallicities of 100 randomly chosen posterior samples from the
photometric fitting, while the cyan line is that of the best-fit parameters.
The red point marks the present location of the object and the error bars are
the standard deviation in the radii and temperatures of the posteriors at the
corresponding ages. Figure 2: Model SEDs for all five objects using the best-
fit parameters from fitting the photometry using MINEsweeper, compared with
the observations. J0107-2827 has a visible third companion whose SED is shown
in gray. The magenta line is the sum of the two luminous stars (the black line
corresponding to the more luminous star lies close to, but under, the magenta
line). Note that while GALEX points are plotted, these were not used in the
photometric fitting to avoid possible contamination from the WD companions.

##  4 Orbital fits

We fit the FEROS and TRES RVs with a Keplerian model using emcee. The free
parameters of the fit are the orbital period $P_{\rm orb}$, periastron time
$T_{p}$, eccentricity $e$, RV semi-amplitude $K_{\star}$, argument of
periastron $\omega$, and center-of-mass RV $\gamma$. In the case of
J2117+0332, where we have RVs from two different instruments, we also fit for
an RV offset between the two instruments as an additional parameter. We set
broad, uniform priors on all parameters. The likelihood function is defined
as:

$\ln L=-\frac{1}{2}\sum_{i}\frac{\left({\rm RV}_{\rm
pred}\left(t_{i}\right)-{\rm RV}_{i}\right)^{2}}{\sigma^{2}_{{\rm RV},i}}$ (2)

where ${\rm RV}_{\rm pred}\left(t_{i}\right)$ and ${\rm RV}_{i}$ are the
predicted and measured RVs at times $t_{i}$, and $\sigma^{2}_{{\rm RV},i}$ are
the errors in the measurements.

The best-fit RV curve for each object is shown in Figure 3. Best-fit values
for $P_{\rm orb}$, $e$, and $K_{\star}$, along with the implied mass functions
$f_{m}$ given these parameters, are reported in Table 4. For comparison, we
also list the mass functions calculated using the same parameters from the
Gaia DR3 SB1 solutions.

We find that all of the systems have a small but non-zero eccentricity. To
confirm that these are significant, we also fit the RVs using a model with
eccentricity and $\omega$ fixed to zero. The residuals from the two models are
plotted on the second and third panels for each object in Figure 3. We see
that the residuals from the model with $e=0$ are obviously larger than those
from the model that fits for $e$, with the possible exception of J2117+0332
(which has $e=0.0007\pm 0.0002$) where the difference is more subtle.

Since eccentricity cannot be negative, observational eccentricities will
result in a positive eccentricity bias for orbits that are nearly circular
(e.g. Hara et al., 2019). For J2117+0332, we generate simulated RVs with the
orbital parameters of the $e=0$ fit at the JDs of our observations, adding to
them Gaussian noise with a standard deviation of $0.05\,\rm km\,s^{-1}$. We
then fit these RVs with a Keplerian model, which yields $e\sim 0.0003\pm
0.0002$. This is comparable to the uncertainty on $e$ we find with the
measured RVs, and smaller than the measured eccentricity. This experiment
provides additional support that the non-zero eccentricity measured for
J2117+0332 is real.

Figure 3: Results of the RV fitting. For each object, the top panel shows the best-fit RV curve over the observed points, the second panel shows the residuals of this fit, the third panel shows the residuals of the fit with eccentricity set to 0, and the final panel shows the implied WD mass as a function of inclination (Section 4.2). For all objects (except possibly J2117+0332), the residuals for the model which fits for the eccentricity are significantly smaller than those with eccentricity fixed to zero. For J1314+3818, we fit the RVs and Gaia astrometry simultaneously (Section 4.1), with the red cross indicating the best-fit $i$ and the corresponding $M_{\rm WD}\sim 1.324\,M_{\odot}$. Name |  $P_{\rm orb}$ [d] | $e$ |  $K_{\star}$ [km/s] |  $f_{m}$ [$M_{\odot}$] |  $f_{m,G}$ [$M_{\odot}$]
---|---|---|---|---|---
J2117+0332 | 17.9239 $\pm$ 0.0001 | 0.0007 $\pm$ 0.0002 | 57.215 $\pm$ 0.011 | 0.3478 $\pm$ 0.0002 | 0.4143 $\pm$ 0.0398
J1111+5515 | 32.1494 $\pm$ 0.0022 | 0.0217 $\pm$ 0.0003 | 49.435 $\pm$ 0.019 | 0.4021 $\pm$ 0.0004 | 0.3981 $\pm$ 0.0052
J1314+3818 | 45.5150 $\pm$ 0.0047 | 0.0503 $\pm$ 0.0003 | 48.468 $\pm$ 0.015 | 0.5349 $\pm$ 0.0005 | –
J2034-5037 | 46.1147 $\pm$ 0.0006 | 0.0079 $\pm$ 0.0002 | 47.299 $\pm$ 0.012 | 0.5056 $\pm$ 0.0004 | 0.6392 $\pm$ 0.0944
J0107-2827 | 49.0063 $\pm$ 0.0008 | 0.0901 $\pm$ 0.0005 | 43.370 $\pm$ 0.011 | 0.4092 $\pm$ 0.0003 | 0.4175 $\pm$ 0.0275
Table 4: Best-fit orbital parameters from the RV fitting. $f_{m}$ is the mass
function given the other three parameters and $f_{m,G}$ is the mass function
given the same parameters from the Gaia solution.

###  4.1 Joint astrometric + RV fitting

One of our targets, J1314+3818, has a Gaia AstroSpectroSB1 solution, meaning
that the Gaia astrometry and RVs were fit with a combined orbital model. This
model has 15 parameters: ra, dec, parallax, pmra, pmdec, a_thiele_innes,
b_thiele_innes, f_thiele_innes, g_thiele_innes, c_thiele_innes,
h_thiele_innes, center_of_mass_velocity, eccentricity, period, t_periastron
(Halbwachs et al., 2023; Gaia Collaboration et al., 2023b). The Thiele-Innes
elements A, B, F, and G describe the astrometric orbit of the photocenter and
are transformations of the Campbell elements. The C and H elements are
constrained from the Gaia RVs of the MS star. In the case of a dark companion,
the photocenter simply traces the MS star.

We fit our RVs and the Gaia constraints simultaneously, using the likelihood
described by El-Badry et al. (2023) which we briefly summarize here. Gaia
stores the correlation matrix of the parameters in a vector corr_vec, from
which, along with the errors on the parameters, we can construct a covariance
matrix. We then construct a log-likelihood function that is a sum of two
terms: one that compares the predicted astrometric parameters and all Thiele-
Innes coefficients to the Gaia constraints, and one that compares the measured
and predicted RVs (Equation 2).

The model we fit has 14 free parameters: ra, dec, parallax, pmra, pmdec,
eccentricity $e$, inclination $i$, angle of the ascending node $\Omega$,
argument of periastron $\omega$, periastron time $T_{p}$, center-of-mass RV
$\gamma$, companion mass $M_{\rm WD}$, and luminous star mass $M_{\star}$. For
$M_{\star}$, we set a Gaussian constraint based on its best fit value and
error obtained from the SED fitting (Table 3). The resulting parameters can be
found in Table 5 and the plots are shown in Figure 3. We find $i\approx
100.0^{\circ}$ (or $\approx 80^{\circ}$) and $M_{\rm WD}=1.306\pm
0.012\,M_{\odot}$.

| Astrometry Only | Astrometry + RV
---|---|---
A [deg] | 198.517 $\pm$ 0.023 | 198.532 $\pm$ 0.010
Dec [deg] | 38.299 $\pm$ 0.013 | 38.307 $\pm$ 0.012
$\pi$ [mas] | 12.446 $\pm$ 0.015 | 12.447 $\pm$ 0.015
PMRA [mas/yr] | 129.523 $\pm$ 0.013 | 129.526 $\pm$ 0.011
PMDEC [mas/yr] | -224.216 $\pm$ 0.013 | -224.217 $\pm$ 0.013
$P_{\rm orb}$ [days] | 45.516 $\pm$ 0.005 | 45.519 $\pm$ 0.000
$e$ | 0.046 $\pm$ 0.006 | 0.0504 $\pm$ 0.0003
$i$ [deg] | 99.834 $\pm$ 0.370 | 99.945 $\pm$ 0.350
$\Omega$ [deg] | 86.017 $\pm$ 0.339 | 85.997 $\pm$ 0.344
$\omega$ [deg] | 99.219 $\pm$ 7.147 | 93.220 $\pm$ 0.307
$t_{\rm peri}$ [days] | -10.996 $\pm$ 0.916 | -11.763 $\pm$ 0.047
$\gamma$ [km/s] | 2.445 $\pm$ 0.209 | 2.811 $\pm$ 0.010
$M_{\rm WD}$ [M⊙] | 1.325 $\pm$ 0.046 | 1.324 $\pm$ 0.037
$M_{\star}$ [M⊙] | 0.713 $\pm$ 0.050 | 0.712 $\pm$ 0.049
Table 5: Best-fit parameters of J1314+3818 from fitting just the astrometric
solution, and when combined with the RVs.

###  4.2 White dwarf masses

From the parameters of the RV fitting, we can calculate the mass function,
$f_{m}$, which provides a constraint on the mass of the unseen (WD) companion
$M_{\rm WD}$:

$f_{m}=\frac{M_{\rm WD}^{3}\sin^{3}i}{\left(M_{\star}+M_{\rm
WD}\right)^{2}}=\frac{P_{\rm orb}K_{\star}^{3}}{2\pi
G}\left(1-e^{2}\right)^{3/2}$ (3)

where $M_{\star}$ is the mass of the luminous star, which we constrained by
fitting the SED (Section 3.6). With just the RVs, the inclination $i$ is not
constrained, meaning that for most of our objects, we can only place a lower
limit on the WD mass which occurs when $i=90^{\circ}$ (i.e. when the orbit is
“edge-on” to our line of sight).

The implied $M_{\rm WD}$ as a function of the inclination is plotted on the
lower panels of Figure 3. We shade the regions for $M_{\star}$ values $\pm
0.05\,M_{\odot}$ above and below the best-fit value from the SED fitting
(Section 3.6.2). The minimum masses $M_{\rm WD,min}$ range from
$1.244^{+0.027}_{-0.027}$ to $1.418^{+0.033}_{-0.033}\,M_{\odot}$. Given the
uncertainties, these values are all consistent with masses close to, but just
below, the Chandrasekhar limit of $\sim 1.4\,M_{\odot}$.

For J1314+3818, we obtain an inclination constraint from astrometry (Section
4.1) and thus a precise value for $M_{\rm WD}$ of $1.323\pm 0.008\,M_{\odot}$,
as opposed to just a lower limit. This point is shown as a red cross on the
plot of the $M_{\rm WD}(i)$ for this object in Figure 3.

The inferred WD masses are summarized in Table 6. At the time of writing,
these WDs are among the most massive WDs known (e.g. Hermes et al., 2013; Curd
et al., 2017; Cognard et al., 2017; Hollands et al., 2020; Caiazzo et al.,
2021; Miller et al., 2022), if they are indeed WDs (see Section 4.3). We note
that most other ultramassive WD candidates have mass estimates that depend on
WD cooling models and mass-radius relations, which are uncertain for the most
massive WDs (e.g. Camisassa et al., 2019; Schwab, 2021). Meanwhile, our
measurements (and similarly, those of Cognard et al., 2017) provide fairly
robust constraints on the mass with minimal assumptions about the WD itself
(though there is still dependence on the stellar models used to infer the mass
of the main-sequence companions).

name | MWD,min [M⊙]
---|---
J2117+0332 | 1.244${}^{+0.027}_{-0.027}$
J1111+5515 | 1.367${}^{+0.028}_{-0.029}$
J1314+3818 | 1.324$\pm$0.037 *
J2034-5037 | 1.418${}^{+0.033}_{-0.033}$
J0107-2827 | 1.271${}^{+0.030}_{-0.031}$
Table 6: Minimum white dwarf masses using orbital parameters from the RV
fitting and $M_{\star}$ from the SED fitting. J1314+3818 is marked in asterisk
(*) to indicate that this value is not a lower bound, but the actual value
given the inclination constraint from its astrometric solution (Table 5).

###  4.3 Nature of the unseen companions

Here we discuss whether the unseen companions could be objects other than WDs.

####  4.3.1 MS binaries or triples

A MS companion with a mass of $\sim 1.3\,M_{\odot}$ would dominate the SEDs of
all the objects in our sample. In this case we would see two sets of lines in
the spectra and changes in the composite line profiles with orbital phase.
Since the spectra of our targets are all well-fit by single-star models, we
can definitively rule out a single MS companion.

A different possibility is that these systems are hierarchical triples
consisting of a close inner binary of two $\sim 0.65\,M_{\odot}$ MS stars
orbiting the primary (e.g. van den Heuvel & Tauris, 2020). Together, the two
would be dimmer than a single $1.3\,M_{\odot}$ MS star. We can estimate the
contribution of such a binary to the overall SED in a similar way we do for
the case of a WD in Appendix B. We once again use pytstelllibs to generate an
SED but for a $0.65\,M_{\odot}$ star on the MS. This mass roughly corresponds
to a K7V star with a radius and temperature of $0.63R_{\odot}$ and $4100$ K
respectively (e.g. Pecaut & Mamajek, 2013). We can then calculate the ratio of
the flux from the two stars to that of the single star which was fitted for in
Section 3.6. At 550 nm, the fractional flux contribution of such an inner
triple would be 4.9, 2.6, 66.9, 13.24, and 5.2%, respectively, for J2117+0332,
J111+5515, J1314+3818, J2034-5037, and J0107-2827. In the infrared, at 3
$\mathrm{\mu m}$, the contribution would be larger, ranging from $\sim 13$ to
$56\%$ for four objects, with the exception of J1314+3818 where the inner
binary would outshine the tertiary. Therefore, for J1314+3818, a hierarchical
triple model is untenable. For the other objects, it is less obvious as there
is a relatively small contribution at optical wavelengths meaning that any
colour difference or spectral contribution is likely not enough to be
distinguished from a single source. We also note that a WD+WD or WD+MS inner
binary would similarly be dim in the optical and difficult to detect, but
forming these in close orbits would be challenging from an evolutionary
standpoint (given the size of the WD progenitor) and thus we do not consider
them further.

We also tested whether an inner binary’s presence could be inferred from the
SED fit. For each system, we constructed a “mock triple” SED by adding
synthetic photometry for the inner binary to synthetic photometry for the
best-fit single star model. We then fit this composite SED with a single-star
model using MINEsweeper and check whether the residuals of the fit worsens
significantly compared to those in Section 3.6. We find that while the median
residual does worsen slightly (at most by a factor of a few), the residuals of
most photometric points still remain within $\lesssim 0.1$ mag. As expected,
the exception is J1314+3818, where the residuals reach 0.2 mag. We conclude
that the worst-case inner binaries could escape detection via a poor SED fit
in all systems except J1314+3818.

We next consider the possible periods of hypothetical inner binaries. There is
a maximum period set by dynamical stability considerations: $P_{\rm
out}/P_{\rm in}\gtrsim 5$ where $P_{\rm out}$ and $P_{\rm in}$ are the outer
and inner orbital periods respectively (we are also taking $e\sim 0$ and the
ratio of the mass of the outer star to that of the inner binary $\sim 1$;
Mardling & Aarseth, 2001; Tokovinin, 2014). Given that our objects have
$P_{\rm orb}\sim 30$ days, this implies $P_{\rm in}\lesssim 6$ days. As for
the minimum period, if the orbit of the inner binary is sufficiently tight, we
may detect ellipsoidal variability due to tidal distortion of the inner
components. We use the code PHysics Of Eclipsing BinariEs (PHOEBE; Prša &
Zwitter, 2005) to generate synthetic light curves of an inner binary of two K
dwarfs for a range of periods from $\sim 0.25$ to 1 days. The amplitude of the
ellipsoidal variability decreases with increasing period. We then add a
fraction of the signal of these synthetic light curves to the observed light
curves (described in Section 3.5). Given that an inner binary contributes
$\lesssim 10$ % of the total light (from above), we set this fraction to 0.1.
We then generate periodograms of these light curves (Astropy Collaboration et
al., 2022) and see whether or not we would be able to detect variability on
half the inner binary’s period. We find that with only $\sim 10$ % of the
light coming from the inner binary, ellipsoidal variability can only be
distinguished from the noise for $P_{\rm inner}\lesssim 0.3$ days. Thus, the
range of possible inner period is $\sim 0.3$ to 6 days.

In summary, with the available observations, we cannot rule out a tight MS
binary in four out of the five systems. However, we emphasize that there are
very few hierarchical triple systems that have outer orbital periods below
$\sim$ 1000 days (this is not a selection effect, see Tokovinin, 2014), while
all our systems have $P_{\rm orb}<50$ days. The few known compact hierarchical
triples that have been found all have significantly more eccentric outer
orbits than our objects, with values ranging from about 0.2 to 0.6. The only
exceptions known are two triples in which the outer tertiary is a giant with a
large radius, which likely circularized their orbits (Rappaport et al., 2022,
2023). All five of our systems have eccentricities close to zero and host MS
stars in orbits that would not circularize in a Hubble time. This fact is
easily understood if the companions are WDs – the binaries would have been
(partially) circularized when the WD progenitor was a red giant. If the
companions were tight MS binaries, there would be no reason to expect circular
outer orbits, and it would be very improbable for all 5 systems to have
$e<0.1$ by chance.

####  4.3.2 Neutron stars

As we report minimum masses that are very close to the Chandrasekhar limit
(Table 6), we also consider the possibility that the companions are neutron
stars (NSs). However, NSs are expected to be born with natal kicks which drive
their orbits to be eccentric (e.g. Hills, 1983; Colpi & Wasserman, 2002;
Podsiadlowski et al., 2005). Thus, we must consider formation mechanisms which
can explain the near zero eccentricities of our objects (Table 4).

In the case of no natal kick and spherically symmetric mass loss forming the
NS (e.g. Blaauw, 1961), the eccentricity acquired (taking an initial
eccentricity of zero) is given by

$e=\frac{\Delta m}{m_{c}+m_{2}}$ (4)

where $\Delta m$ is the mass lost, $m_{c}$ is the remaining core/NS mass, and
$m_{2}$ is the mass of the companion (Hills, 1983). In the case where a
$8\,M_{\odot}$ star explodes by a core-collapse supernova (SN) to form a
$1.3\,M_{\odot}$ NS around a $1\,M_{\odot}$ MS companion, we see that $e>1$
and the system will be unbound. In reality, the massive progenitor may lose a
significant amount of mass through winds or binary interactins prior to the
explosion in which case the eccentricity will be smaller and may allow the
binary to survive, though likely in an eccentric orbit. Even in the case of an
ultra-stripped SN explosion with $\sim 0.3\,M_{\odot}$ of ejecta (De et al.,
2018; Yao et al., 2020), we expect $e\sim 0.1$, significantly larger than the
majority of our systems (though even lower ejecta masses are possible; Tauris
et al. 2013, 2015). Moreover, if the SN is asymmetric (in its ejecta or
neutrino emission), a strong kick can be imparted on the NS, which will likely
result in large eccentricities, if it does not unbind the NS (e.g. Fryer et
al., 1999; Tauris et al., 1999). Thus, a NS formed in this way is unlikely to
be found in very circular orbits like our targets. Alternatively, a NS may
form from a massive WD accreting up to the Chandrasekhar limit through
accretion-induced collapse (for a recent review on this topic, see Wang & Liu,
2020). Here, the ejecta masses are expected to be significantly smaller,
though quite uncertain, with values ranging from $1\times 10^{-3}$ to
$0.05\,M_{\odot}$ (Darbha et al., 2010; Fryer et al., 1999). Such ejecta
masses could correspond to $e\lesssim 0.02$ which are consistent with the
eccentricities of some of our objects. However, it is difficult to explain how
the progenitor WD would have accreted the necessary mass to begin with. Our
systems all contain MS star companions which do not have strong winds, so
there should be no significant wind accretion. Moreover, our objects are in
orbits that are too wide for there to have been MT from the MS star through
RLOF. Thus, while accretion-induced collapse could produce NSs in circular
orbits, it struggles to do so in our systems where there are no obvious MT
mechanisms.

Therefore, we conclude that these alternative scenarios are improbable (though
we emphasize that they are possible) and we proceed under the assumption that
the unseen companions are WDs.

##  5 Comparison to other binary populations

Here we compare the properties of our targets to other related classes of
binaries, including other WD+MS PCEBs and WD + millisecond pulsar (MSP)
binaries.

###  5.1 Literature PCEBs

The Sloan Digital Sky Survey (SDSS; Abazajian et al., 2009) detected large
numbers of close WD+MS binaries. Rebassa-Mansergas et al. (2007) identified 37
new PCEBs from the SDSS PCEB survey (described in Section 7.3), and combining
this with 25 that were previously known, Zorotovic et al. (2010) compiled a
total of 62 PCEB systems. In addition, the “white dwarf binary pathways
survey” also identified several PCEBs with AFGK companions using TESS light
curves (diamond markers; Hernandez et al., 2021, 2022a, 2022b). In Figure 4,
we plot the minimum orbital separation $a_{\rm peri}$ against $M_{\rm WD}$ for
these literature PCEBs, the five objects from our sample, and several self-
lensing WD+MS binaries discovered by the Kepler survey (SLBs + KOI-3278;
Kawahara et al., 2018; Kruse & Agol, 2014). With the exception of J1314+3818
where the precise value of $M_{\rm WD}$ was obtained using astrometry (Section
4.1), we have plotted $M_{\rm WD,min}$ for our objects which is indicated with
arrows. The colours of the points represent the mass of the luminous (MS)
companion, $M_{\star}$. The gray dashed line comes from the $P_{\rm
orb}-M_{\rm WD}$ relation derived in Rappaport et al. (1995) for stable MT
(with a spread in orbital period of a factor of $\sim 2.4$), where $P_{\rm
orb}$ has been converted to separation assuming a $1\,M_{\odot}$ MS star and
circular ($e=0$) orbit (assuming instead a $0.1\,M_{\odot}$ M dwarf star only
shifts this downwards by a small amount). The fact that all of these objects
lie below this relation means that they are unlikely to have formed through
stable RLOF, with the possible exception of the SLBs, which are not too far
below the relation.

The blue line in Figure 4 shows the maximum radius of the WD progenitor, if it
evolved in isolation. Binaries located below this line must have interacted at
some point in their evolution. This approximate relation was obtained by first
calculating the progenitor mass $M_{\rm init}$ using the WD Initial-Final Mass
Relation (IFMR) derived in Williams et al. (2009): $M_{\rm
final}=0.339+0.129M_{\rm init}$, then generating MIST evolutionary tracks
(Dotter, 2016; Choi et al., 2016) to identify the maximum radius reached by a
star with a given $M_{\rm init}$. Note that this is a conservative limit since
the RLOF would begin before the giant touches the companion.

####  5.1.1 A population of PCEBs in wide orbits?

IK Peg was previously isolated in its region of the $P_{\rm orb}-M_{\rm WD}$
parameter space: being in a wider orbit with a period of 22 days and hosting a
more massive WD of $\sim 1.2\,M_{\odot}$ (Wonnacott et al., 1993) than the
vast majority of SDSS PCEBs. Our five targets fall in the same region as IK
Peg. Their current orbits are far too tight for the binaries to have escaped
interaction when the WD progenitors were red giants or AGB stars, strongly
suggesting that these objects are indeed PCEBs. As we show in Section 6, the
current orbits can only be understood as an outcome of CEE if additional
energy sources (besides liberated orbital energy) helped unbind the common
envelope.

The SLBs occupy a different isolated region in this space, with normal WD
masses but at separations even wider than our systems and IK Peg. Like our
targets, they contain solar-type MS stars, which are more massive than the M
dwarfs in the SDSS PCEB sample. Kruse & Agol (2014) initially interpreted
KOI-3278 as a “normal” PCEB, but Zorotovic et al. (2014) subsequently showed
that the system’s wide orbit requires an extra source of energy, beyond
orbital energy, to contribute to the CE ejection. The three SLBs identified by
Kawahara et al. (2018) have even wider orbits than KOI-3278. Those authors
interpreted the systems as having formed through stable MT, but Figure 4 shows
that SLB 2 and 3 fall well below the Rappaport et al. (1995) prediction.
Formation through stable MT thus seems tenable only if these WDs have
significantly overestimated masses.

We also distinguish the objects discovered by Hernandez et al. (2021, 2022a,
2022b) (plotted with diamond markers) from the SDSS PCEBs as they host higher-
mass mass ($\sim 1\,M_{\odot}$) MS stars. Compared to our objects, these have
shorter orbital periods ($P_{\rm orb}\sim 1-2$ days) and can therefore be
explained with just the liberated orbital energy, without needing to invoke
additional sources. This tells us that having an intermediate-mass MS star as
a companion does not necessarily lead to a wide post-CE orbit.

We note that while Figure 4 may suggest that there are two distinct groups of
PCEBs in wide orbits with different WD masses (self-lensing binaries vs. IK-
Peg analogs), this is not necessarily the case. We remind the reader that here
we specifically targeted very massive companions which would at least partly
explain why we did not find any that are less massive. A search for more
objects located in currently sparse regions on the plot would be useful.

Figure 4: Periastron separation, $a_{\rm peri}$, vs. WD mass, $M_{\rm WD}$,
for a sample of PCEBs from Zorotovic et al. (2011, circle markers) and the
five objects from this work (stars; the arrows indicate lower limits). IK Peg
is distinguished from the other known PCEBs with a triangle marker as it lies
very close to our objects. We also plot the PCEBs from the “white dwarf binary
pathways survey" (diamond markers; Hernandez et al., 2021, 2022a, 2022b).
Finally, we plot self-lensing binaries (SLBs) discovered by Kawahara et al.
(2018) as well as KOI-3278 (Kruse & Agol, 2014, also a SLB), which were all
detected by Kepler. The colours of the points represents the masses of the
luminous MS companions. The dashed gray shows the prediction for stable MT
from Rappaport et al. (1995) and the blue line indicates the maximum radius
reached by the WD progenitor for a range of WD masses.

###  5.2 Eccentricities and comparison to MSP binaries

On Figure 5, we plot the eccentricity $e$ against $P_{\rm orb}$ and compare
our objects to MSP + WD binaries. We plot objects primarily from the Australia
Telescope National Facility (ATNF) catalogue (Manchester et al., 2005) (taking
the version analyzed by Hui et al. 2018) and distinguish those with minimum WD
companion masses above and below $0.45\,M_{\odot}$ (approximate upper limit to
the mass of a He WD). We also plot the theoretical relation derived by Phinney
(1992) for MSP pulsars with He WDs formed through stable MT.

MSPs form through “recycling”, where an old NS is spun up to short periods by
the transfer of mass and angular momentum from a companion (e.g. Alpar et al.,
1982; Radhakrishnan & Srinivasan, 1982; Bhattacharya & van den Heuvel, 1991).
Tides are almost expected to circularize MSP + WD binaries, but a very low
orbital eccentricity remains because convection in the WD progenitor produces
a time-varying quadrupole moment, leading to perturbations and eccentricity
excursions that are larger in longer-period systems, which hosted larger
giants (Phinney, 1992). To date, the period–eccentricity relation has mainly
been tested with MSP+WD binaries, because their eccentricities can be easily
measured with high precision. However, a similar process should operate in
MS+WD binaries, if MT occurs over a long enough period for tidal
circularization to occur.

Figure 5 shows that in general, the MSP binaries with more massive (CO/ONeMg)
WDs tend to have higher eccentricities at fixed period than those with low
mass (He) WDs. The standard interpretation is that the systems with He WDs
formed via stable MT, while those with more massive WDs formed through CEE
(e.g. van den Heuvel, 1994; Tauris et al., 2012). Although NS + red giant
orbits are expected to be circularized prior to the onset of MT, eccentricity
is produced during the dynamical plunge-in phase of CEE (e.g. Ivanova et al.,
2013), and there is likely insufficient time for this eccentricity to be fully
damped between the end of CEE and the formation of the WD (e.g. Glanz &
Perets, 2021).

The objects in our sample have periods and eccentricities similar to the
longest-period MSP + CO/ONeMg binaries, perhaps pointing to a similar
formation history. In particular, the MSP binaries J1727-2946 (Lorimer et al.,
2015) and J2045+3633 (Berezina et al., 2017) (circled in green) are located
very close to our objects in Figure 5. Both of these systems contain mildly
recycled pulsars (with spin periods of 27 and 32 ms) and massive WDs ($M_{\rm
WD}>0.8\,M_{\odot}$), and they have eccentricities of 0.045 and 0.017 at orbit
periods of 40 and 32 days. A common envelope origin has also been proposed for
J2045+3633 (McKee et al., 2020). There are also stable MT scenarios that have
been proposed to explain the formation of MSP + CO/ONeMg WD binaries (e.g.
Tauris et al., 2000), which would likely predict lower eccentricities. We
refer readers to Tauris (2011) for a concise overview of this topic.

It is worth mentioning that several systems with low-mass WDs in anomalously
eccentric orbits have been discovered, called eccentric MSPs (eMSPs; e.g.
Bailes, 2010). These points are circled in red in Figure 5. Interestingly,
these eMSPs occupy a narrow range in orbital periods with eccentricities that
are comparable to the two eccentric MSP + CO WD binaries (as pointed out by
Berezina et al., 2017). This is unexpected, as binaries with He WDs are
thought to have formed through distinct evolutionary paths, involving long
periods of stable MT in low mass X-ray binaries (Tauris, 2011). There have
been multiple proposed mechanisms to explain the large eccentricities, some of
which may be applicable to one or more of the systems just discussed (MSP + He
WD, MSP + CO WD, MS + WD PCEB). These include MSPs being directly formed from
the accretion-induced collapse of a super-Chandrasekhar mass ONeMg WD (Freire
& Tauris, 2014), interaction with a circumbinary disk formed from material
lost from the WD during H shell flashes (Antoniadis, 2014), or a circumbinary
disk formed from the ejected envelope after a CE phase (Dermine et al., 2013).

Figure 5: Eccentricity vs. orbital period. We plot the sample of MSP + WD
binaries mainly from the ATNF catalogue (Manchester et al. 2005, with a few
others, compiled by Hui et al. 2018), differentiating between those with
minimum WD mass above and below $0.45\,M_{\odot}$ (orange triangles and blue
circles, respectively). We show our objects with magenta stars. The gray line
is the theoretical relation for MSP + He WD binaries formed through stable MT
derived by Phinney (1992). The points circled in red are eMSPs (e.g., see
Stovall et al., 2019, for a description of the five plotted here). We also
circle in green two binaries with massive CO WDs with large eccentricities
(Lorimer et al., 2015; Berezina et al., 2017).

##  6 Feasibility of formation through common envelope evolution

To test whether our targets could have formed via CEE, we ran MESA binary
evolution models (Paxton et al., 2011; Paxton et al., 2013, 2015, 2018) of
progenitor stars to the WDs in our systems, evolving intermediate-mass stars
up to the asymptotic giant branch (AGB).

In the $\alpha$-formalism, the CEE ends when the loss of orbital energy from
the spiral-in exceeds the binding energy of the envelope $E_{\rm bind}$,
resulting in the envelope being ejected:

$E_{\rm bind}=\alpha_{\rm CE}\left(-\frac{GM_{\rm
WD}M_{\star}}{2a_{f}}+\frac{GM_{i}M_{\star}}{2a_{i}}\right)$ (5)

where $M_{i}$ is the mass of the WD progenitor, $a_{i}$ is the initial
separation at the onset of CE, $a_{f}$ is the separation at the end of the CE
phase, and $\alpha_{\rm CE}$ is the fraction of the liberated orbital energy
that goes into unbinding the envelope. Thus, $E_{\rm bind}$ determines the
final separation of the system for a given initial separation.

In the simplest case, $E_{\rm bind}$ is just the gravitational binding energy
of the envelope. However, previous works using BPS have found that in this
case, no values of $\alpha_{\rm CE}<1$ can reproduce the relatively wide
orbits of IK Peg and KOI-3278 (Davis et al., 2010; Zorotovic et al., 2010,
2014; Parsons et al., 2023), which have comparable separations to our systems.
Moreover, it is clear that additional energy exists within the stellar
envelope, which can potentially help to unbind it. Here, we consider the
inclusion of internal energy which is defined in MESA as the sum of thermal
and recombination energy (Paxton et al., 2018). This makes the binding energy
less negative (possibly even positive), corresponding to an envelope that is
less bound. Recombination of H and He can occur if there is some process (in
this case, the binary interaction through the CE phase) that causes the
envelope to expand and cool down, which will release energy (e.g. Ivanova et
al., 2013).

From initial-final mass relations of WDs (e.g. Williams et al., 2009; Cummings
et al., 2018; El-Badry et al., 2018), we expect the MS progenitor masses of
our ultra-massive WDs to be in the range of $6-9\,M_{\odot}$. Thus, we run
MESA models of a $7\,M_{\odot}$ star, following its evolution up to the AGB.
Our inlists are based on those of Farmer et al. (2015) (same wind and mixing
prescriptions; For more details, we refer readers to Section 2 of Farmer et
al., 2015), although we have updated them for the more recent MESA version
r22.05.1. We only run non-rotating models.

The evolution of the $7\,M_{\odot}$ model, from pre-main sequence to
termination at the tip of the AGB, is shown on the HR diagram in the leftmost
panel of Figure 6. In the following, we consider progenitors on the red giant
branch (RGB), including the sub-giant branch (SGB), and the AGB which are
highlighted in blue and red on the diagram respectively. We define these
phases following the convention used in the MIST project (See Section 2.1 of
Dotter, 2016). We note that our model terminates before core carbon burning,
but as the envelope binding energy becomes positive before this point (see
below), our main conclusions should not be largely affected by this.

At each timestep, we calculate the binding energy of the envelope as a sum of
the gravitational ($E_{\rm grav}$) and internal ($E_{\rm int}$) components:

$\displaystyle E_{\rm bind}$ $\displaystyle=E_{\rm grav}+E_{\rm int}$ (6)
$\displaystyle=\int_{M_{\rm core}}^{M_{\rm tot}}-\frac{Gm}{r(m)}+U(m)dm$ (7)

where $E_{\rm grav}$ and $E_{\rm int}$ correspond to the first and second
parts of the integrand respectively, $m$ is the mass enclosed within a radius
$r$, and $U(m)$ is the internal energy per unit mass. The integral is taken
from the He core boundary to the surface of the star. We use the default
definition of the He core in MESA which is where the hydrogen mass fraction
$X_{H}<0.1$ and helium mass fraction $X_{\rm He}>0.1$. We have tested changing
these boundaries to 0.01 and find no significant change. Here, we assume that
all of the internal energy contributes to the envelope binding energy. We
discuss the consequences of relaxing this assumption in the following section.

Using equation 5, we can solve for $a_{f}$ for Roche lobe overflow occurring
at different points on the RGB and AGB. We calculate the initial separation
$a_{i}$ using the Eggleton formula (Eggleton, 1983), assuming the giant fills
its Roche lobe at the onset of the CEE. The mass of the WD remaining after the
envelope ejection, $M_{\rm WD}$, is assumed to be equal to the core mass of
the giant $M_{\rm core}$. We also assume $M_{\star}=1\,M_{\odot}$, which is
close to the median value for our objects (Table 3).

The predicted $a_{f}$ is shown for a range of $a_{i}$ in Figure 6. In the
central panel, the gray dashed line marks $a_{f}=0.15$ AU, which is slightly
smaller than the minimum $a_{\rm peri}$ of our objects at $\sim 0.18$ AU
(Figure 4). The red dashed line marks $a_{f}=0.01$ AU ($\sim 2R_{\odot}$),
below which a $\sim 1\,M_{\odot}$ MS star cannot fit inside the orbit and thus
a PCEB would not form (a merger, or perhaps stable MT of the MS onto the WD,
may occur instead). The orange line shows the case where we only consider the
gravitational binding energy ($E_{\rm bind}=E_{\rm grav}$) and set
$\alpha_{\rm CE}=1$ (i.e. 100% of the orbital energy loss goes into envelope
ejection). We see that it never crosses the $a_{f}=0.15$ AU mark, meaning that
even in this optimistic case, there is not enough energy to unbind the
envelope and produce orbits as wide as our observed systems. In the remaining
three cases, we include the internal energy ($E_{\rm bind}=E_{\rm grav}+E_{\rm
int}$, Equation 6) and let $\alpha_{\rm CE}=0.3$ (the “standard" value), 0.6,
and 0.9. We see that for each case, there is a region where $a_{f}$ exceeds
0.15 AU. On the right panel, we zoom into this region and find $a_{i}$ ranges
from $\sim$ 3.5 - 4.4 AU across all values of $\alpha_{\rm CE}$, with a
narrower range for lower values of $\alpha_{\rm CE}$.

In Figure 6, we exclude models in regions of parameter space where the
$\alpha$-formalism does not make clear predictions, namely, those in which
$E_{\rm bind}>0$ (the envelope is unbound when recombination energy is
included). This case is likely still relevant for producing wide PCEBs, but it
is unclear what the final separation should be: if the MS star does not
penetrate deep into the giant’s envelope, it is unlikely to trigger the
release of much recombination energy. For the model shown in Figure 6, these
conditions are realized late in the AGB evolution, corresponding to radii of
$\sim 520-888R_{\odot}$ and initial separations of $\sim 4.4-7.7$AU.

Figure 6: Left: HR diagram showing the evolution of a $7\,M_{\odot}$ star
starting from pre-MS to the AGB. The blue sections indicate what we refer to
as the RGB (but which also includes the SGB) and the red sections represents
the AGB. Center: Plots of the final separation $a_{f}$ (i.e. birth period at
the end of CEE) over a range of initial separations $a_{i}$. $a_{i}$ is taken
to be the orbital semi-major axis when the giant (WD progenitor) fills its
Roche lobe. We mark $a_{f}=$ 0.01 AU $\sim 2R_{\odot}$ (red dashed line) below
which the MS star would not fit in the orbit and a PCEB cannot form. The
orange dashed line is the case where only the gravitational binding energy is
considered and $\alpha_{\rm CE}=1$. We see that in this case, no values of
$a_{i}$ result in $a_{f}>0.15$ AU (gray dashed line) which is approximately
the minimum separation of our objects. The other three lines are the cases
where internal energy is added to the binding energy for $\alpha_{\rm
CE}=0.3$, 0.6, and 0.9. We see that these lines lie above the dashed line for
some range of $a_{i}$. Right: Zoom in on the region where $a_{f}>0.15$ AU. We
see that overall, $a_{i}\sim 3.5-4.4$ AU can result in the wide orbital
separations of our systems.

####  6.0.1 Efficiency of recombination energy

The calculation above consider two kinds of internal energy: thermal energy
and recombination energy. While thermal energy is commonly considered as an
“extra” energy source, it is closely related to the gravitational potential
energy through the virial theorem, and all the the thermal energy should be
included in the binding energy by default (e.g. Ivanova et al., 2013). The
inclusion of recombination energy is more uncertain.

Some works have argued that most of the energy released by recombination will
quickly be transported to the photosphere through radiation and/or convection
and then radiated away (Sabach et al., 2017; Grichener et al., 2018).
Meanwhile, Ivanova (2018) found that such energy transport is inefficient in
typical AGB stars and that recombination is in fact a significant source of
additional energy. Given this uncertainty, our previous assumption that all of
the internal energy contributes to unbinding the envelope may be too
optimistic.

We can assess the sensitivity of our results to this assumption by splitting
the internal energy into two components:

$E_{\rm bind}=E_{\rm grav}+\alpha_{\rm th}E_{\rm th}+\alpha_{\rm rec}E_{\rm
rec}$ (8)

where $E_{\rm th}$ and $E_{\rm rec}$ are the thermal and recombination
energies, and $\alpha_{\rm th}$ and $\alpha_{\rm rec}$ are the respective
efficiencies. As described above, we set $\alpha_{\rm th}=1$.

MESA does not provide a simple way to individually track the thermal and
recombination energies. We can, however, approximate the thermal energy using
the ideal gas law, which is a reasonable approximation in the envelopes of AGB
and RGB stars. In this case, the energy density is $(3/2)P/\rho$, where $P$ is
the pressure and $\rho$ is the mass density. We subtract this from the total
internal energy output by MESA to get an estimate of the recombination energy
and see what value of $\alpha_{\rm rec}$ would be required to produce wide
orbits.

As shown in Figure 7, we find that for a canonical value of $\alpha_{\rm
CE}=0.3$, we require $\alpha_{\rm rec}\gtrsim 0.5$, and for $\alpha_{\rm
CE}=1$, we require $\alpha_{\rm rec}\gtrsim 0.25$ for $a_{f}\gtrsim 0.15$ AU.
Thus, our models point towards a relatively large fraction of recombination
energy being needed to produce wide PCEBs. This is not unexpected, as
recombination energy dominates the internal energy of the envelopes of cool
stars. For typical stellar compositions, recombination energy dominates over
internal energy for temperatures below $\approx 2\times 10^{5}$ K (Ivanova et
al., 2013). This boundary lies deep in the envelope of our models on the AGB,
at roughly 20% of the stars’ radii.

If a large fraction of recombination energy actually escapes, our results
would imply that other sources of energy must be invoked to produce wide PCEBs
(e.g. jets from the accreting star; Sabach et al., 2017; Moreno Méndez et al.,
2017). We defer more detailed calculations and further discussion of this
topic to future work.

Figure 7: The same plot as in the central panel of Figure 6 but separating out
the thermal and recombination components of the internal energy $E_{\rm int}$.
As in Figure 6, we plot the case where only $E_{\rm grav}$ is considered with
$\alpha_{\rm CE}=1$ (orange dashed line). The pink and green lines show cases
where $\alpha_{\rm CE}=$ 0.3 and 1, and $\alpha_{\rm rec}=$ 0.5 and 0.25
respectively, where $a_{f}$ just exceeds the 0.15 AU mark.

###  6.1 Models for lower- and higher-mass giants

We performed similar calculations for a $1\,M_{\odot}$ red giant, which might
be a progenitor to a wide PCEB hosting a lower-mass ($0.5-0.6\,M_{\odot}$) WD.
We evolved a $1\,M_{\odot}$, solar-metallicity star using inlists from the
1M_pre_ms_to_wd calculation in the MESA test suite. The default wind
parameters in that calculation are set so that the mass loss is unusually
efficient on the AGB (Blocker_scaling_factor = 0.7). This speeds up the
calculation by preventing the star from evolving far up the AGB and
encountering thermal pulses, but it terminates the AGB phase unrealistically
early. We instead set Blocker_scaling_factor = 0.05 following Farmer et al.
(2015).

The same plots as in Figure 6 but for the $1\,M_{\odot}$ model are shown in
Figure 8. Even in the case where we consider only the gravitational binding
energy with $\alpha_{\rm CE}=1$, there is a range of initial separations
($\sim 1.5-4.5$ AU) for which $a_{f}>0.15$ AU. This only occurs on the
thermally-pulsating phase of the AGB (TP-AGB) where the envelope becomes very
loosely bound. Wide separations can also result for much smaller values of
$\alpha_{\rm CE}\sim 0.1$ if MT starts at the tip of the AGB. This suggests
that it is possible to produce PCEBs in wide orbits containing
$0.5-0.6\,M_{\odot}$ WDs (such as the self-lensing binaries) without the need
to invoke additional energy sources, but only if the MT begins on the TP-AGB.
It should, however, be kept in mind that the TP-AGB phase is a particularly
difficult phase to model (e.g. Marigo & Girardi, 2007; Girardi & Marigo,
2007), so this conclusion may depend on the adopted stellar models.

With the inclusion of internal energy, wide PCEB orbits are produced for a
broad range of initial separations. The full range of initial separations for
which the calculations predict $a_{f}$ > 0.15 AU is $\sim$ 0.9 to 4.5 AU. For
$a_{i}\sim 3-4.5$ AU, the envelope’s binding energy is positive when
recombination energy is included, so the $\alpha$ formalism does not
straightforwardly predict a final separation. For separations $a_{i}\lesssim
2.5$ AU, CEE would likely commence on the RGB, preventing the wide PCEB
outcome from being realized in practice. Overall, this calculation suggests
that efficient envelope ejection can similarly produce wide PCEBs with both
low and high mass WDs.

We also ran models of more massive stars, with initial masses of 12 and
$20\,M_{\odot}$, that become red supergiants and will leave behind neutron
stars or black holes. We find that the envelopes of these stars are
significantly more bound than that of the $7\,M_{\odot}$ super-AGB star. This
implies that it is difficult to form wide BH/NS + solar-type stars via CEE.
Similar conclusions have been reached by other studies (e.g. Kalogera, 1999;
Kiel & Hurley, 2006; Giacobbo & Mapelli, 2018; Fragos et al., 2019; El-Badry
et al., 2023).

Figure 8: Same as Figure 6 but for a $1\,M_{\odot}$ model. We omit the
$\alpha_{\rm CE}=0.6$ case to avoid over-crowding. Even in the case where just
$E_{\rm grav}$ is considered (orange), there is a range of initial separations
for which which wide PCEBs can be produced with $\alpha_{\rm CE}=1$. When
$E_{\rm int}$ is included, there is a wide range of initial separations for
which wide PCEBs can be produced, even with $\alpha_{\rm CE}<1$ (green and
purple). The models producing wide PCEBs are those in which MT begins during
the TP-AGB phase, when envelope is very weakly bound.

##  7 Discussion

###  7.1 Formation through stable MT?

Given that the binaries in our sample have wider orbits than traditional
PCEBs, it is natural to wonder whether they could have formed through stable
MT instead of CEE. We briefly discuss this possibility here.

The onset of dynamically unstable MT is determined by the donor star’s
adiabatic response to mass loss. In general, a system in which the donor is
more massive than the accretor will tend to be more unstable. The critical
mass ratio above which mass transfer is unstable, $q_{\rm crit}$, depends on
the stellar structure and evolutionary state of the donor. But in the case of
our systems with an ultramassive WD progenitor of mass $\sim 6-7\,M_{\odot}$
and an intermediate MS companion of $\sim 1\,M_{\odot}$, the mass ratio is
large enough to exceed even conservative values of $q_{\rm crit}\sim 3-4$
(donor/accretor; e.g. Hjellming & Webbink, 1987; Ge et al., 2010; Temmink et
al., 2023).

The non-zero eccentricities of our observed systems also points towards CEE,
which is expected to be less efficient at tidal circularization than stable
RLOF, and may actually drive small eccentricities during the plunge-in phase
or via torques from circumbinary material (e.g. Ivanova et al., 2013). These
arguments suggest that stable mass transfer is unlikely to have formed the
systems in our sample.

We would be remiss here to not mention the “$\gamma$”-formalism, another
commonly used prescription of CEE. This was originally invoked to model the
formation of double CO WD binaries, which were thought to require a widening
of the orbit after the first CE phase, which cannot occur in the
$\alpha$-formalism (Nelemans et al., 2000; Nelemans & Tout, 2005). The
parameter $\gamma$ can be understood as the ratio of the angular momentum lost
per mass of ejected material to the average angular momentum per unit mass of
the initial binary (Paczynski, 1976). While this formalism can produce the
wide orbits seen in our systems, it should be emphasized that it was designed
precisely for this purpose and does not fundamentally solve the issues
associated with energy conservation, which must still hold. It has also be
argued that the $\gamma$-formalism does not actually describe the CEE – the
result of unstable MT – but instead a phase of stable, non-conservative MT.
See Section 5 of Ivanova et al. (2013) for further discussion on this
formalism.

###  7.2 Relative frequency of wide and close PCEBs

The small number of wide PCEBs discovered so far raises the question of
whether they are intrinsically rarer than close PCEBs, or just more difficult
to detect. Here we describe the selection biases against wide PCEBs in
previous surveys. Given the complex and thus far poorly understood selection
function of the Gaia DR3 binary sample, we do not attempt to infer the space
density of wide PCEBs here. Instead, we compare the distances to various
samples of PCEBs as a rough diagnostic of their relative frequencies.

We cross-match the sample of literature PCEBs compiled by Zorotovic et al.
(2010, also shown in our Figure 4) to Gaia DR3 to obtain their parallaxes. We
find that the median distance to SDSS PCEBs within that sample is 328 pc,
which is significantly farther than the median distance of 108 pc for non-SDSS
PCEBs in the sample. This likely reflects the fact that most of the non-SDSS
objects were discovered serendipitously from all-sky studies of bright stars,
in many cases having been recognized as binaries via photometric variability.
In contrast, the SDSS objects were discovered spectroscopically from a parent
sample that is deep but only observed a small fraction of all stars. Our
targets have distances ranging from 80 to 510 pc, with a median of 308 pc
(Table 1). At 80 pc, J1314+3818 is nearer than any of the SDSS PCEBs. IK Peg,
another wide PCEB, is at 46 pc which is nearer than the majority of the close
PCEBs in the literature.

A particularly interesting case to consider is the binary G 203-47 (Delfosse
et al., 1999). That system contains a $\approx 0.27\,M_{\odot}$ MS star
orbiting a dark companion that is almost certainly a WD in a period of 14.7
days, similar to the wide PCEBs studied here. At a distance of only 7.5 pc, G
203-47 is one of the 10 nearest known WDs, and probably the nearest PCEB! It
is 3 times nearer – corresponding to a 27 times smaller search volume – than
the nearest short-period PCEB, RR Cae, but has been largely overlooked by
works attempting to constrain CE physics with PCEBs. While it is dangerous to
draw population-level conclusions from a single object, this strongly suggests
that wide PCEBs are quite common.

###  7.3 Comparison to other surveys

While the SDSS survey for PCEBs (Rebassa-Mansergas et al., 2007) was highly
effective at finding WD + M dwarf PCEBs in tight orbits ($\lesssim 1$ day), it
was biased against finding systems like the ones presented here. This is
because the sample was selected based on RV variations detected in low-
resolution BOSS spectra, which are more easily detected in close binaries with
short orbital periods. Furthermore, the SDSS PCEB survey identified candidates
by searching for sources with composite spectra in which contributions of both
the WD and the MS companion were detectable. This leads to a strong bias in
favor of low-mass (M dwarf) main-sequence companions.

The White Dwarf Binary Pathways Survey conducted a search for WD + AFGK PCEBs.
They first selected AFGK MS stars from the RAVE and LAMOST surveys, and then
cross-matched them to GALEX, identifying objects with UV excess as candidates
for having a WD companion (Parsons et al., 2016). From these WD+MS binary
candidates, they selected PCEB candidates as those binaries with RV variations
detectable in their low-resolution multi-epoch spectra, mainly from LAMOST
(Rebassa-Mansergas et al., 2017). This also leads to a strong bias in favor of
short periods.

The White Dwarf Binary Pathways Survey did find three binaries with orbital
periods of several weeks, but Lagos et al. (2022) concluded that they were
likely contaminants. The binaries in question had significant eccentricities
($e=0.266-0.497$; see Table 1 of Lagos et al., 2022), atypical of PCEBs. Based
on HST spectra and high contrast imaging, they concluded that at least two are
hierarchical triples in which the WD is a distant tertiary. Our objects would
likely not have been found by their search because they have negligible UV
excess.

##  8 Conclusions

We presented five post-common envelope binaries (PCEBs) containing ultra-
massive WD candidates and intermediate mass MS stars with long orbital periods
(18 - 49 days). These were discovered as part of a broader search for compact
object binaries from the Gaia DR3 NSS catalog. Previous surveys identified
PCEBs using a combination of RV variability, photometric variability, and UV
excess, which made them biased towards finding PCEBs with M dwarfs in close
orbits. Systems like the ones presented here pose a potential challenge in
simplified models of common envelope evolution (CEE) as their formation
requires loosely bound donor envelopes which can be quickly ejected, leaving
them in wide orbits with non-zero eccentricities. Our main findings are as
follows:

  1. 1.

Nature of the unseen companions: The companions are dark objects with masses
of $1.2-1.4\,M_{\odot}$ – more massive than the solar-type stars orbiting
them. The simplest explanation is that they are WDs. We consider two possible
alternatives: (1) a tight binary containing two $\sim 0.65\,M_{\odot}$ MS
stars. In the most pessimistic case, such an inner binary could escape
detection in 4 of our 5 targets. However, the near-circular orbits we observe
– which would be a natural consequence of tidal circularization if the
companions are WDs – are not expected in this hierarchical triple scenario. No
tight hierarchical triples with outer MS stars and circular outer orbits are
known, and very few triples are known with outer periods below 1000 days. (2)
a neutron star (NS). This is also unlikely due to the circular orbits of our
systems, as NSs are expected to be born with natal kicks that send them to
highly eccentric orbits. Given these considerations, we proceed under the
assumption that the unseen companions are WDs.

  2. 2.

WD masses: From RVs (Section 4), we measure orbital solutions and mass
functions. Combining this with the masses of the luminous components obtained
from SED fitting (Section 3.6), we obtain minimum WD masses. These range from
to $1.244^{+0.027}_{-0.027}$ to $1.418^{+0.033}_{-0.033}\,M_{\odot}$, all
consistent with masses just below the Chandrasekhar limit. One object,
J1314+3818, has a Gaia astrometric solution, which we fit simultaneously with
the RVs to constrain the inclination. For this object, we obtain a precise
mass of $1.324\pm 0.037\,M_{\odot}$. Assuming the dark companions are in fact
WDs, they are among the most massive WDs known.

  3. 3.

Comparison to other PCEBs: Our newly discovered systems have longer periods
and host more massive WDs and MS stars than most known PCEBs (Figure 4). The
only similar system previously known is IK Peg. However, it is important to
note that selection effects of most previous searches strongly favored short-
period PCEBs with low-mass MS stars.

  4. 4.

Comparison to MSP + WD binaries: We find that our objects have large
eccentricities relative to the bulk of the MSP+WD binaries (Figure 5).
However, they have similar eccentricities at fixed orbital period to MSP + CO
WD binaries, which likely also formed through CEE.

  5. 5.

Evolutionary models of the WD progenitors: We ran MESA models of a
$7\,M_{\odot}$ MS star up the RGB and AGB (Section 6), following its internal
energy (thermal + recombination) and calculating the expected final separation
according to the $\alpha$ formalism. We find that there is no point in the
evolutionary phase of the star where final separations comparable to those of
our objects ($\gtrsim 0.15$ AU) are predicted if internal energy does not aid
in unbinding the envelope. In the case where internal energy is included,
there is a range in initial separations ($\sim 3.5-4.5$ AU) where final
separations exceed $\sim 0.15$ AU. For initial separations wider than 4.5 AU,
the binding energy of the envelope is positive when CEE begins, such that a
range of (wide) final separations are plausible.

  6. 6.

Space density: We compare distances of literature PCEBs to those of our
objects and a few others in wider orbits. At $\sim 80$ pc, J1314+3818 is
nearer than any of the SDSS PCEBs (Zorotovic et al., 2010). The median
distance of objects in our sample is comparable to that of all literature
PCEBs. The nearest known PCEB, G 203-47, has a period of 15 days, much wider
than typical PCEBs in the literature with $P<1$ day. A detailed estimate of
the space density of wide PCEBs will have to wait for a better characterized
Gaia selection function, but these early discoveries suggest that it is
comparable to or larger than that of close PCEBs.

## Acknowledgements

We would like to thank Matthias Schreiber and Monica Zorotovic for providing
feedback that helped to improve this manuscript, and Thomas Tauris for
enlightening discussions. We also thank Thomas Masseron and Keith Hawkins for
assistance with BACCHUS. Finally, we thank Hans-Walter Rix and Eleonora Zari
for observing some of our objects under their FEROS programs.

N.Y. and K.E. were supported in part by NSF grant AST-2307232. This paper
includes data gathered with the 6.5 meter Magellan Telescopes located at Las
Campanas Observatory, Chile.

This work has made use of the VALD database, operated at Uppsala University,
the Institute of Astronomy RAS in Moscow, and the University of Vienna.

## Data Availability

The data underlying this article are available upon reasonable request to the
corresponding author.

## References

  * Abazajian et al. (2009) Abazajian K. N., et al., 2009, ApJS, 182, 543
  * Alpar et al. (1982) Alpar M. A., Cheng A. F., Ruderman M. A., Shaham J., 1982, Nature, 300, 728
  * Andrae et al. (2023) Andrae R., Rix H.-W., Chandra V., 2023, ApJS, 267, 8
  * Antoniadis (2014) Antoniadis J., 2014, ApJ, 797, L24
  * Astropy Collaboration et al. (2022) Astropy Collaboration et al., 2022, ApJ, 935, 167
  * Bailes (2010) Bailes M., 2010, New Astron. Rev., 54, 80
  * Bashi et al. (2022) Bashi D., Shahaf S., Mazeh T., Faigler S., Dong S., El-Badry K., Rix H. W., Jorissen A., 2022, MNRAS, 517, 3888
  * Belczynski et al. (2016) Belczynski K., Holz D. E., Bulik T., O’Shaughnessy R., 2016, Nature, 534, 512
  * Berezina et al. (2017) Berezina M., et al., 2017, MNRAS, 470, 4421
  * Bernstein et al. (2003) Bernstein R., Shectman S. A., Gunnels S. M., Mochnacki S., Athey A. E., 2003, in Iye M., Moorwood A. F. M., eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 4841, Instrument Design and Performance for Optical/Infrared Ground-based Telescopes. pp 1694–1704, doi:10.1117/12.461502
  * Bhattacharya & van den Heuvel (1991) Bhattacharya D., van den Heuvel E. P. J., 1991, Phys. Rep., 203, 1
  * Blaauw (1961) Blaauw A., 1961, Bull. Astron. Inst. Netherlands, 15, 265
  * Bohlin et al. (2017) Bohlin R. C., Mészáros S., Fleming S. W., Gordon K. D., Koekemoer A. M., Kovács J., 2017, AJ, 153, 234
  * Brahm et al. (2017) Brahm R., Jordán A., Espinoza N., 2017, PASP, 129, 034002
  * Buchhave et al. (2010) Buchhave L. A., et al., 2010, ApJ, 720, 1118
  * Buchhave et al. (2012) Buchhave L. A., et al., 2012, Nature, 486, 375
  * Caiazzo et al. (2021) Caiazzo I., et al., 2021, Nature, 595, 39
  * Camacho et al. (2014) Camacho J., Torres S., García-Berro E., Zorotovic M., Schreiber M. R., Rebassa-Mansergas A., Nebot Gómez-Morán A., Gänsicke B. T., 2014, A&A, 566, A86
  * Camisassa et al. (2019) Camisassa M. E., et al., 2019, A&A, 625, A87
  * Cardelli et al. (1989) Cardelli J. A., Clayton G. C., Mathis J. S., 1989, ApJ, 345, 245
  * Cargile et al. (2020) Cargile P. A., Conroy C., Johnson B. D., Ting Y.-S., Bonaca A., Dotter A., Speagle J. S., 2020, ApJ, 900, 28
  * Choi et al. (2016) Choi J., Dotter A., Conroy C., Cantiello M., Paxton B., Johnson B. D., 2016, ApJ, 823, 102
  * Cognard et al. (2017) Cognard I., et al., 2017, ApJ, 844, 128
  * Colpi & Wasserman (2002) Colpi M., Wasserman I., 2002, ApJ, 581, 1271
  * Cummings et al. (2018) Cummings J. D., Kalirai J. S., Tremblay P. E., Ramirez-Ruiz E., Choi J., 2018, ApJ, 866, 21
  * Curd et al. (2017) Curd B., Gianninas A., Bell K. J., Kilic M., Romero A. D., Allende Prieto C., Winget D. E., Winget K. I., 2017, MNRAS, 468, 239
  * Darbha et al. (2010) Darbha S., Metzger B. D., Quataert E., Kasen D., Nugent P., Thomas R., 2010, MNRAS, 409, 846
  * Davis et al. (2010) Davis P. J., Kolb U., Willems B., 2010, MNRAS, 403, 179
  * Davis et al. (2012) Davis P. J., Kolb U., Knigge C., 2012, MNRAS, 419, 287
  * De Marco et al. (2011) De Marco O., Passy J.-C., Moe M., Herwig F., Mac Low M.-M., Paxton B., 2011, MNRAS, 411, 2277
  * De et al. (2018) De K., et al., 2018, Science, 362, 201
  * Delfosse et al. (1999) Delfosse X., Forveille T., Beuzit J. L., Udry S., Mayor M., Perrier C., 1999, A&A, 344, 897
  * Dermine et al. (2013) Dermine T., Izzard R. G., Jorissen A., Van Winckel H., 2013, A&A, 551, A50
  * Dotter (2016) Dotter A., 2016, ApJS, 222, 8
  * Dotter et al. (2017) Dotter A., Conroy C., Cargile P., Asplund M., 2017, ApJ, 840, 99
  * Eggleton (1983) Eggleton P. P., 1983, ApJ, 268, 368
  * El-Badry et al. (2018) El-Badry K., Rix H.-W., Weisz D. R., 2018, ApJ, 860, L17
  * El-Badry et al. (2021) El-Badry K., Rix H.-W., Heintz T. M., 2021, MNRAS, 506, 2269
  * El-Badry et al. (2023) El-Badry K., et al., 2023, MNRAS, 518, 1057
  * Fűrész (2008) Fűrész G., 2008, PhD thesis, University of Szeged, Hungary
  * Farmer et al. (2015) Farmer R., Fields C. E., Timmes F. X., 2015, ApJ, 807, 184
  * Foreman-Mackey et al. (2013) Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306
  * Fragos et al. (2019) Fragos T., Andrews J. J., Ramirez-Ruiz E., Meynet G., Kalogera V., Taam R. E., Zezas A., 2019, ApJ, 883, L45
  * Freire & Tauris (2014) Freire P. C. C., Tauris T. M., 2014, MNRAS, 438, L86
  * Fryer et al. (1999) Fryer C., Benz W., Herant M., Colgate S. A., 1999, ApJ, 516, 892
  * Furlan et al. (2018) Furlan E., et al., 2018, ApJ, 861, 149
  * Gaia Collaboration et al. (2022) Gaia Collaboration et al., 2022, arXiv e-prints, p. arXiv:2206.06215
  * Gaia Collaboration et al. (2023a) Gaia Collaboration et al., 2023a, A&A, 674, A1
  * Gaia Collaboration et al. (2023b) Gaia Collaboration et al., 2023b, A&A, 674, A34
  * Ge et al. (2010) Ge H., Hjellming M. S., Webbink R. F., Chen X., Han Z., 2010, ApJ, 717, 724
  * Giacobbo & Mapelli (2018) Giacobbo N., Mapelli M., 2018, MNRAS, 480, 2011
  * Girardi & Marigo (2007) Girardi L., Marigo P., 2007, A&A, 462, 237
  * Glanz & Perets (2021) Glanz H., Perets H. B., 2021, MNRAS, 507, 2659
  * Green et al. (2019) Green G. M., Schlafly E., Zucker C., Speagle J. S., Finkbeiner D., 2019, ApJ, 887, 93
  * Grichener et al. (2018) Grichener A., Sabach E., Soker N., 2018, MNRAS, 478, 1818
  * Halbwachs et al. (2023) Halbwachs J.-L., et al., 2023, A&A, 674, A9
  * Han et al. (2002) Han Z., Podsiadlowski P., Maxted P. F. L., Marsh T. R., Ivanova N., 2002, MNRAS, 336, 449
  * Han et al. (2003) Han Z., Podsiadlowski P., Maxted P. F. L., Marsh T. R., 2003, MNRAS, 341, 669
  * Hara et al. (2019) Hara N. C., Boué G., Laskar J., Delisle J. B., Unger N., 2019, MNRAS, 489, 738
  * Hayes et al. (2022) Hayes C. R., et al., 2022, ApJS, 262, 34
  * Hermes et al. (2013) Hermes J. J., Kepler S. O., Castanheira B. G., Gianninas A., Winget D. E., Montgomery M. H., Brown W. R., Harrold S. T., 2013, ApJ, 771, L2
  * Hernandez et al. (2021) Hernandez M. S., et al., 2021, MNRAS, 501, 1677
  * Hernandez et al. (2022a) Hernandez M. S., et al., 2022a, MNRAS, 512, 1843
  * Hernandez et al. (2022b) Hernandez M. S., et al., 2022b, MNRAS, 517, 2867
  * Hills (1983) Hills J. G., 1983, ApJ, 267, 322
  * Hjellming & Webbink (1987) Hjellming M. S., Webbink R. F., 1987, ApJ, 318, 794
  * Hollands et al. (2020) Hollands M. A., et al., 2020, Nature Astronomy, 4, 663
  * Hui et al. (2018) Hui C. Y., Wu K., Han Q., Kong A. K. H., Tam P. H. T., 2018, ApJ, 864, 30
  * Hurley et al. (2002) Hurley J. R., Tout C. A., Pols O. R., 2002, MNRAS, 329, 897
  * Ivanova (2018) Ivanova N., 2018, ApJ, 858, L24
  * Ivanova et al. (2013) Ivanova N., et al., 2013, A&ARv, 21, 59
  * Ivanova et al. (2015) Ivanova N., Justham S., Podsiadlowski P., 2015, MNRAS, 447, 2181
  * Jayasinghe et al. (2023) Jayasinghe T., Rowan D. M., Thompson T. A., Kochanek C. S., Stanek K. Z., 2023, MNRAS, 521, 5927
  * Kalogera (1999) Kalogera V., 1999, ApJ, 521, 723
  * Kalogera & Webbink (1998) Kalogera V., Webbink R. F., 1998, ApJ, 493, 351
  * Kaufer et al. (1999) Kaufer A., Stahl O., Tubbesing S., Nørregaard P., Avila G., Francois P., Pasquini L., Pizzella A., 1999, The Messenger, 95, 8
  * Kawahara et al. (2018) Kawahara H., Masuda K., MacLeod M., Latham D. W., Bieryla A., Benomar O., 2018, AJ, 155, 144
  * Kelson (2003) Kelson D. D., 2003, PASP, 115, 688
  * Kelson et al. (2000) Kelson D. D., Illingworth G. D., van Dokkum P. G., Franx M., 2000, ApJ, 531, 159
  * Kiel & Hurley (2006) Kiel P. D., Hurley J. R., 2006, MNRAS, 369, 1152
  * Kochanek et al. (2017) Kochanek C. S., et al., 2017, PASP, 129, 104502
  * Koester (2010) Koester D., 2010, Mem. Soc. Astron. Italiana, 81, 921
  * Kruse & Agol (2014) Kruse E., Agol E., 2014, Science, 344, 275
  * Lagos et al. (2022) Lagos F., Schreiber M. R., Parsons S. G., Toloza O., Gänsicke B. T., Hernandez M. S., Schmidtobreick L., Belloni D., 2022, MNRAS, 512, 2625
  * Lallement et al. (2022) Lallement R., Vergely J. L., Babusiaux C., Cox N. L. J., 2022, A&A, 661, A147
  * Lindegren (2018) Lindegren L., 2018, Gaia Data Processing and Analysis Consortium (DPAC) technical note GAIA-C3-TN-LU-LL-124,
  * Livio & Soker (1988) Livio M., Soker N., 1988, ApJ, 329, 764
  * Lomb (1976) Lomb N. R., 1976, Ap&SS, 39, 447
  * Lorimer et al. (2015) Lorimer D. R., et al., 2015, MNRAS, 450, 2185
  * Manchester et al. (2005) Manchester R. N., Hobbs G. B., Teoh A., Hobbs M., 2005, AJ, 129, 1993
  * Marchant et al. (2021) Marchant P., Pappas K. M. W., Gallegos-Garcia M., Berry C. P. L., Taam R. E., Kalogera V., Podsiadlowski P., 2021, A&A, 650, A107
  * Mardling & Aarseth (2001) Mardling R. A., Aarseth S. J., 2001, MNRAS, 321, 398
  * Marigo & Girardi (2007) Marigo P., Girardi L., 2007, A&A, 469, 239
  * Masseron et al. (2016) Masseron T., Merle T., Hawkins K., 2016, BACCHUS: Brussels Automatic Code for Characterizing High accUracy Spectra, Astrophysics Source Code Library, record ascl:1605.004 (ascl:1605.004)
  * McKee et al. (2020) McKee J. W., et al., 2020, MNRAS, 499, 4082
  * Meng & Podsiadlowski (2017) Meng X., Podsiadlowski P., 2017, MNRAS, 469, 4763
  * Meyer & Meyer-Hofmeister (1979) Meyer F., Meyer-Hofmeister E., 1979, A&A, 78, 167
  * Miller et al. (2022) Miller D. R., Caiazzo I., Heyl J., Richer H. B., Tremblay P.-E., 2022, ApJ, 926, L24
  * Moreno Méndez et al. (2017) Moreno Méndez E., López-Cámara D., De Colle F., 2017, MNRAS, 470, 2929
  * Nelemans & Tout (2005) Nelemans G., Tout C. A., 2005, MNRAS, 356, 753
  * Nelemans et al. (2000) Nelemans G., Verbunt F., Yungelson L. R., Portegies Zwart S. F., 2000, A&A, 360, 1011
  * Paczynski (1976) Paczynski B., 1976, in Eggleton P., Mitton S., Whelan J., eds, Vol. 73, Structure and Evolution of Close Binary Systems. p. 75
  * Padmanabhan et al. (2008) Padmanabhan N., et al., 2008, ApJ, 674, 1217
  * Parsons et al. (2016) Parsons S. G., Rebassa-Mansergas A., Schreiber M. R., Gänsicke B. T., Zorotovic M., Ren J. J., 2016, MNRAS, 463, 2125
  * Parsons et al. (2023) Parsons S. G., et al., 2023, MNRAS, 518, 4579
  * Paxton et al. (2011) Paxton B., Bildsten L., Dotter A., Herwig F., Lesaffre P., Timmes F., 2011, ApJS, 192, 3
  * Paxton et al. (2013) Paxton B., et al., 2013, ApJS, 208, 4
  * Paxton et al. (2015) Paxton B., et al., 2015, ApJS, 220, 15
  * Paxton et al. (2018) Paxton B., et al., 2018, ApJS, 234, 34
  * Pecaut & Mamajek (2013) Pecaut M. J., Mamajek E. E., 2013, ApJS, 208, 9
  * Phinney (1992) Phinney E. S., 1992, Philosophical Transactions of the Royal Society of London Series A, 341, 39
  * Piskunov et al. (1995) Piskunov N. E., Kupka F., Ryabchikova T. A., Weiss W. W., Jeffery C. S., 1995, A&AS, 112, 525
  * Podsiadlowski et al. (2005) Podsiadlowski P., Pfahl E., Rappaport S., 2005, in Rasio F. A., Stairs I. H., eds, Astronomical Society of the Pacific Conference Series Vol. 328, Binary Radio Pulsars. p. 327
  * Prša & Zwitter (2005) Prša A., Zwitter T., 2005, ApJ, 628, 426
  * Radhakrishnan & Srinivasan (1982) Radhakrishnan V., Srinivasan G., 1982, Current Science, 51, 1096
  * Rappaport et al. (1995) Rappaport S., Podsiadlowski P., Joss P. C., Di Stefano R., Han Z., 1995, MNRAS, 273, 731
  * Rappaport et al. (2022) Rappaport S. A., et al., 2022, MNRAS, 513, 4341
  * Rappaport et al. (2023) Rappaport S. A., et al., 2023, MNRAS, 521, 558
  * Rebassa-Mansergas et al. (2007) Rebassa-Mansergas A., Gänsicke B. T., Rodríguez-Gil P., Schreiber M. R., Koester D., 2007, MNRAS, 382, 1377
  * Rebassa-Mansergas et al. (2017) Rebassa-Mansergas A., et al., 2017, MNRAS, 472, 4193
  * Ryabchikova et al. (2015) Ryabchikova T., Piskunov N., Kurucz R. L., Stempels H. C., Heiter U., Pakhomov Y., Barklem P. S., 2015, Phys. Scr., 90, 054005
  * Sabach et al. (2017) Sabach E., Hillel S., Schreier R., Soker N., 2017, MNRAS, 472, 4361
  * Scargle (1982) Scargle J. D., 1982, ApJ, 263, 835
  * Scherbak & Fuller (2023) Scherbak P., Fuller J., 2023, MNRAS, 518, 3966
  * Schlafly & Finkbeiner (2011) Schlafly E. F., Finkbeiner D. P., 2011, ApJ, 737, 103
  * Schwab (2021) Schwab J., 2021, ApJ, 916, 119
  * Shahaf et al. (2023) Shahaf S., Bashi D., Mazeh T., Faigler S., Arenou F., El-Badry K., Rix H. W., 2023, MNRAS, 518, 2991
  * Shappee et al. (2014) Shappee B. J., et al., 2014, ApJ, 788, 48
  * Skrutskie et al. (2006) Skrutskie M. F., et al., 2006, AJ, 131, 1163
  * Stovall et al. (2019) Stovall K., et al., 2019, ApJ, 870, 74
  * Tauris (2011) Tauris T. M., 2011, in Schmidtobreick L., Schreiber M. R., Tappert C., eds, Astronomical Society of the Pacific Conference Series Vol. 447, Evolution of Compact Binaries. p. 285 (arXiv:1106.0897), doi:10.48550/arXiv.1106.0897
  * Tauris et al. (1999) Tauris T. M., Fender R. P., van den Heuvel E. P. J., Johnston H. M., Wu K., 1999, MNRAS, 310, 1165
  * Tauris et al. (2000) Tauris T. M., van den Heuvel E. P. J., Savonije G. J., 2000, ApJ, 530, L93
  * Tauris et al. (2012) Tauris T. M., Langer N., Kramer M., 2012, MNRAS, 425, 1601
  * Tauris et al. (2013) Tauris T. M., Langer N., Moriya T. J., Podsiadlowski P., Yoon S. C., Blinnikov S. I., 2013, ApJ, 778, L23
  * Tauris et al. (2015) Tauris T. M., Langer N., Podsiadlowski P., 2015, MNRAS, 451, 2123
  * Temmink et al. (2023) Temmink K. D., Pols O. R., Justham S., Istrate A. G., Toonen S., 2023, A&A, 669, A45
  * Tokovinin (2014) Tokovinin A., 2014, AJ, 147, 87
  * Toonen & Nelemans (2013) Toonen S., Nelemans G., 2013, A&A, 557, A87
  * Tout et al. (1997) Tout C. A., Aarseth S. J., Pols O. R., Eggleton P. P., 1997, MNRAS, 291, 732
  * Tremblay & Bergeron (2009) Tremblay P. E., Bergeron P., 2009, ApJ, 696, 1755
  * Wang & Liu (2020) Wang B., Liu D., 2020, Research in Astronomy and Astrophysics, 20, 135
  * Webbink (1984) Webbink R. F., 1984, ApJ, 277, 355
  * Webbink (2008) Webbink R. F., 2008, in Milone E. F., Leahy D. A., Hobill D. W., eds, Astrophysics and Space Science Library Vol. 352, Astrophysics and Space Science Library. p. 233 (arXiv:0704.0280), doi:10.1007/978-1-4020-6544-6_13
  * Willems & Kolb (2004) Willems B., Kolb U., 2004, A&A, 419, 1057
  * Williams et al. (2009) Williams K. A., Bolte M., Koester D., 2009, ApJ, 693, 355
  * Wonnacott et al. (1993) Wonnacott D., Kellett B. J., Stickland D. J., 1993, MNRAS, 262, 277
  * Wright et al. (2010) Wright E. L., et al., 2010, AJ, 140, 1868
  * Yao et al. (2020) Yao Y., et al., 2020, ApJ, 900, 46
  * Zorotovic & Schreiber (2022) Zorotovic M., Schreiber M., 2022, MNRAS, 513, 3587
  * Zorotovic et al. (2010) Zorotovic M., Schreiber M. R., Gänsicke B. T., Nebot Gómez-Morán A., 2010, A&A, 520, A86
  * Zorotovic et al. (2011) Zorotovic M., et al., 2011, A&A, 536, L3
  * Zorotovic et al. (2014) Zorotovic M., Schreiber M. R., Parsons S. G., 2014, A&A, 568, L9
  * van den Heuvel (1994) van den Heuvel E. P. J., 1994, A&A, 291, L39
  * van den Heuvel & Tauris (2020) van den Heuvel E. P. J., Tauris T. M., 2020, Science, 368, eaba3282

##  Appendix A Radial velocities

The radial velocities for each of our objects are listed in Tables 7 to 11.

HJD | RV [km/s] | Instrument
---|---|---
2459753.9514 | -100.36 $\pm$ 0.04 | TRES
2459814.7098 | -13.62 $\pm$ 0.04 | FEROS
2459817.6290 | 11.59 $\pm$ 0.04 | FEROS
2459818.6918 | 6.58 $\pm$ 0.04 | FEROS
2459819.7567 | -5.42 $\pm$ 0.04 | FEROS
2459820.7507 | -21.86 $\pm$ 0.06 | FEROS
2459821.6934 | -39.98 $\pm$ 0.04 | FEROS
2459822.6943 | -59.91 $\pm$ 0.12 | FEROS
2459823.7031 | -78.14 $\pm$ 0.12 | FEROS
2459824.6181 | -91.25 $\pm$ 0.05 | FEROS
2459825.7693 | -101.13 $\pm$ 0.09 | FEROS
2459826.7605 | -102.44 $\pm$ 0.06 | FEROS
2459827.6655 | -97.77 $\pm$ 0.05 | FEROS
2459828.7288 | -85.67 $\pm$ 0.05 | FEROS
2459889.6855 | 10.91 $\pm$ 0.04 | TRES
2459903.5700 | -27.14 $\pm$ 0.06 | FEROS
2460074.9628 | -83.22 $\pm$ 0.04 | TRES
2460080.8559 | -65.56 $\pm$ 0.05 | FEROS
2460091.9633 | -67.47 $\pm$ 0.06 | TRES
2460093.9142 | -95.82 $\pm$ 0.05 | TRES
2460095.9101 | -101.34 $\pm$ 0.04 | TRES
2460097.9193 | -80.48 $\pm$ 0.03 | TRES
2460100.9396 | -23.21 $\pm$ 0.04 | TRES
2460102.9528 | 5.82 $\pm$ 0.04 | TRES
2460106.9231 | -11.16 $\pm$ 0.04 | TRES
2460107.9399 | -29.31 $\pm$ 0.04 | TRES
2460108.9243 | -48.74 $\pm$ 0.04 | TRES
2460111.8596 | -96.28 $\pm$ 0.04 | FEROS
2460140.7176 | 10.41 $\pm$ 0.06 | FEROS
Table 7: RVs for J2117+0332. HJD | RV [km/s] | Instrument
---|---|---
2459900.0100 | -61.22 $\pm$ 0.04 | TRES
2459912.9834 | -21.44 $\pm$ 0.04 | TRES
2459924.9851 | 3.19 $\pm$ 0.04 | TRES
2459934.9399 | -76.47 $\pm$ 0.02 | TRES
2459951.9817 | 19.79 $\pm$ 0.04 | TRES
2459954.9640 | 15.24 $\pm$ 0.04 | TRES
2459958.9732 | -11.88 $\pm$ 0.03 | TRES
2459963.9729 | -58.62 $\pm$ 0.05 | TRES
2459970.9312 | -73.19 $\pm$ 0.05 | TRES
2459972.9333 | -60.89 $\pm$ 0.04 | TRES
2459974.9553 | -43.56 $\pm$ 0.04 | TRES
Table 8: RVs for J1111+5515 HJD | RV [km/s] | Instrument
---|---|---
2459925.0123 | 9.44 $\pm$ 0.04 | TRES
2459936.0571 | -45.70 $\pm$ 0.04 | TRES
2459951.9913 | 22.62 $\pm$ 0.03 | TRES
2459970.9937 | 5.95 $\pm$ 0.04 | TRES
2459975.0058 | -22.33 $\pm$ 0.03 | TRES
2459979.0157 | -41.33 $\pm$ 0.03 | TRES
2459982.0153 | -45.77 $\pm$ 0.04 | TRES
2459988.0460 | -31.05 $\pm$ 0.03 | TRES
2459991.8859 | -10.68 $\pm$ 0.03 | TRES
2460000.8443 | 39.02 $\pm$ 0.03 | TRES
2460007.9507 | 49.84 $\pm$ 0.03 | TRES
Table 9: RVs for J1314+3818 HJD | RV [km/s] | Instrument
---|---|---
2459813.7415 | 24.14 $\pm$ 0.04 | FEROS
2459817.6103 | 48.93 $\pm$ 0.04 | FEROS
2459818.8187 | 55.94 $\pm$ 0.06 | FEROS
2459819.7716 | 61.08 $\pm$ 0.03 | FEROS
2459820.7658 | 65.91 $\pm$ 0.04 | FEROS
2459822.7247 | 73.31 $\pm$ 0.05 | FEROS
2459823.7338 | 76.15 $\pm$ 0.05 | FEROS
2459824.7085 | 77.91 $\pm$ 0.05 | FEROS
2459825.7838 | 78.86 $\pm$ 0.06 | FEROS
2459826.7799 | 78.97 $\pm$ 0.03 | FEROS
2459827.6420 | 78.25 $\pm$ 0.03 | FEROS
2459830.6480 | 70.78 $\pm$ 0.04 | FEROS
2459832.6634 | 61.90 $\pm$ 0.05 | FEROS
2459898.5129 | -11.85 $\pm$ 0.03 | FEROS
2459904.5650 | 15.34 $\pm$ 0.04 | FEROS
2459915.5389 | 75.04 $\pm$ 0.07 | FEROS
2459920.5358 | 77.22 $\pm$ 0.04 | FEROS
2460050.8438 | 63.57 $\pm$ 0.03 | FEROS
2460076.8618 | -11.34 $\pm$ 0.03 | FEROS
2460077.8044 | -13.52 $\pm$ 0.03 | FEROS
2460113.7077 | 36.52 $\pm$ 0.02 | FEROS
2460139.7119 | 44.55 $\pm$ 0.02 | FEROS
Table 10: RVs for J2034-5037 HJD | RV [km/s] | Instrument
---|---|---
2459813.9113 | 14.98 $\pm$ 0.03 | FEROS
2459817.7473 | 35.35 $\pm$ 0.04 | FEROS
2459818.8898 | 41.02 $\pm$ 0.03 | FEROS
2459819.8547 | 45.46 $\pm$ 0.03 | FEROS
2459820.7829 | 49.27 $\pm$ 0.04 | FEROS
2459821.7682 | 52.86 $\pm$ 0.03 | FEROS
2459822.7691 | 55.67 $\pm$ 0.07 | FEROS
2459823.7667 | 57.84 $\pm$ 0.07 | FEROS
2459824.8124 | 58.92 $\pm$ 0.04 | FEROS
2459825.8175 | 58.90 $\pm$ 0.06 | FEROS
2459826.8445 | 58.08 $\pm$ 0.03 | FEROS
2459827.8081 | 56.18 $\pm$ 0.04 | FEROS
2459828.7738 | 53.31 $\pm$ 0.03 | FEROS
2459829.7258 | 49.85 $\pm$ 0.04 | FEROS
2459898.7135 | -27.52 $\pm$ 0.03 | FEROS
2459904.6857 | -16.73 $\pm$ 0.03 | FEROS
2459922.5500 | 58.70 $\pm$ 0.03 | FEROS
2460112.8601 | 40.72 $\pm$ 0.03 | FEROS
2460140.8231 | -26.44 $\pm$ 0.03 | FEROS
2460158.8008 | 24.81 $\pm$ 0.03 | FEROS
Table 11: RVs for J0107-2827

##  Appendix B SED contribution from WDs

As a test of whether the WD is contributing to the observed photometry, we
generated Koester WD models (Tremblay & Bergeron, 2009; Koester, 2010) and
added them to the fitted SEDs from Section 3.6. The photometry with and
without this addition was calculated and compared. For the WD models, we set
log(g) = 9.0 given the large masses of our WDs, and tested a range of
temperatures, $T_{\rm eff,WD}=10000\ -\ 60000$ K.

For J1314+3818, which is the most intrinsically faint object amongst the five,
we find that the SDSS u band photometry changes by more than 0.2 mags if
$T_{\rm eff,WD}\gtrsim 30000$ K. This is large compared to the typical errors
on the the photometry which is about 0.02 mags. This is shown on Figure 9 for
$T_{\rm eff,WD}=40000$ K where the bottom panel shows that the difference is
$\sim 0.25$ mags in this case. Therefore, to avoid possible contamination,
this point was excluded in our final SED fitting described in Section 3.6.

For all other objects, we find no significant contribution, with the SDSS u
band photometry changing by less than 0.03 mags for $T_{\rm eff,WD}=60000$ K.
Furthermore, we note that the WD cooling timescale gets longer at lower
temperatures meaning that they are more likely to be cool. All this
considered, we have included the SDSS u band in our SED fitting for these
objects.

Figure 9: Plot of the fitted SED for J1314+3818 as shown in Figure 2 (gray
solid line), along with the model SED for a WD of $T_{\rm eff,WD}=40,000$
(gray dashed lines). The sum of these two SEDs are shown in blue. The lower
plot shows the difference in photometry calculated using the gray solid and
the blue solid lines. We see that this difference is $\sim 0.2$ mags for the
SDSS u band.