Wide post-common envelope binaries containing ultramassive white dwarfs: evidence for efficient envelope ejection in massive asymptotic giant branch stars

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_orb ≲ 1 d), implying relatively inefficient harnessing of binaries’ orbital energy for envelope expulsion. Here, we present follow-up observations of five binaries from 3rd data release of Gaia mission containing solar-type MS stars and probable ultramassive WDs (⁠|$M\gtrsim 1.2\ {\rm M}_{\odot}$|⁠) with significantly wider orbits than previously known PCEBs, P_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. Our one-dimensional stellar models including recombination energy confirm prior predictions that a wide range of PCEB orbital periods, extending up to months or years, can potentially result from Roche lobe overflow of a luminous asymptotic giant branch (AGB) star. This evolutionary scenario may also explain the formation of several wide WD + MS binaries discovered via self-lensing, as well as a significant fraction of post-AGB binaries and barium stars.

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, spiralling inward on a dynamical or thermal time-scale. 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. Modelling 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 (NSs; 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, Tout & Pols 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 (‘α’) 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 modelled 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 (WDs) or hot subdwarfs in tight orbits (P_orb ≲ 1 d) with a main-sequence (MS) star, having been produced when the MS star spiralled 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 α-formalism (Han et al. 2002, 2003; Camacho et al. 2014). Such modelling makes it possible to empirically constrain α, and several calculations have found that the observations can best be reproduced by models assuming α ≈ 0.3, meaning that |$\approx\!30{{\ \rm per\,cent}}$| of the orbital energy liberated during inspiral goes into ejecting the envelope (Zorotovic et al. 2010; Davis, Kolb & Knigge 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 d and hosts an unusually massive WD, with mass |$M\approx 1.2\ {\rm M}_{\odot}$| (Wonnacott, Kellett & Stickland 1993). The system’s wide orbit (a ≈ 0.2 au) means that less orbital energy was liberated during the MS star’s inspiral than in typical PCEBs with a ≈ 0.01 au. Zorotovic et al. (2010) found that in the α-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 α values) to explain (Zorotovic, Schreiber & Parsons 2014). Energy released by H and He recombination in the expanding envelope of the WD progenitor is a prime suspect for supplying the additional energy (first explored by Paczyński & Ziółkowski 1968, and later studied by e.g. Webbink 2008; Ivanova, Justham & Podsiadlowski 2015; Ivanova 2018)

Most PCEBs studied to date were identified via their composite spectra and radial velocity (RV) variability detectable with low-resolution spectra (Rebassa-Mansergas et al. 2007, 2017; Lagos et al. 2022). This leads to strong selection effects in favour 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 (DR3) of the Gaia mission (Gaia Collaboration 2023a) contains orbital solutions for more than 10⁵ astrometric binaries, and for more than 10⁵ single-lined spectroscopic binaries identified from medium-resolution spectra (Gaia Collaboration 2023c). This data set 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 MS stars and probable ultramassive WD candidates. Section 2 describes our identification of wide PCEB candidates from the Gaia DR3 catalogue. Section 3 describes follow-up spectroscopic observations to obtain RVs, spectral analysis to calculate metallicities, and fitting to the broad-band spectral energy distributions (SEDs) 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 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) catalogue (Gaia Collaboration 2023c). We selected promising candidates for further follow-up based on their mass functions, colour–magnitude diagram (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 per cent 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 five promising sources that are single-lined and whose Gaia-reported orbits were validated by our follow-up. All five 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. The names, Gaia DR3 source IDs, and basic information of these five objects are summarized in Table 1. 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 (SNR) applied to astrometric solutions with short periods (see Gaia Collaboration 2023c). 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. (2023b) on the basis of its large astrometric mass ratio function (also see Shahaf et al. 2023a). Another object in our sample, J2034−5037, was previously identified by Jayasinghe et al. (2023) as a candidate NS + MS binary.

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 SEDs. A log of our observations and measured RVs can be found in Appendix  A.

3.1 Fibrefed Extended Range Optical Spectrograph

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

We reduced the data using the CERES pipeline (Brahm, Jordán & Espinoza 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 fibre. 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, with log(g) = 4.5 and solar metallicity.

3.2 Tillinghast Reflector Echelle Spectrograph

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 atop Mount Hopkins, Arizona. TRES is a fibrefed echelle spectrograph with a wavelength range of 390–910 nm and spectral resolution R ≈ 44 000 (1 × 1 binning). Exposure times ranged from 1800 to 3600 s. 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 Magellan Inamori Kyocera Echelle

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 600 s. This yielded a spectral resolution R ∼ 35 000 (2 × 2 binning) and typical SNR of approximately 35 and 16 per pixel on the red and blue side, respectively. The total wavelength coverage was ∼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. We co-added the two spectra obtained for J2034−5037 across two nights (HJD 2460092.7715 and 2460118.7732) in the same way. J2117+0332 was observed once (HJD 2460118.8313).

3.4 Metallicities

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

3.4.1 Stellar Parameter Classification

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–5360 Å, centred on the Mg I b triplet. It then fits the peaks of the cross-correlation function with a three dimensional (3D) third-order polynomial to return best-fitting values of effective temperature T_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_eff values are ∼0.08 dex and ∼50 K, respectively (See also Furlan et al. 2018). We report these floors in Table 2.

3.4.2 Brussels Automatic Code for Characterizing High accUracy Spectra

For the MIKE and FEROS spectra, we used the Brussels Automatic Code for Characterizing High accUracy Spectra (BACCHUS, Masseron, Merle & Hawkins 2016; Hayes et al. 2022). This code performs 1D local thermodynamic equilibrium (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 with 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–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 with those calculated by Andrae, Rix & Chandra (2023) using the Gaia XP very low-resolution spectra. These authors derive T_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, Rix & Chandra (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 metallicities 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, Rix & Chandra (2023) within 1σ. 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-fitting [Fe/H] values from our SED fitting (Section 3.6), which uses the SPC and BACCHUS metallicities as a prior.

3.5 Light curves

We retrieved observed light curves for our objects from the All-Sky Automated Survey for Supernovae (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 ∼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 ∼0.002–0.003 for all objects except J1314+3813, where they have amplitude ∼0.006.

3.6 SED fitting

We constructed broad-band SEDs of our targets using synthetic ugriz Sloan Digital Sky Survey (SDSS) photometry calculated from Gaia XP spectra (Gaia Collaboration 2023b; 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₁W₂W₃ 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, Clayton & Mathis (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₀ 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_eff = 60 000 K would not significantly contribute to the photometry of all but one of our targets. For the one exception, J1314+3818, we find that a WD with T_eff ≳ 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 modelling of stellar photometry and spectra. We only use the code’s photometric modelling 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_⋆, 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 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_⋆, [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. While values for T_eff are also obtained from spectral analysis (Table 2), given the degeneracy that can exist between T_eff and log g in spectroscopic fits and the high quality of the SED fits, we do not use them to constrain the outputs here.

Putting everything together, the final likelihood function is

where ‘mag’ stands for apparent magnitudes and the summation is over the appropriate photometric filters for each object, σ_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 σ_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_⋆, EEP, and [Fe/H]_init are the parameters directly fitted by MINEsweeper, while [Fe/H], T_eff, and R are calculated from the isochrones corresponding to the fitted parameters. The fit to parallax and the reported errors are described in Section 3.6.2. We also list constraints on [Fe/H] and T_eff for comparison with the values measured from spectroscopy (Table 2).

Fig. 1 shows MESA Isochrones & Stellar Tracks (MIST) isochrones corresponding to the stellar parameters of 100 random posterior samples (grey). The cyan lines show the best-fitting parameters and the red point marks the present inferred parameters of the MS stars. The labels also indicate the stellar ages, which range from 1.84 to 11.57 Gyr. Two systems, J111+5515 and J0107−2827, host stars that have slightly evolved off the MS. This is likely to be the result of selection bias, as evolved stars are brighter and thus overrepresented 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 Fig. 2. The model SEDs plotted were generated using pytstelllibs¹ with the best-fitting 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. As mentioned above, for J1314+3818, we found that WDs with |$T_{\rm eff} \gtrsim 30\, 000$| K would significantly contribute to the SDSS u-band photometry (Appendix  B), so we excluded this point in fitting. GALEX near-ultraviolet (NUV) observations are shown in Fig. 2, but these were also excluded from the fitting for all targets to avoid potential contamination from the WD companion. We investigate the expected contributions of WD companions to the NUV photometry in Appendix  B. There we show that the observed NUV photometry is consistent with no contributions from the WDs for all our targets. This places an upper limit on the effective temperatures of the WDs, ranging from |$T_{\rm eff,\, WD}\ \lt\ 7750\ {\rm K}$| for J1314+3818 to |$T_{\rm eff,\, WD}\ \lt\ 60\,000\ {\rm K}$| for J1111+5515.

3.6.1 A wide tertiary

One system, J0107−2827, has a resolved tertiary separated by a distance of 2.21 arcsec, corresponding to a projected physical separation of 1033 au (Gaia DR3 ID 5033 197 892 724 532 608). 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, Rix & Heintz 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 MS, its G-band absolute magnitude of ∼7.3 (calculated using the reported apparent magnitude of 15.7) corresponds to a mass of approximately |$0.67\ {\rm M}_{\odot}$|⁠. We assume solar metallicity (consistent with the initial metallicity we infer for the primary) and an age of ∼6 Gyr. Using these parameters, we generated photometry for this third star (grey line in Fig. 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, ϖ, free, allowing us to propagate parallaxes uncertainties through to the stellar parameters. From Table 1, we see that the Renormalized 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 catalogue 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 2023c; Halbwachs et al. 2023). From these, we select those with phot_g_mean_mag < 13 and RUWE values comparable with our targets. We then calculate the standard deviation of the difference between the parallaxes reported from single-star solutions (gaia_source) and binary solution (NSS catalogue). 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-fitting 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.1, we also set an uncertainty floor of |$\pm 0.05\ {\rm M}_{\odot}$| on M_⋆ to obtain conservative errors on the inferred masses of the WDs.

4. ORBITAL FITS AND CONSTRAINTS ON THE UNSEEN COMPANION

We fit the FEROS and TRES RVs with a Keplerian model using emcee. The free parameters of the fit are the orbital period P_orb, periastron time T_p, eccentricity e, RV semi-amplitude K_⋆, argument of periastron ω, and centre-of-mass RV γ. 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

where RV_pred(t_i) and 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-fitting RV curve for each object is shown in Fig. 3. Best-fitting values for P_orb, e, and K_⋆, 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 ω fixed to zero. The residuals from the two models are plotted on the second and third panels for each object in Fig. 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 ± 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 ∼ 0.0003 ± 0.0002. This is comparable with 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.

4.1 Masses of the unseen companion

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 companion M_WD. (Note that this notation implies that the unseen companion is a WD which we argue is the most likely scenario in Section 4.3. We keep this notation here for consistency throughout the paper):

where M_⋆ 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° (i.e. when the orbit is ‘edge-on’ to our line of sight).

The implied M_WD as a function of the inclination is shown in Fig. 4. We shade the regions for M_⋆ values |$\pm 0.05\ {\rm M}_{\odot}$| above and below the best-fitting value from the SED fitting (Section 3.6.2). The minimum masses M_{WD, min} range from 1.244 ± 0.027 to |$1.418\pm 0.033\ {\rm M}_{\odot}$|⁠. Given the uncertainties, these values are all consistent with masses close to, but just below, the Chandrasekhar limit of |$\sim\!1.4\ {\rm M}_{\odot}$|⁠.

For J1314+3818, we obtain an inclination constraint from astrometry (Section 4.2) and thus a precise value for M_WD of |$1.324 \pm 0.037\ {\rm M}_{\odot}$|⁠, as opposed to just a lower limit. This point is shown as a red cross on the plot of the M_WD(i) for this object in Fig. 3.

The inferred WD masses are summarized in Table 5. At the time of writing, these WDs are among the most massive WDs known (e.g. Hermes et al. 2013; Cognard et al. 2017; Curd et al. 2017; Hollands et al. 2020; Caiazzo et al. 2021; Miller et al. 2022), if they are indeed WDs (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 MS companions).

4.2 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, centre_of_mass_velocity, eccentricity, period, and t_periastron (Gaia Collaboration 2023c; Halbwachs et al. 2023). The Thiele–Innes elements A, B, F, and G describe the astrometric orbit of the photocentre 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 photocentre 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 Ω, argument of periastron ω, periastron time T_p, centre-of-mass RV γ, companion mass M_WD, and luminous star mass M_⋆. For M_⋆, we set a Gaussian constraint based on its best-fitting value and error obtained from the SED fitting (Table 3). The resulting parameters can be found in Table 6 and the plots are shown in Fig. 3. We find i ≈ 99.9° and |$M_{\rm WD} = 1.324 \pm 0.037\ {\rm M}_{\odot}$|⁠.

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\ {\rm 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\ {\rm 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\ {\rm 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\ {\rm M}_{\odot}$| star on the MS. This mass roughly corresponds to a K7V star with a radius and temperature of 0.63R_⊙ 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 per cent, respectively, for J2117+0332, J111+5515, J1314+3818, J2034−5037, and J0107−2827. In the infrared, at 3 |$\mu$|m, the contribution would be larger, ranging from ∼13 to |$56{{\ \rm per\,cent}}$| 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). The same challenges would apply for a massive WD formed from a WD + WD merger. Therefore, we do not consider these options 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-fitting 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 ≲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_out/P_in ≳ 5 where P_out and P_in are the outer and inner orbital periods, respectively (we are also taking e ∼ 0 and the ratio of the mass of the outer star to that of the inner binary ∼1, Mardling & Aarseth 2001; Tokovinin 2014). Given that our objects have P_orb ∼ 30 d, this implies P_in ≲ 6 d. 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 (Prša & Zwitter 2005) to generate synthetic light curves of an inner binary of two K dwarfs for a range of periods from ∼0.25 to 1 d. 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 ≲10 per cent 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 ∼10 per cent of the light coming from the inner binary, ellipsoidal variability can only be distinguished from the noise for P_inner ≲ 0.3 d. Thus, the range of possible inner period is ∼0.3–6 d.

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 ∼ 1000 d (this is not a selection effect, see Tokovinin 2014), while all our systems have P_orb < 50 d. 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 five 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 5), we also consider the possibility that the companions are 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, Pfahl & Rappaport 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

where Δm is the mass lost, m_c is the remaining core/NS mass, and m₂ is the mass of the companion (Hills 1983). In the case where a 8-|${\rm M}_{\odot}$| star explodes by a core-collapse supernova (SN) to form a 1.3-|${\rm M}_{\odot}$| NS around a 1-|${\rm 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 interactions 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 ultrastripped SN explosion with |$\sim\!0.3\ {\rm M}_{\odot}$| of ejecta (De et al. 2018; Yao et al. 2020), we expect e ∼ 0.1, significantly larger than the majority of our systems (though even lower ejecta masses are possible, Tauris et al. 2013; Tauris, Langer & Podsiadlowski 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, an 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 × 10⁻³ to |$0.05\ {\rm M}_{\odot}$| (Fryer et al. 1999; Darbha et al. 2010). Such ejecta masses could correspond to e ≲ 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 Roche lobe overflow (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 with other related classes of binaries, including other WD + MS PCEBs and WD + millisecond pulsar (MSP) binaries.

5.1 Literature PCEBs

The 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 via ultraviolet (UV) excess and RV variability detected from the RAVE and LAMOST surveys (diamond markers, Hernandez et al. 2021, 2022a, b). In Fig. 5, we plot the minimum orbital separation a_peri against M_WD for these literature PCEBs, the five objects from our sample, and several self-lensing binaries (SLBs, i.e. WD + MS binaries) discovered by the Kepler survey (SLBs + KOI-3278, Kruse & Agol 2014; Kawahara et al. 2018). With the exception of J1314+3818 where the precise value of M_WD was obtained using astrometry (Section 4.2), we have plotted M_{WD, min} for our objects which is indicated with arrows. The colours of the points represent the mass of the luminous (MS) companion, M_⋆. The grey dashed line comes from the P_orb–M_WD relation derived in Rappaport et al. (1995) for stable MT (with a spread in orbital period of a factor of ∼2.4), where P_orb has been converted to separation assuming a 1-|${\rm M}_{\odot}$| MS star and circular (e = 0) orbit (assuming instead a 0.1-|${\rm 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 Fig. 5 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_init using the WD Initial–Final Mass Relation (IFMR) derived in Williams, Bolte & Koester (2009): M_final = 0.339 + 0.129M_init, then generating MIST evolutionary tracks (Choi et al. 2016; Dotter 2016) to identify the maximum radius reached by a star with a given M_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_orb–M_WD parameter space: being in a wider orbit with a period of 22 d and hosting a more massive WD of |$\sim\!1.2\ {\rm M}_{\odot}$| (Wonnacott, Kellett & Stickland 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 asymptotic giant branch (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, Schreiber & Parsons (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 Fig. 5 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, b) (plotted with diamond markers) from the SDSS PCEBs as they host higher mass (⁠|$\sim\!1\ {\rm M}_{\odot}$|⁠) MS stars. Compared with our objects, these have shorter orbital periods (P_orb ∼ 1–2 d) 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 Fig. 5 may suggest that there are two distinct groups of PCEBs in wide orbits with different WD masses (SLBs versus IK-Peg analogues), 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.

5.2 Eccentricities and comparison with MSP binaries

On Fig. 6, we plot the eccentricity e against P_orb and compare our objects with MSP + WD binaries. We plot objects primarily from the Australia Telescope National Facility (ATNF) catalogue (Manchester et al. 2005) (taking the version analysed by Hui et al. 2018) and distinguish those with minimum WD companion masses above and below |$0.45\ {\rm M}_{\odot}$| (this is the 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 are formed 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.

Fig. 6 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, Langer & Kramer 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 Fig. 6. Both of these systems contain mildly recycled pulsars (with spin periods of 27 and 32 ms) and massive WDs (⁠|$M_{\rm WD}\ \gt\ 0.8\ {\rm M}_{\odot}$|⁠), and they have eccentricities of 0.045 and 0.017 at orbit periods of 40 and 32 d. 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, van den Heuvel & Savonije 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 Fig. 6. Interestingly, these eMSPs occupy a narrow range in orbital periods with eccentricities that are comparable with the two eMSP + 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, or 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 disc formed from material lost from the WD during H shell flashes (Antoniadis 2014), or a circumbinary disc formed from the ejected envelope after a ce phase (Dermine et al. 2013).

5.3 Other related populations

Several other binary populations have been identified as likely products of CEE but have been largely neglected by works inferring CEE parameters from observed PCEB populations. We discuss them and their possible relation to the objects in our sample below.

First, the ‘post-AGB’ binaries (van Winckel 2003; Oomen et al. 2018). Post-AGB stars are short-lived transitional objects that have lost most all of their envelopes during the AGB phase and are transitioning towards hotter temperatures at constant luminosity on their way to the WD cooling track. Most of the known post-AGB stars have temperatures of T_eff = 5000–8000 K and radii of |$20{-}50\, R_{\odot }$|⁠. About 50 post-AGB stars are known to be in binaries, with periods ranging from approximately hundred to a few thousand days. These periods are too short for the binaries to have avoided interaction during the AGB phase, and thus post-AGB binaries are generally modelled as having formed through CEE (e.g. Izzard & Jermyn 2018). Like the objects in our sample, their wide orbits would require very efficient envelope ejection to explain.

It has also been established that a significant fraction (a lower limit of |$10{-}20{{\ \rm per\,cent}}$|⁠, Bond 2000; Miszalski et al. 2009) of the central stars of planetary nebulae (CSPNe) are in binaries (for a recent review, see Jones & Boffin 2017). The majority of the known systems were discovered by their photometric variability and have orbital periods ≲1 d (e.g. Miszalski et al. 2009; Jacoby et al. 2021). Given the short lifetimes of PNe, these should have undergone the CEE very recently with little subsequent evolution. Post-AGB stars are immediate progenitors of PNe, so the striking mismatch between the long periods of post-AGB binaries and short periods of binary CSPNe is puzzling. We suspect that it owes mostly to observational selection effects: Post-AGB stars are too large and puffy to fit inside short-period orbits, while longer period binary CSPNe are difficult to detect. If this explanation holds, a significant population of intermediate- and long-period binary CSPNe should still await detection. Indeed, a few binary CSPNe have been found with longer orbital periods of ∼10–30 d (e.g. Brown et al. 2019; Jacoby et al. 2021), comparable with those of IK Peg and our five systems, and with long-term RV monitoring, a few with orbital periods of thousands of days have also been detected (Van Winckel et al. 2014; Jones et al. 2017).

Another related class of binaries is the barium stars, which are red giants found in binaries that show an enhancement in s-process elements, such as barium. Their peculiar abundances are thought to be the result of the star accreting material from an AGB star companion though winds (Boffin & Jorissen 1988) or RLOF (McClure 1983). The AGB star then evolves into a WD, forming a PN in the process and leaving behind the barium stars in relatively wide orbits at the centres of PNe. A few such systems have indeed been found (e.g. Miszalski et al. 2012, 2013), while many more have been found with similar orbital periods that are not in PNe (e.g. Jorissen et al. 2019). Most known barium stars have periods of 200–5000 d: comparable with the post-AGB binaries, but wider than the five objects in our sample. Here too, the fact that the polluted stars are giants likely introduces a significant selection bias in favour of long-period systems, and short-period analogues hosting MS stars have recently been identified (Roulston et al. 2021).

Finally, a subset of the low-eccentricity binaries with periods of ∼200–1500 d discovered via phase modulation of δ scuti pulsations (Murphy et al. 2014, 2018) are likely to contain WD companions, many of which likely formed through a similar channel to the barium and post-AGB stars. Since the amplitude of phase modulations depends on the physical size of the orbit, selection effects in this sample likely also favour long periods.

To summarize, the orbital periods of the binaries in our sample fall between the tight binaries that have previously been used to constrain CEE models (WD + MS PCEBs and binary CSPNe), and the long-period post-AGB and barium star binaries, which may have formed through either CEE or wind accretion. Given the complex selection effects affecting all the observed post-interaction binaries, more work is required to understand how these populations are related and to infer the intrinsic post-CEE period distribution.

6. FEASIBILITY OF FORMATION THROUGH CEE

To test whether our targets could have formed via CEE, we ran Modules for Experiments in Stellar Astrophysics (MESA) models (Paxton et al. 2011, 2013, 2015, 2018) of the progenitor star to the WD in our systems, evolving an intermediate-mass star up to the AGB. We emphasize that we are not using binary models to trace the evolution during the ce phase, which is beyond the capabilities of MESA. We are simply constructing a realistic model of the giant at the onset of MT to calculate the energy budget as described below.

In the α-formalism, the CEE ends when the loss of orbital energy from the spiral-in exceeds the binding energy of the envelope E_bind, resulting in the envelope being ejected:

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 {\small CE}}$| is the fraction of the liberated orbital energy that goes into unbinding the envelope. Thus, E_bind determines the final separation of the system for a given initial separation.

In the simplest case, E_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 {\small CE}} \lt 1$| can reproduce the relatively wide orbits of IK Peg and KOI-3278 (Davis, Kolb & Willems 2010; Zorotovic et al. 2010; Zorotovic, Schreiber & Parsons 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. Paczyński & Ziółkowski 1968; Ivanova et al. 2013).

From IFMRs of WDs (e.g. Williams, Bolte & Koester 2009; Cummings et al. 2018; El-Badry, Rix & Weisz 2018), we expect the MS progenitor masses of our ultramassive WDs to be in the range of |$6{-}9\ {\rm M}_{\odot}$|⁠. Thus, we run MESA models of a 7-|${\rm M}_{\odot}$| star, following its evolution up to the AGB. Our inlists are based on those of Farmer, Fields & Timmes (2015) (same wind and mixing prescriptions; for more details, we refer readers to section 2 of Farmer, Fields & Timmes 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-|${\rm M}_{\odot}$| model, from pre-MS to termination at the tip of the AGB, is shown on the Hertzsprung-Russell (HR) diagram in the leftmost panel of Fig. 7. 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 time-step, we calculate the binding energy of the envelope as a sum of the gravitational (E_grav) and internal (E_int) components:

where E_grav and E_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_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, both the thermal and recombination components, 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 RLOF 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_WD, is taken to be equal to the helium core mass of the giant, M_core. This quantity grows from ∼1.2 to |$1.7\ {\rm M}_{\odot}$| on the RGB and falls back down to |$\sim\!1.1\ {\rm M}_{\odot}$| on the AGB, when some of the helium is mixed back into the envelope during second dredge-up (e.g. Busso, Gallino & Wasserburg 1999). We assume |$M_{\star } = 1\ {\rm 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 Fig. 7. In the central panel, the grey dashed line marks a_f = 0.15 au, which is slightly smaller than the minimum a_peri of our objects at ∼0.18 au (Fig. 5). The red dashed line marks a_f = 0.01 au (⁠|$\sim\!2\, R_{\odot }$|⁠), below which a |$\sim\!1$|-M_⊙ 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_bind = E_grav) and set |$\alpha _{\rm {\small CE}} = 1$| (i.e. 100 per cent 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_bind = E_grav + E_int, equation 6) and let |$\alpha _{\rm {\small 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 ∼ 3.5–4.4 au across all values of |$\alpha _{\rm {\small CE}}$|⁠, with a narrower range for lower values of |$\alpha _{\rm {\small CE}}$|⁠.

We also ran a 6-M_⊙ WD progenitor model and found qualitatively similar results but with the a_i range for which a_f > 0.15 au shifted to ∼3.2–4.1 au. Given the simplified treatment here, these small quantitative differences should not be overinterpreted (and hence we do not explore them further) but it does tell us that there is likely a broader range of initial separations for which wide PCEBs can be formed than that implied by a single model.

In Fig. 7, we exclude models in regions of parameter space where the α-formalism does not make clear predictions, namely, those in which E_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 Fig. 7, these conditions are realized late in the AGB evolution, corresponding to radii of ∼520–888R_⊙ and initial separations of ∼4.4–7.7 au.

6.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 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, Sabach & Soker 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. The effectiveness of recombination energy in widening the final orbit has also been explored in hydrodynamic simulations, with a range of results (e.g. Reichardt et al. 2020; González-Bolívar et al. 2022). 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:

where E_th and E_rec are the thermal and recombination energies, and α_th and α_rec are the respective efficiencies. As described above, we set α_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/ρ, where P is the pressure and ρ 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 α_rec would be required to produce wide orbits.

As shown in Fig. 8, we find that for a canonical value of |$\alpha _{\rm {\small CE}} = 0.3$|⁠, we require α_rec ≳ 0.5, and for |$\alpha _{\rm {\small CE}} = 1$|⁠, we require α_rec ≳ 0.25 for a_f ≳ 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 ≈2 × 10⁵ K (Ivanova et al. 2013). This boundary lies deep in the envelope of our models on the AGB, at roughly 20 per cent 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, Moreno Méndez, López-Cámara & De Colle 2017; Sabach et al. 2017). We defer more detailed calculations and further discussion of this topic to future work.

6.2 Models for lower and higher mass giants

We performed similar calculations for a 1-M_⊙ red giant, which might be a progenitor to a wide PCEB hosting a lower mass (⁠|$0.5{-}0.6\ {\rm M}_{\odot}$|⁠) WD. We evolved a |$1\ {\rm 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, Fields & Timmes (2015).

The same plots as in Fig. 7 but for the 1-M_⊙ model are shown in Fig. 9. Even in the case where we consider only the gravitational binding energy with |$\alpha _{\rm {\small CE}} = 1$|⁠, there is a range of initial separations (∼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 {\small 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\ {\rm M}_{\odot}$| WDs (such as the SLBs) 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. Girardi & Marigo 2007; Marigo & Girardi 2007), so this conclusion may depend on the adopted stellar models.

With the inclusion of all of the 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 ∼0.9–4.5 au. For a_i ∼ 3–4.5 au, the envelope’s binding energy is positive when recombination energy is included, so the α-formalism does not straightforwardly predict a final separation. For separations a_i ≲ 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\ {\rm M}_{\odot}$|⁠, that become red supergiants and will leave behind NSs or black holes. We find that the envelopes of these stars are significantly more bound than that of the 7-M_⊙ 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).

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 MT is unstable, q_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\ {\rm M}_{\odot}$| and an intermediate MS companion of |$\sim\!1\ {\rm M}_{\odot}$|⁠, the mass ratio is large enough to exceed even conservative values of q_crit ∼ 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 MT is unlikely to have formed the systems in our sample.

We would be remiss here to not mention the ‘γ’-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 α-formalism (Nelemans et al. 2000; Nelemans & Tout 2005). The parameter γ 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 γ-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 with 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 Fig. 5) 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_⊙ MS star orbiting a dark companion that is almost certainly a WD in a period of 14.7 d, 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 (≲1 d), 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 favour of low-mass (M dwarf) MS 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 favour 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 PCEBs containing ultramassive WD candidates and intermediate mass MS stars with long orbital periods (18–49 d). These objects were discovered as part of a broader search for compact object binaries from the Gaia DR3 NSS catalogue. 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 short-period orbits. Systems like the ones presented here pose a potential challenge in simplified models of 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:

• Nature of the unseen companions: The companions are dark objects with masses of |$1.2{-}1.4\ {\rm M}_{\odot}$| – more massive than the solar-type stars orbiting them. The simplest explanation is that they are WDs. We consider two possible alternatives: (i) A tight binary containing two |$\sim\!0.65$|-M_⊙ MS stars. In the most pessimistic case, such an inner binary could escape detection in four of our five 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 d. (ii) An 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.

• 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 to |$1.418\pm 0.033\ {\rm 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\ {\rm M}_{\odot}$|⁠. Assuming the dark companions are in fact WDs, they are among the most massive WDs known.

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

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

• Evolutionary models of the WD progenitors: We ran MESA models of a |$7$|-M_⊙ MS star up the RGB and AGB (Section 6), following its internal energy (thermal + recombination) and calculating the expected final separation according to the α-formalism. We find that there is no point in the evolutionary phase of the star where final separations comparable to those of our objects (≳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 (∼3.5–4.5 au) where final separations exceed ∼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.

• Space density: We compare distances of literature PCEBs to those of our objects and a few other wide PCEBs. At ∼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 d, much wider than typical PCEBs in the literature with P < 1 d. 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 with or larger than that of close PCEBs.

  We also note that several binary populations not typically included in PCEB studies, such as post-AGB binaries and short-period barium and carbon stars, also contain systems with wide orbits that can most readily be explained by efficient ce ejection on the AGB (see Section 5.3). For future work, modelling all of these populations together may pave the way to a more complete understanding of MT processes resulting in the formation of close to intermediate WD + MS binaries.

ACKNOWLEDGEMENTS

We would like to thank Matthias Schreiber, Monica Zorotovic, Sterl Phinney, and the anonymous referee for providing detailed 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, as well as Hans-Walter Rix and Eleonora Zari for observing some of our objects under their FEROS programmes. Finally, we thank Matthias Kruckow, Shi-Jie Gao, Orsola De Marco, and Lev Yungelson for pointing us to relevant and interesting references.

NY and KE were supported in part by NSF grant AST-2307232.

This work presents results from the European Space Agency (ESA) space mission Gaia. Gaia data are being processed by the Gaia Data Processing and Analysis Consortium (DPAC). Funding for the DPAC is provided by national institutions, in particular the institutions participating in the Gaia Multi-Lateral Agreement (MLA). The Gaia mission website is https://www.cosmos.esa.int/gaia. The Gaia Archive website is http://archives.esac.esa.int/gaia.

This paper includes data gathered with the 6.5-m Magellan telescopes located at Las Campanas Observatory, Chile.

This work has made use of the VALD data base, 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.

Footnotes

References

APPENDIX A: RADIAL VELOCITIES

The RVs for each of our objects are listed in Tables A1–A5.

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} = 10\,000{-}60\,000$| K. We used the pyphot² package to calculate the synthetic photometry in different filters.

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 mag if T_{eff, WD} ≳ 30 000 K. This is large compared to the typical errors on the photometry which is about 0.02 mag. This is shown on Fig. B1 for T_{eff, WD} = 40 000 K where the bottom panel shows that the difference is ∼0.25 mag 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 mag for T_{eff, WD} = 60 000 K. Furthermore, we note that the WD cooling time-scale 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.

We performed the same exercise for the GALEX NUV photometry. In Table B1, we summarize the minimum T_{eff, WD} above which flux contributions of a WD would change the NUV photometry by more than 0.1 mag. We see that WD companions could appreciably change the predicted NUV flux above T_{eff, WD} ∼ 7750 K for J1314+3818. The other targets contain MS stars that are brighter in the UV, and so the WDs would be detectable only if they were significantly hotter, with effective temperatures above T_{eff, WD} ∼ 19 750–60 000 K. According to the WD cooling models from Bédard et al. (2020)³, these limits correspond to minimum cooling ages of roughly 3 Gyr for J1314+3818, and 10–700 Myr for the other objects.

To avoid potential contamination from the WDs, we exclude the GALEX NUV points in fitting the SEDs. However, as described in Section 3.6, the observed NUV photometry is nevertheless in good agreement with predictions of the best-fitting single-star model, placing a lower limit on the ages of the WDs. We note that in Fig. 2, J1314+3818 may appear to have UV excess, but this owes mostly to the red leak in the GALEX NUV bandpass: Calculating synthetic photometry for the best-fitting single-star model, we find that the excess is only 0.13 mag, which is similar to the uncertainty in the observed photometry of 0.12 mag. If this excess is real, it could be explained by a 1.25-|${\rm M}_{\odot}$| WD with T_{eff, WD} ∼ 7750 K. Such a WD would have an FUV magnitude of 27.3, which is well below the GALEX detection limit. Deeper UV observations (e.g. with the Hubble Space Telescope) would be needed to assess whether a significant excess is present. Similarly, J0107−2827 has a GALEX excess that is insignificant given the uncertainty.