Wide post-common envelope binaries from Gaia: orbit validation and formation models

Astrometry from {\it Gaia} DR3 has enabled the discovery of a sample of 3000+ binaries containing white dwarfs (WD) and main-sequence (MS) stars in relatively wide orbits, with orbital periods $P_{\rm orb} = (100-1000)$ d. This population was not predicted by binary population synthesis models before {\it Gaia} and -- if the {\it Gaia} orbits are robust -- likely requires very efficient envelope ejection during common envelope evolution (CEE). To assess the reliability of the {\it Gaia} solutions, we measured multi-epoch radial velocities (RVs) of 31 WD+MS binary candidates with $P_{\rm orb} = (40-300)$ d and \texttt{AstroSpectroSB1} orbital solutions. We jointly fit the RVs and astrometry, allowing us to validate the {\it Gaia} solutions and tighten constraints on component masses. We find a high success rate for the {\it Gaia} solutions, with only 2 out of the 31 systems showing significant discrepancies between their {\it Gaia} orbital solutions and our RVs. Joint fitting of RVs and astrometry allows us to directly constrain the secondary-to-primary flux ratio $S$, and we find $S\lesssim 0.02$ for most objects, confirming the companions are indeed WDs. We tighten constraints on the binaries' eccentricities, finding a median $e\approx 0.1$. These eccentricities are much lower than those of normal MS+MS binaries at similar periods, but much higher than predicted for binaries formed via stable mass transfer. We present MESA single and binary evolution models to explore how the binaries may have formed. The orbits of most binaries in the sample can be produced through CEE that begins when the WD progenitor is an AGB star, corresponding to initial separations of $2-5$ au. Roughly 50\% of all post-common envelope binaries are predicted to have first interacted on the AGB, ending up in wide orbits like these systems.


INTRODUCTION
Close binaries containing a main sequence (MS) star and a white dwarf (WD) are the end products of binary interactions that must have occurred when the WD progenitor was a giant.These systems are therefore important tools to constrain the physics of mass transfer (MT) processes.
The most uncertain mode of binary interaction is common envelope evolution (CEE), in which a companion is engulfed by the envelope of an evolved star and spirals inward on a short timescale.If enough orbital energy is liberated in this process, the envelope is unbound and ejected, leaving behind a post-common envelope binary (PCEB) consisting of the donor's core and the companion star in a close orbit.Otherwise, the result is a stellar merger.CEE is critical in the formation of a wide variety of exotic systems.
In the past few decades, several surveys have been successful in the search for PCEBs (Rebassa-Mansergas et al. 2007;Parsons et al. 2016) and the population they have discovered has played an important role in calibrating CEE models (e.g.Zorotovic et al. 2010).However, in part due to selection biases, the majority of these known PCEBs have short orbital periods of less than a day.Their formation can be relatively well-explained by simply energy conservation arguments (the "α-formalism"; e.g.Livio & Soker 1988;Tout et al. 1997;De Marco et al. 2011), where a fraction of the orbital energy lost during the binary's inspiral is enough to overcome the gravitational binding energy of the giant donor's envelope.arXiv:2405.06020v2[astro-ph.SR] 15 Aug 2024 2 A handful of binaries have been discovered which challenge these simplified CEE models (Wonnacott et al. 1993;Kruse & Agol 2014;Kawahara et al. 2018;Yamaguchi et al. 2024b).These have longer orbital periods of months to a few years, making them too close to have avoided interaction when the WD progenitor was a giant, but much wider than naively expected for orbital energy losses to have been sufficient for successful envelope ejection.If these are in fact PCEBs, they may indicate that additional energy sources are needed to overcome the envelope binding energy and several such sources have been put forward and debated in the literature (e.g.recombination energy Paczyński & Zió lkowski 1968;Webbink 2008;Ivanova 2018;Belloni et al. 2024c,b; jets from the accretor Sabach et al. 2017;Moreno Méndez et al. 2017).
The third data release from the Gaia mission (Gaia Collaboration et al. 2023a) provided orbital solutions for over 169,000 astrometric binaries (Gaia Collaboration et al. 2023b).From these, Shahaf et al. (2024) identified several thousand MS + WD binary candidates, most of which have orbital periods that range between of ∼ 100 − 1000 d and WD masses higher than predicted for binaries formed by stable mass transfer at these periods.This sample tells us that such systems are relatively common and that they were largely missed by previous searches, which were sensitive only to short periods.Therefore, this sample highlights missing pieces in our understanding of binary interactions, motivating followup observations to validate the Gaia solutions and put better constraints on their orbital and stellar parameters.
In this paper, we present follow-up spectroscopic observations of 31 objects from the Shahaf et al. (2024) sample of MS + WD binaries with orbital periods less than 300 days.Section 2 summarizes the original selection carried out by Shahaf et al. (2024) and the properties of our sub-sample.In Section 3, we describe our spectroscopic observations to obtain radial velocities (RVs) and describe fits to the broadband spectral energy distributions (SEDs) to constrain the MS star mass.In Section 4, we jointly fit the RVs and the Gaia astrometry to constrain orbital parameters and flux ratios.We then discuss the resulting flux ratio distribution, inferred success rate of the Gaia astrometric solutions, and comparison of our systems to other populations in Section 5.In Section 6, we present MESA models testing formation channels through CEE as well as stable MT.We briefly discuss the completeness of our sample in Section 7. Lastly, in Section 8, we summarize our main results and conclude.

SAMPLE SELECTION
We are carrying out an RV survey of WD + MS binary candidates selected by Shahaf et al. (2024).Their selection was based on the astrometric mass ratio function (AMRF) which was presented in Shahaf et al. (2019) and defined in terms of observational parameters: where α is the angular semi-major axis of the photocentric orbit, ϖ is the parallax, M 1 is the mass of the primary (i.e. the more luminous star), and P is the orbital period.This can also be written in terms of the mass ratio q = M 2 /M 1 and the flux ratio S = F 2 /F 1 , where F 1 and F 2 are the G-band fluxes of the primary and secondary: Given a relation S(q) for different possible companions (e.g. a single MS star, a tight binary containing two MS stars, or a dark object), one can define regions in A -M 1 space in which each of these possibilities are most likely.In particular, above some critical AMRF for a given M 1 , the system is very likely to host a dark companion/compact object (Shahaf et al. 2023).Shahaf et al. (2024) selected ∼9000 astrometric binaries with AMRF values that are too large to be explained by a single luminous companion.These systems could either host WD companions, or be triples, in which the companion is a tight MS + MS binary.They then removed systems that fall above the main sequence in a color-magnitude diagram, which were suspected to be triples (For a detailed description, we refer readers to Shahaf et al. 2024).
For this paper, we selected 31 objects from this sample with orbital periods less than 300 days to obtain follow-up spectroscopy.Basic information for these are summarized in Table 1.These all have "AstroSpec-troSB1" solutions, meaning that both astrometry and RVs from Gaia were fit with a single orbital solution.The full Shahaf et al. (2024) sample also contains objects with "Orbital" solutions, which are fitted with just astrometry.Some of these are part of a larger sample for which we are also carrying out spectroscopic follow-up and which we will describe in future work.On the left panel of Figure 1, we show the distribution of binaries from the full sample across eccentricity and orbital period.We mark the 31 objects from this work with stars, and those from our larger sample (for which we have not yet obtained sufficient RVs for a well-constrained orbit) with diamonds.We see that these are roughly evenly sampled across orbital periods and cover a wide range of possible eccentricities.On the right panel, we plot the same objects on a Gaia extinction-corrected color-magnitude diagram (CMD)(for ∼ 45, 000 sources from gaiadr3.gaia source with ag gspphot values and parallax over error > 5).The luminous stars are main-sequence G and K stars.They lie on the left half of the main sequence, because Shahaf et al. (2024) discarded redder objects as possible triples.

FEROS
We obtained 135 spectra with the Fiber-fed Extended Range Optical Spectrograph (FEROS; Kaufer et al. 1999) mounted on the MPG/ESO 2.2 m telescope at the La Silla Observatory.All observations used 1×1 binning and the resulting spectra have resolution R ≈ 50, 000.Exposure times ranged from 600 to 1200 seconds, depending on the brightness of the target.The data were reduced with 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 with a second fiber.

GALAH and APOGEE RVs
For a few objects, we also include RVs from the publicly available GALAH DR3 (Buder et al. 2021) and APOGEE DR17 (Abdurro'uf et al. 2022) catalogues.From the GALAH DR3 Main and Value-Added calatalogues, we took the dates, RVs, and RV errors from the columns heliocentricJD, rv galah, and e rv galah, respectively.From the APOGEE DR17 allVisit file, we used the columns JD, VHELIO, and VRELERR.
We obtained FEROS spectra at one epoch for five calibration stars with both GALAH and APOGEE RVs to obtain RV offsets between the three RV scales.We calculate median offsets from FEROS of +0.19 km s −1 and −0.31 km s −1 for GALAH and APOGEE RVs respectively, which we apply when fitting the RVs in Section 4.
All RVs obtained for our objects are listed in Tables 2 to 4 in Appendix A.

SED fitting
We constructed the broadband SED of each source using synthetic SDSS ugriz photometry from Gaia XP spectra (Gaia Collaboration et al. 2023c), 2MASS JHK photometry (Skrutskie et al. 2006), and WISE W 1 W 2 W 3 photometry (Wright et al. 2010).
For declinations below −30 • , we obtained E(B − V ) values from the Lallement et al. (2022) 3D dust map, and for declinations above −30 • , we used the Bayestar2019 (Green et al. 2019) 3D dust map.We take R V = 3.1 assuming a Cardelli et al. (1989) extinction law.We take E(g − r) provided by the Bayestar2019 map to be approximately equal to E(B − V ) (Schlafly & Finkbeiner 2011) and take A 0 (extinction evaluated at 500 nm) provided by the Lallement et al. (2022) map to be A V .
The SED fitting procedure is identical to that described in Yamaguchi et al. (2024b) to which we refer the readers for details, but we summarize it here.We use the MINEsweeper code (Cargile et al. 2020) which was designed to jointly model photometry and spectra, but in this work, we only use the code to interpolate between MIST evolutionary models (Dotter 2016;Choi et al. 2016) and predict synthetic photometry.The free parameters are the parallax, mass M ⋆ , initial metallicity [Fe/H] init , and Equivalent Evolutionary Phase (EEP, a monotonic function of age; see Dotter 2016).The code uses neural network interpolation to swiftly calculate the SEDs and photometry for a given set of parameters.We use emcee, a Python Markov chain Monte Carlo (MCMC) sampler (Foreman-Mackey et al. 2013), to sample from the posterior.For the metallicities, we placed Gaussian constraints based on the values calculated using the Gaia XP very low-resolution spectra as described in Andrae et al. (2023).However, we do not use the values as published in Andrae et al. (2023), which used "single-star" parallaxes from the gaia source table.Instead, we use values re-calculated using parallaxes from the nss two body orbit table which takes into account of the binary motion (Gaia Collaboration et al. 2023b;Halbwachs et al. 2023).We also use a parallax prior based on the Gaia astrometric orbital solution.
We do not fit for the contribution from the WD companion.This is justified as these objects were selected on the basis of their AMRF values which suggested a dark companion.And for a typical WD companion with T eff ∼ 10, 000 K, the expected flux contribution in the optical is small (< 0.1 mag; see discussion in Appendix B of Yamaguchi et al. 2024b).
The results from the SED fitting for two example objects are shown in Figure 2. The fits for all other objects can be found in Appendix B. Note that about half of the objects have GALEX far/near ultraviolet (NUV/FUV) photometry which are plotted on top of the SEDs, but we did not use them in the SED fitting as hot WDs may contribute significantly in the UV.In fact, there are a few NUV and FUV observations where significant excess is present, one being J0947-2018 shown in Figure 2. The medians and standard deviations of the marginalized posterior distributions for each parameter for all 31 objects are given in Table 5.The masses range from 0.64 to 1.09 M ⊙ , with a median at 0.85 M ⊙ .Note that although some of the mass uncertainties are formally very small, we adopt a minimum uncertainty of 0.03 M ⊙ in the later fitting to account for possible systematics in the stellar models and/or photometry (e.g.bias in u-band magnitudes from the Gaia XP spectra; Gaia Collaboration et al. 2023c).

JOINT GAIA + RV FITTING
All of the objects in this sample have Astro-SpectroSB1 solutions with bit indexes of 65535, meaning that the Gaia astrometry and RVs were fitted with a single orbital model containing 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, G are transformation of the Campbell elements (photocenter semi-major axis, inclination, argument of periapsis, and the position angle of the ascending node) which describe the photocenteric orbit from astrometry.The C and H elements describe the motion of the primary as inferred from spectroscopic RVs.
For each object, we simultaneously fit the Gaia astrometric solution and the follow-up RVs.The primary advantage of astrometry over just RVs in fitting orbits is that it constrains the inclination and thus the precise mass of the companion, as opposed to a lower limit.The benefit of having follow-up RVs is that they validate the Gaia solution and, given their small uncertainties, allow us to place much tighter constraints on the parameters.In addition, as we discuss below, RVs allows us to constrain the flux ratio of the components.
We use the likelihood function as defined in El-Badry et al. (2023), which is the product (sum in log space) of two components.The first is the Gaia term.For each model with a vector of the 15 astrometric quantities θ Gaia , the likelihood is calculated as (3) where µ Gaia is the vector of best-fit values from Gaia and Σ Gaia is the corresponding covariance matrix, which can be constructed with the correlation vector provided in the gaiadr3.nsstwo body orbit table.
The second term comes from fitting the RVs.We use a standard Keplerian model and compare its predictions to our observations: Table 1.Basic information from Gaia DR3 of all objects studied 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 arcminutes.G is the G-band mean magnitude, GoF is the Gaia goodness of fit, P orb is the orbital period reported by Gaia, and ϖ is the parallax.Last column lists the number of RVs we measured for each object.
where RV pred (t i ) and RV i are the predicted and measured RVs at a time t i and σ RV,i is the error in RV i .
Our model fits 15 parameters: ra, dec, parallax ϖ, pmra, pmdec, period P orb , eccentricity e, inclination i, argument of periapsis ω, angle of the ascending node Ω, periastron time T p , center of mass velocity γ, luminous star mass M 1 , WD mass M 2 , and finally, the flux ratio S. We place a Gaussian prior on M 1 using the best-fit value and uncertainty obtained from the SED fitting (with an uncertainty floor of 0.03 M ⊙ ; Section 3.3).We also run a model using just the Gaia constraints (equation 3) to assess the relative importance of the Gaia data and our RVs in constraining the orbit.
The flux ratio S comes in because the angular semimajor axis of the photocenter α is related to that of the primary α 1 by By definition, S should be 0 or greater.It is zero in the case of a completely dark companion, where α = α 1 .We emphasize that astrometry alone does not directly constrain this S -RV measurements of the luminous primary are necessary.Since all objects presented in this work have combined astrometric and single-lined spectroscopic solutions (AstroSpectroSB1), the Gaia solutions alone provide some constraint on this value (through the C and H elements) which will be tightened by the follow-up RVs.In the case of purely astrometric solutions (Orbital) however, follow-up RVs are necessary.
As described in Section 2, our objects were selected to have AMRF values and CMD positions indicative of dark companions.We therefore expect S ≪ 1 if the Gaia solutions are robust (this is true, even for those with Orbital solutions).This means that if a lower bound of zero is set, the posterior distributions will tend to run into this lower bound, and since S cannot be negative, this will lead to a bias toward higher M 2 .Note also that if the Gaia solution is exactly correct, S = 0 should give the minimum predicted RV variation (because in the absence of noise or systematic errors, the photocenter follows the primary).But if the Gaia solution is not exactly correct, it is possible for the true RV variation to be smaller than it predicts, which corresponds to negative S. To avoid introducing a bias toward higher M 2 , in the fitting, we allow S to range from -0.5 to 1 (none of the solutions run into the lower bound).
The left panel of Figure 3 compares RV curves from 50 random draws of the posterior using just the Gaia solution (cyan) to those from jointly fitting the RVs and the Gaia solution (black), for two example objects, J0602-1747 and J1834+1525.These objects were chosen to represent a typical case with a reliable Gaia orbital solution that agrees with our measured RVs, and a case with an unreliable Gaia solution.Similar plots of the RV curves for all 31 objects can be found in Appendix C. The right panel shows the corner plots for several key parameters.For J1951-0305, we see that the inclusion of our follow-up RVs tightens the constraints on the parameters (particularly P orb and M 2 ) but the best-fit values from both fits are consistent with each other to within 1σ.Meanwhile, for J1834+1525, we see that there are significant disagreements between the solutions from the two fits.As we discuss further in Section 5.2, this is one of the few objects for which the Gaia solutions have underestimated errors, despite having goodness of fit values indicative of an unproblematic solution.

Flux ratio
Figure 5 shows a histogram of the distribution of S. As expected, the majority of objects have S ≲ 0 which is indicative of a dark companion/WD.However, this Figure 3. Results of joint Gaia + RV fitting for two example objects, J0602-1747 (top) and J1834+1525 (bottom).These represent cases with reliable and unreliable Gaia orbital solutions, respectively.Left: RV curves using orbital parameters from 50 random draws of the posterior resulting from the fitting of the Gaia solution only (cyan) and the joint fitting with our follow-up RVs (black).Right: The corresponding corner plots showing some key parameters for the same objects.
is not conclusive, because depending on the masses of both components, a given flux ratio may or may not be large enough to be consistent with a triple scenario.
Equating equations 1 and 2, one can solve for the secondary mass M 2 as a function of the flux ratio S, given all of the other parameters from the orbital fitting (Section 4).On Figure 4, we plot this curve as a black line for several exemplar objects.The best-fit values of M 2 and S must be a point that lies somewhere on this line (black marker).
We use synthetic photometry of zero-age MS stars in the Gaia band from MIST to compute theoretical curves on this plot for a single MS companion (blue dotted) and an inner binary (orange dashed) of two equal-mass MS stars.
Based on this plot, we can classify each object into one of two groups: • Group 1: Inconsistent with being a triple system.
The best-fit point (black marker) and its error bars (1-σ) do not cross the theoretical curve for an inner binary companion (orange dashed line).
• Group 2: Cannot rule out a triple.All other cases -the best-fit point is consistent with the orange dashed line, to within 1-σ.Note that objects in group 2 could still host WD companions -the data just do not rule out a triple.
In Figure 4, we show examples of objects that belong to each of these two groups.In total, 26 out of the 31 objects are in group 1.This means the large majority of all objects are inconsistent with being in triples and are thus likely to host WDs.As for the remaining 5 objects, while we cannot rule out triples, we note that all of them are also consistent with hosting dark companions.
Modeling of the SEDs with triple models could likely set more stringent limits on S for objects with only a few RVs, but here we only include constraints based on the joint fitting of the Gaia solutions and RVs.
The right panel of Figure 5 plots the eccentricity against S for all objects, with the markers indicating their "groups" as described above.We notice two key features: (1) all eccentricities are relatively low.This is expected for binaries hosting WD companions, but unexpected if the systems are primarily hierarchical triples (note also that those in group 2 tend to have larger un-certainties in S); (2) even among the systems in group 1, most of the eccentricities are significantly different from 0. This indicates any tidal circularization process has not been very efficient.We compare these eccentricities to those of millisecond pulsar + WD binaries in Section 5.4.

Success rate of Gaia solutions
One purpose of obtaining RV follow-up was to validate the Gaia astrometric solutions -in this case, specifically those with AstroSpectroSB1 solutions.In Figure 6, we plot the best-fit values of e and M 2 obtained from the fitting with and without the inclusion of RVs.If they are in agreement, the points would lie along the gray dashed line.We see that for the majority of objects, the change in the best-fit values are minimal, with the primary effect of the RVs being to reduce the uncertainties.However, there are two (five) cases with > 2σ (> 1σ) discrepancy in e between the two solutions.The two sources are J1834+1525 and J1922-4624.Figures 19 and 20 in Appendix C show that the discrepancy between measured and predicted RVs is clearly visible.These objects have goodness-of-fit values of 1.7 and 2.7 respectively, which are not enhanced compared to the rest of the objects.These sources would therefore not have been identified to have spurious solutions without RV follow-up.
Thus, at least 2 out of the 31 objects in our sample ended up having unreliable Gaia solutions.While this suggests that most AstroSpectroSB1 solutions are reliable, there is a non-negligible amount which are not, so RV follow-up is important for objects that appear unusual based on their Gaia-only solutions.Note also that we pre-selected sources that did not have high goodnessof-fit values (≲ 4.8).The fraction of spurious solutions also appears to be higher for purely astrometric solutions, as we will describe in future work.

Comparison to literature PCEBs
In Figure 7, we plot orbital period, P orb , against WD mass, M 2 , for the objects in our sample as well as literature PCEBs.The colors of the points represent the mass of the luminous companion, M 1 .We also mark the corresponding minimum orbital separation, a peri = a(1−e), on the right axis for a 0.6M ⊙ WD and a 1M ⊙ companion in a circular orbit (this is the default assumption that we make for any conversion between P orb and a peri ).Binaries formed through stable Roche lobe overflow are expected to evolve along the gray dashed line derived from the P orb −M 2 relation from Rappaport et al. (1995), with the shaded region showing the spread in P orb by a factor of 2.4 (as a result of different giant masses and compositions).The blue solid line marks the maximum radius reached by an isolated WD progenitor.We obtained this using MIST evolutionary tracks for stars with a range of masses (Dotter 2016;Choi et al. 2016) which were converted to WD masses using the initial-final mass relation from Williams et al. (2009).Any system lying below this line must have interacted at some point during its evolution.
The best-studied sample of literature PCEBs comes from the Sloan Digital Sky Survey (SDSS) which identified many MS + WD binaries.From these, Rebassa-Mansergas et al. (2007) discovered 37 new PCEBs.On Figure 7, we plot these, along with 25 previously known PCEBs, for a total 62 systems (compiled by Zorotovic et al. 2010).We also plot the PCEBs from the "white dwarf binary pathways survey" with more massive AFGK companions (Hernandez et al. 2021(Hernandez et al. , 2022a,b),b).All of these systems have orbital periods of ≲ 1 d, and their formation can be explained by simple energy conservation arguments, formally known as the α-formalism.In the simplest case, this is the balance between a fraction α of the orbital energy loss and the gravitational binding energy of the donor's envelope.Various models of these systems have found that a value of α ∼ 0.3 can reproduce the observed population (Zorotovic et al. 2010;Davis et al. 2010;Toonen & Nelemans 2013;Camacho et al. 2014;Zorotovic & Schreiber 2022;Scherbak & Fuller 2023).
The above systems are labelled "short-period PCEBs" in Figure 7 to distinguish them from the "wide PCEBs" that have been found over the years.These have much longer orbital periods of ∼ months, and the systems that have been characterized so far host unusually massive (≳ 1.2 M ⊙ ) WDs.These include IK Peg (Wonnacott et al. 1993) as well as five systems from Gaia recently identified by Yamaguchi et al. (2024b).Several self-lensing binaries (SLBs) discovered with Kepler are also plotted.These host less massive WDs and have even longer orbital periods from ∼ 100 − 700 d (Kruse & Agol 2014;Kawahara et al. 2018;Masuda et al. 2019;Yamaguchi et al. 2024a).All of these systems (with the possible exception of one SLB) lie firmly below both the blue solid line and the gray region, meaning that they must have interacted at some point but not via stable mass transfer, pointing towards CEE.However, the wide orbits of these systems imply less orbital shrinkage during the common envelope phase than is released in the formation of close PCEBs.This has to lead to works arguing that additional sources of energy are required to unbind the envelope (e.g.Davis et al. 2010;Zorotovic et al. 2010;Parsons et al. 2016;Scherbak & Fuller 2023;Yamaguchi et al. 2024b), while others circumvent this with mass transfer from a donor on the TP-AGB  1) calculated using the best-fit orbital parameters derived in Section 4; black point marks the best-fit values of M2 and S. The blue solid line and orange dashed lines correspond to theoretical curves for single MS star and equal-mass inner binary companions.On the top row are "group 1" objects (as described in Section 5.1) which must host dark companions: the constraints on M2 and S are inconsistent with an inner MS+MS binary for any S. Bottom row shows "group 2" objects, for which the flux ratio constraint cannot rule out an inner MS+MS binary.
The objects from this work are plotted with star markers.These have orbital periods ranging from ∼ 50 − 300 d, intermediate between the wide PCEBs and Kepler SLBs.They all lie in a region of this parameter space expected to be sparsely populated according to standard binary evolution models (Shahaf et al. 2024) -having more massive WDs than expected at their orbital periods if they evolved through stable MT, while also having significantly longer orbital periods than expected for PCEBs.We discuss possible formation scenarios of these systems in Section 6.

Period-eccentricity relation
On Figure 8, we plot the period against eccentricity of our objects, along with literature millisecond pulsars (MSPs) + WD binaries from the Australia Telescope National Facility (ATNF) pulsar catalogue (Manchester et al. 2005; we plot the sample compiled by Hui et al. 2018).We distinguish systems with WD masses above and below 0.45 M ⊙ .This is the approximate boundary between He WDs and more massive CO or ONeMg WDs, which correspond to systems thought to have formed from stable and unstable MT.The gray line is the theoretical relation expected for systems containing He WDs formed from stable MT derived in Phinney 0.08 0.04 0.00 0.04 0.08 0.12 0.16 0.20 0  (1992).Cohen et al. (2024) recently extended this relation to CO WDs with AGB progenitors and found an upper limit in the eccentricity of ∼ 3 × 10 −3 , which is marked with a dashed horizontal line.We see that our objects, plotted with pink stars, all have significantly higher eccentricities than the vast majority of the MSP + WD binaries at similar periods.We also plot the 5 systems from Yamaguchi et al. (2024b) with green diamonds, four of which lie at comparable eccentricities at shorter orbital periods.None of these MS + WD binaries have eccentricities consistent with zero nor are they consistent with the Phinney (1992) relation.All of them also lie above the Cohen et al. ( 2024) limit.Since it is unlikely that all of these systems are in triples, this suggests that Kozai-Lidov oscillations are unlikely to be responsible for their eccentricities.Moreover, since many MSP + WD binaries -which must have undergone stable MT in order to spin up the MSP -do lie close to the Phinney (1992) relation, this suggests that these systems are not stable MT products.We note also that several recent 1D and 3D models have shown that initially eccentric systems that undergo CEE can retain eccentricities of ∼ 0.01−0.1 post-inspiral, comparable to those of our systems (Sand et al. 2020;Glanz & Perets 2021;Bronner et al. 2024).

Formation through CEE
Using Modules for Experiments in Stellar Astrophysics (MESA; Paxton et al. 2011Paxton et al. , 2013Paxton et al. , 2015Paxton et al. , 2018) models of the WD progenitor, we can calculate the range of orbital periods that can result from CEE.We performed similar calculations in Yamaguchi et al. (2024b) and refer the reader to that work for more details.In short, we use a single star model of the WD progenitor (i.e. the donor) to trace its properties with time and use the α-formalism to predict the final orbital separation of the PCEB as a function of the separation at the time when mass transfer began.We use MESA version 22.05.1.
The α-formalism is a statement of energy conservation between the binding energy of the star's envelope E bind and a fraction α CE of the orbital energy released during CEE: where M WD is the WD mass, M i is the WD progenitor mass, M ⋆ is the companion mass, a i is the initial separation at the onset of mass transfer (calculated using the  Phinney (1992).We see that our objects have significantly higher eccentricities than predicted by this relation.
tant PCEB).Given the short timescale of the CEE, we assume M ⋆ remains constant.
We define E bind as the sum of the gravitational binding energy of the envelope E grav and a fraction α int of its internal energy E int .The internal energy makes the total binding energy of the envelope less negative, meaning it is less bound and thus easier to eject (see e.g.Paczyński & Zió lkowski 1968;Ivanova et al. 2013).This term was introduced to allow for the formation of wide PCEBs (P orb ∼ months -years) where the orbital shrinkage during the CEE may not be enough to overcome the gravitational binding energy (e.g.Scherbak & Fuller 2023;Yamaguchi et al. 2024b).MESA defines the internal energy as the sum of the "thermodynamic" (thermal + radiation energy) and recombination energy (Paxton et al. 2018).It does not provide an easy way to extract the two components separately, so we estimate the specific thermal energy as 3P/(2ρ), and radiation energy as 3P/ρ.The total binding energy of the enve-lope is: (7) where we separately define the efficiencies α th and α rec corresponding to thermodynamic and recombination energies respectively.In the following calculations, we always take α th = 1.
For different combinations of α CE and α rec , and taking observed values of M ⋆ , we can rearrange equation 6 to solve for the final separation a f as a function of the donor model parameters when mass transfer starts.For all of our calculations, we use M ⋆ = 0.85 M ⊙ which is the median mass of luminous companion in our systems.
The WD masses for the objects in our sample range from around 0.5 to 0.8 M ⊙ (excluding a single object J2006+1234 at 1 M ⊙ ), with a median at ∼ 0.6 M ⊙ which corresponds to progenitor masses of 1 to 3 M ⊙ , with a median at ∼ 1.5 M ⊙ (e.g.Williams et al. 2009;Cummings et al. 2018;El-Badry et al. 2018).Therefore, we ran 1, 1.5, and 3 M ⊙ MESA models.
The inlists are based on those from the 1M pre ms to wd test suite, though several changes and additions have been made following (Rees & Izzard 2024).
Recently, Belloni et al. (2024a,b,c) emphasized the importance of mass transfer from a donor in the thermally pulsating AGB (TP-AGB) phase to produce wide PCEB orbits.This phase is difficult to model, with several instabilities that arise as the envelope mass is reduced by mass loss, causing it to expand and cool.These can terminate models prematurely or lead to prohibitively slow convergence, and are especially problematic for more massive stars which are the progenitors to ultra-massive WDs.In Yamaguchi et al. (2024b), we terminated models at the onset of the TP-AGB phase for this reason.In this work, we run our models until the end of the TP-AGB, with help from inlists and custom MESA routines created by Rees & Izzard (2024) which prevent instabilities and ensure convergence (we refer readers to their paper for details; their models are also discussed in Rees et al. 2024).Note that in their work, they use the mass loss scheme of Vassiliadis & Wood (1993) on the AGB, as opposed to that of Bloecker (1995).We briefly discuss the effect of the chosen wind prescription in Section 6.1.2.Models are terminated when most of the star's en-velope has been stripped off (envelope mass limit = 0.01) or if the minimum timestep (min timestep limit = 1d-6) is reached.The latter only occurs for the 3 M ⊙ model, which has an envelope mass of 0.28 M ⊙ at termination.
6.1.1.Definition of the envelope in the calculation of binding energy In calculating the binding energy of the envelope, it has been standard practice in the literature to consider the entire envelope, from the boundary of the donor's core to the photosphere (e.g.Marchant et al. 2021;Scherbak & Fuller 2023;Yamaguchi et al. 2024b;Belloni et al. 2024b).However, in the formation of a wide PCEB, the companion never reaches the vicinity of the donor's core.It is thus not clear that the donor's orbital energy plays any significant role in ejecting material in the inner envelope.This means that rather than integrating energies from the core, it may be sufficient to integrate from the observed orbital separation to the surface.In practice, the star will likely expand when its outer envelope is stripped off, lifting more material from the inner envelope to the orbit of the secondary.However, in this case, most of the work to eject the inner envelope is done by the core, not the companion.
The binding energy of the outer envelope can be significantly less than that of the whole envelope.Moreover, the relative importance of the thermal and recombination energies varies inside the envelope, with recombination energy dominating in the outer layers and thermal energy dominating deeper in the interior.In the calculations below, we trace the binding energies of both (a) the entire envelope above the helium core, and (b) only the material beyond 0.7 AU, which is the typical final separation of the binaries in our sample.

Results of CEE models
Starting with the 1.0 M ⊙ model, the leftmost panel of Figure 9 shows the evolutionary track of this star on a Hertzsprung-Russel (HR) diagram.On the middle columns, we plot a f against a i with and without the inclusion of recombination energy, and for various values of α CE .We show the full tracks but in practice, only the sections with the thicker line weight is accessible, where a i (or R) is larger than it was at any earlier time.Otherwise, CEE would have been triggered prior to reaching that point in the donor's evolution.The red dashed lines at 0.01 AU mark the approximate separation below which the MS companion would not fit in the orbit.The gray and blue dotted lines at 0.15 and 0.7 AU respectively are the typical separations of previously discovered wide PCEBs including IK Peg and the systems from Yamaguchi et al. (2024b), and the 31 objects from ), but where no recombination energy is included (αrec = 0).The three solid lines are cases for αCE = 0.3, 0.6 and 0.9, with the inclusion of 50% of the available recombination energy (αrec = 0.5).The horizontal dashed lines at 0.01, 0.15, and 0.7 AU represent the minimum separation below which the luminous primary would not fit in the orbit, the minimum separation of the "wide PCEBs", and the typical separation of the objects from this work.Black arrows indicate initial separations where the addition of recombination energy has resulted in the overall binding energy being positive, where the envelope is formally unbound.The upper and lower rows correspond to cases where the binding energy is calculated over the entire envelope vs. just above the final orbital separation of 0.7 AU, respectively.Right: Zoom-in of the region above 0.15 AU of the middle plot in the same row, isolating the regions with thicker line weight which can be reached by the donor in practice (i.e.where the radius is larger than at any previous point in the star's evolution).Wide orbits can form without the addition of any recombination energy if αCE is sufficiently large.This is true with both definitions of the envelope binding energy, though the range of initial separations over which they form significantly increased when just the outer envelope is considered.
this work, respectively.The right columns zoom in on the parameter space that produces wide PCEBs.
First, we examine the result of considering the binding energy of the entire envelope (upper row).We find that for a 1 M ⊙ donor, the orbits of wide PCEBs (both ≳ 0.15 AU and ≳ 0.7 AU) can form without the need to invoke any recombination energy (α rec = 0; orange dashed line).This is true even for relatively modest values of α CE ∼ 0.4.As seen on the rightmost panel, even wider orbital separations can result with the inclusion of 50% of the available recombination energy.
Next, we move on to the same 1 M ⊙ model but considering the binding energy of the envelope above 0.7 AU (lower row).We see that the results change significantly.Small initial separations corresponding to the donor radii being less than ∼ 0.7 AU are, of course, ex-cluded.We see that in this case, even larger separations of ∼ 2 AU can be reached with α rec = 0, for a wider range of initial separations.
The situation changes slightly for the 1.5 M ⊙ case, shown in Figure 10.Once again considering the binding energy of the entire envelope first (upper row), we find that for α rec = 0, the maximum separation predicted by our models is just below ∼ 0.7 AU for α CE = 1.We note however that the exact results are dependent on model assumptions (e.g. using the Vassiliadis & Wood (1993) wind prescription results in a small range of a i that can produce a f > 0.7 AU).As expected, with α rec = 0.5, wider orbits can form.For the case of the binding energy above 0.7 AU (lower row), there is a small range of initial separations resulting in a f ≳ 0.7 AU (and a much larger region for a f ≳ 0.15 AU), even with α rec = 0.  9, but for a 1.5 M⊙ model.For this model, we see that wide orbits like those of our systems (a f > 0.7 AU) do not form without recombination energy (orange dashed line) if the binding energy of the entire envelope is considered.Meanwhile, there is a small range of initial separations over which they can form if only the outer envelope is considered.
Lastly, for the 3 M ⊙ model (Figure 11), wide orbits (a f > 0.7 AU) can form without added recombination energy for α CE ≳ 0.4, similar to the 1 M ⊙ model.This suggests that there is a non-monotonic relationship between the donor mass and post-common envelope orbital separation, which may be worth exploring with finer model grids.Unlike the 1 M ⊙ case, the consideration of just the outer envelope in the 3 M ⊙ model does not significantly change the resulting final orbits nor the range of initial separations over which wide orbits can form, which reflects the change in envelope structure between the two models.
In these calculations, we exclude regions where the total binding energy becomes positive due to the addition of recombination energy (indicated by the black arrows in the figures), where the α-formalism is no longer appropriate.In these regions, the envelope is formally unbound and therefore, even wider orbits with little to no shrinkage may be produced.
In Figure 12, we plot the envelope mass against initial orbital separation for the three models.Comparing to Figures 9, 10, and 11, we see that wide orbits with a f ≳ 0.7 AU are formed predominantly when mass transfer begins during the thermal pulses.This corre-sponds to envelope masses below roughly 0.2, 0.4, and 0.8 M ⊙ for the 1, 1.5, and 3 M ⊙ models respectively.In panel h of Figure 13, we see that for the 1.5 M ⊙ model, this point is reached during the last 10 5 years.This highlights the importance of modeling the TP-AGB phase in its entirety, beyond the first thermal pulse and until a significant amount of the envelope is lost.
The predicted core masses corresponding to wide final orbits range from ∼ 0.55−0.65 M ⊙ across the 1, 1.5, and 3 M ⊙ models.These values are similar to the masses of most of the WDs in our sample.6.1.3.The importance of thermal pulses to forming wide PCEBs Belloni et al. (2024a,b) ran a grid of MESA models to explain the formation of KOI-3278 (Kruse & Agol 2014) and FN Sgr (Brandi et al. 2005) and concluded that such systems can form through CEE if the donor is on the TP-AGB without the need for any recombination energy.Despite some differences between their MESA models and ours (single star vs. binary models, various parameter choices), the final orbits are predicted using the α-formalism and their results are in reasonably good agreement with ours.  .Envelope mass against separation at initial Roche lobe overflow for the 1, 1.5, and 3 M⊙ models.As in Figure 9, regions where ai is larger than at any previous point have been plotted with a thicker line weight.The black arrow indicates the direction of evolution.We see that the range of initial separations needed to form the orbits of our systems in each of these cases (see Figures 9,10,and 11) occur primarily during the thermal pulses.
In Figure 13, we plot several stellar properties against time for the 1.5 M ⊙ model on the TP-AGB (the energies plotted are of the entire envelope).From panels a -c, we see that as expected, the internal energy increases with the addition of recombination energy, and correspondingly, the total BE of the envelope becomes less negative.Unlike the models in Belloni et al. (2024a), we find that the total binding energy of our star without recombination energy approaches zero but always remains negative.Note also that the contribution of the radiation pressure to the total thermodynamic energy is sub-dominant to that of the gas pressure throughout the star's evolution until the final few thermal pulses (panel d).The predicted evolution is similar for the other models.
In panel (e), we plot the structural binding energy parameter, λ, defined as where M d is the total mass of the donor, M d,env is its envelope mass, and R d is its radius.de Kool (1990) first introduced λ to account for differences in envelope structure in calculating the binding energy, E bind (generalizing the approximation from Webbink 1984).Tabulated values of λ have been widely used in the context of the αformalism to estimate E bind without the use of detailed stellar models like MESA.Without including recombination energy and in between thermal pulses, we see that λ takes typically assumed values of ∼ 0.5.However, during pulses, the envelope is less tightly bound and it increases to ∼ 1.With recombination energy included, E bind can change signs, which causes λ to blow up and then become negative as E bind becomes positive during thermal pulses.Even between pulses, λ can get as large as ∼ 5.This illustrates the very large range of possible λ values for TP-AGB stars, building on previous works on the dependence of λ and a star's evolutionary phase (e.g.Dewi & Tauris 2000;Tauris & Dewi 2001;Xu & Li 2010).
On the AGB, CEE will almost always be initiated during a thermal pulse, since it is during thermal pulses that stars reach their largest radii (panel f) and will first encounter their companions (and have the least bound envelopes; panels c and e).Similar conclusions were made by González-Bolívar et al. (2022) who carried out 3D smoothed particle hydrodynamic simulations of CEE from an AGB star during a thermal pulse to a low mass companion.It is critical that future stellar model grids used in population synthesis codes capture the properties of AGB stars during the thermal pulses, not between them where they spend most of their time.
We emphasize that while recombination is not required to explain wide PCEB orbits based on energy conservation arguments, it dominates the internal energy in the outer envelope of a TP-AGB donor (on average, E rec /E int ∼ 0.8 integrated over R > 0.7 AU).Therefore, recombination likely does play a role in the expansion and ejection of the envelope and its effect should not be neglected in making predictions about the orbits of PCEBs.Indeed, it has been shown to be important in more detailed models of the CE phase itself, both in 1D and 3D (e.g.Sand et al. 2020;González-Bolívar et al. 2022;Bronner et al. 2024).

Expected PCEB period distribution
Figure 14 shows the distribution of the orbital period and separation of solar-type MS + MS binaries (Raghavan et al. 2010).Systems in the red region of the plot (a i ∼ 2 − 6 AU) will undergo the first common envelope event when the donor (i.e. the more massive MS star) becomes an AGB star.As we showed with our MESA models, these can form wide PCEBs (a f ≳ 0.15 AU).Those in tighter orbits (a i ∼ 0.2−2 AU), shaded in blue, will interact as the donor goes on the first giant branch and will likely form close PCEBs.Those in even smaller orbits will likely merge while those in orbits ≳ 6 AU will likely never interact.The red and blue regions are comparable in area (∼ 10 %), suggesting that we expect similar number of PCEBs in close (P orb ≲ 1 d) and wide (P orb ∼ 10 − 1000d) orbits.We emphasize that this is an over-simplification -we see from Figure 7 that there are several short-period PCEBs with massive WDs that likely formed from AGB progenitors.Thus, not all systems in the red region of Figure 14 will form wide orbits.One possibility is that CEE only leads to wide orbits if interaction begins during the most extended phase of a thermal pulse.Despite this caveat, our analysis suggests that a significant fraction of binaries that begin mass transfer when the donor is on the AGB produce wide PCEBs, and thus that wide and close PCEBs may form with comparable birth rates.

Formation through stable MT
Given the very wide orbits of our systems compared to previously identified PCEBs, we also consider the possibility of formation through stable mass transfer, as suggested by Hallakoun et al. (2023).There are several different tests for this possibility which we describe briefly here.
The first one, which has already been discussed, is whether or not systems fall on the theoretical P orb − M 2 relation along which systems undergoing stable Roche lobe overflow are expected to evolve.We plot this relation from Rappaport et al. (1995) on Figure 7.We note that this specific relation was derived for giants with masses between 0.8 to 2 M ⊙ , which have degenerate cores on the RGB, and thus should be appropriate for the majority of our objects.The fact that none of our objects lie within the uncertainties of the relation suggests that their orbits are inconsistent with having formed from stable mass transfer.
A different commonly used criterion of mass transfer stability is whether mass transfer is expected to selfregulate or lead to a runaway.This boundary is often defined by a critical mass ratio, q c .Using the initial mass -final mass relation for WDs, we can calculate that the initial mass ratios of our systems (q = M donor /M accretor ) range from ∼ 1 − 7, with a median value of ∼ 2 − 3. Note that here we are taking the mass of the accretor to be equal to the mass of the luminous companion today, making these values of q upper limits.Values of q c above which dynamically mass transfer occur vary between works, with typical values of ∼ 1 to the more conservative values of ∼ 3−4 (e.g.Hjellming & Webbink 1987;Chen & Han 2008;Ge et al. 2010Ge et al. , 2015;;Temmink et al. 2023).Therefore, the single object with q ∼ 7 lies firmly in the unstable regime.Meanwhile, the others have q within the range of literature q c values.
To contextualize these mass ratios, we perform another experiment which involves comparing the donor's response to mass loss, Figure 13.Stellar properties as a function of age for a 1.5 M⊙ star on the TP-AGB.On the left column, we plot the gravitational binding energy (a), internal energy with and without recombination energy (b), total binding energy (c; Egrav,env + Eint,env), and ratio of radiation to thermal energy (d).All energies are of the entire envelope.In panel (c), points where the total binding energy becomes positive (i.e. the envelope is formally unbound) are plotted with circle markers -this only occurs if the recombination energy is considered.On the right column, the structural binding energy parameter λ (e; see text), radius (f), luminosity (g), and envelope mass (h).In blue is the region corresponding to systems where the first interaction will occur when one star becomes a red giant.If such binaries undergo CEE, these are expected to form close PCEBs.In red are systems where the first interaction will occur when one star is on the AGB, which is expected to form wide PCEBs.These systems each make up ∼ 10% of the total population.man et al. 1997).In the case where the donor expands faster than the Roche lobe, the mass transfer becomes increasingly more rapid, and thus unstable.To this end, we take the 1.5 M ⊙ model evolved regularly up to a certain point (Section 6.1) and from there, remove its mass by increasing the wind strength (setting Reimers scaling factor = 5.0 and Blocker scaling factor = 0.5).This corresponds to mass loss rates ranging from ∼ 10 −10 to 10 −6 M ⊙ yr −1 .We start mass loss at the start of the sub-giant branch (SGB), roughly halfway up the RGB, and the start of the TP-AGB.The left panel of Figure 15 shows the change in the radius as a function of mass, from which we can get ζ d .We can then compare this to ζ L which is just a function of the mass ratio q and the parameter α (not to be confused with α CE ) which quantifies the mass lost as an isotropic wind from the accretor (along with several other parameters for other modes of MT; Soberman et al. 1997).On the right panel, we plot ζ d and ζ L as a function of mass ratio where the mass of the accretor is once again taken to be 0.85 M ⊙ and assumed to be constant (i.e.fully non-conservative, β = 1).We see that for MT that begins on the SGB and RGB, ζ d < 0 and ζ L > 0 for most values of q (i.e. the donor expands as the donor loses mass while the Roche lobe shrinks).Therefore, we expect mass transfer occurring at these mass ratios to become unstable and remain that way.The behaviour for MT beginning on the TP-AGB is complex due to the occurrence of thermal pulses, making this stability criterion not well suited for this case.We note also that the donor's response can vary depending on the mass loss rate (e.g.Passy et al. 2012;Ge et al. 2020) so the outcomes of this analysis do not necessarily encompass all possible scenarios.Lastly, we run several binary models of a 1.5 M ⊙ donor and a 0.85 M ⊙ point-mass companion.We take an isolated model of the donor star, as described in Section 6.1, load it into a binary model at a wide enough separation for there to be no initial interaction, and let it evolve until it overfills its Roche lobe and mass transfer begins.For the MT scheme, we use the optically thick overflow from Kolb & Ritter (1990) and assume that it is fully non-conservative (mass transfer beta = 1; Rappaport et al. 1982;Tauris & van den Heuvel 2006).Models are terminated when the envelope is stripped off.In Figure 16, we plot the orbital separation, mass loss rate, and Roche lobe filling factor as a function of time since the start of mass transfer from a donor on the SGB, RGB, and TP-AGB (left, middle, and right columns, respectively).We see that for the SGB donor, the mass transfer remains stable -there is an initial rise in the mass loss rate which peaks at ∼ 10 −8 M ⊙ yr −1 and then gradually declines, and the donor is contained within its Roche lobe throughout (i.e.filling factor ≲ 1).Therefore, at these mass ratios, it is possible to have stable MT if the donor is on the SGB, though this would not produce WDs with masses that we observe.Meanwhile, for the RGB and TP-AGB donor, the donor's envelope is fully convective, and we see that the mass loss rates peak at much larger values of ∼ 10 −2 M ⊙ yr −1 and the filling factor exceeds 1.These point to the onset of instability.While the models do not terminate, the subsequent evolution is likely unreliable -specifically, the accretor will likely be unable to accept mass at the rate it is being delivered and thermally adjust while remaining in equilibrium, so it will expand itself, leading to CEE.Moreover, once a star greatly fills its Roche lobe, mass loss could occur through the outer Langrangian points (L2/L3) which can carry away more angular momentum and promote unstable mass transfer (e.g.Shu et al. 1979;Reichardt et al. 2019;Marchant et al. 2021).A caveat to this analysis is that we did not explore the different parameters to tune mass transfer efficiency (α, β, δ, γ; see Soberman et al. 1997;Tauris & van den Heuvel 2006).For simplicity, we assumed that all mass is lost from the vicinity of the accreting star (i.e.β = 1) and thus we cannot exclude the possibility of stable nonconservative mass transfer in all cases.This may be a topic of interest to explore in future work.

Takeaway from the MESA models
Figure 15.Left: Change in the radius of an initially 1.5 M⊙ donor as mass is stripped away from it at different points in its evolution.The black arrow indicates the direction of evolution.We see that up until the final stages of envelope stripping, its radius consistently increases as it loses mass (i.e.ζ d < 0).Right: ζ d against mass ratio, assuming Ma = 0.85 M⊙.We also plot curves above which there are conservative and fully non-conservative (β = 1) stable mass transfer (Soberman et al. 1997).We see that our models lie firmly below these lines when M d /Ma > 1, suggesting that mass transfer will be unstable at its onset.
In summary, these considerations suggest that the mass transfer in the progenitors of our WD+MS binaries likely would have become unstable, making them wide PCEBs (though note the caveat above about nonconservative mass transfer).This is in contrast with the predictions of most population synthesis models, which predict no binaries to exist with the periods and masses characteristic of the sample (e.g.Shahaf et al. 2024).However, detailed stellar models with realistic envelope structures like those of this work and Belloni et al. (2024b) suggest that it is possible to produce these PCEBs as long as mass transfer begins when the donor is on the TP-AGB where its envelope is very loosely bound.What is still needed is a population model to test whether or not the sheer abundance of the observed systems in this region of the parameter space matches what is expected, given these stellar models and a realistic initial binary population.
We also mention that there are related alternative theories that could produce wider orbit systems that we have not discussed in this work, with one being the grazing envelope evolution (GEE) where jets from the accretor prevents its complete engulfment in a common envelope (for details, see Soker 2015).GEE has been proposed to explain post-AGB binaries in intermediate orbits ∼ 1 AU (e.g.Abu-Backer et al. 2018) comparable to those of the systems in this work.

COMPLETENESS OF THE SAMPLE
Looking at Figure 8, there is a lack of systems with lower WD masses that lie in the stable MT region.This is likely at least partly a selection effect.The systems were selected from the Shahaf et al. (2024) sample of probable MS + WD binaries based on their AMRF, A, values.Given a primary mass, this requires that they contain a WD that is massive enough such that no single MS companion could produce as large an AMRF as observed.This boundary is shown in Thus, less massive WDs can only be detected in systems that host less massive primaries.However, as described in (Shahaf et al. 2024), the sample is biased against the least massive MS companions because they are dimmer and and therefore would not be the photocentric primary if they host a sufficiently hot WD.This is similarly discussed by Garbutt et al. (2024) whose sample of WD + FGK binaries were selected on the basis of their UV excess, making them most sensitive to lowmass WD companions and complementary to this sample.Moreover, as described in Halbwachs et al. (2023), AstroSpectroSB1 solutions include RV measurements which are only made for bright objects (G RVS < 12).These factors combined means that at these periods (∼ 100 − 1000 d), there is a bias against finding stable MT systems with the expected WD masses.It may be possible to identify such systems in future work on the basis of their near-zero eccentriciities.
A detailed discussion of the selection function will be left for future work.Here we simply emphasize that Figure 16.Orbital separation, mass loss rate, and Roche lobe filling factor as a function of time for binary MESA models where mass transfer from a 1.5 M⊙ donor begins on the SGB, RGB and TP-AGB.The start of the Roche lobe overflow (RLOF) is defined to be the point where the filling factor, R1/RL, first reaches 1 (gray dashed line).We see that mass transfer from a donor on the SGB remains stable, while in the other two cases, instability seems to be triggered (mass loss rate reaches very high values, log| Ṁ1| > −2, and R1/RL significantly exceeds 1).
while WD+MS binaries formed by stable mass transfer are expected to be underrepresented in the sample, most of the objects in the sample are inconsistent with formation through stable mass transfer and are likely products of CEE.

CONCLUSIONS
In this work, we presented follow-up spectroscopic observations of 31 MS + WD binaries with orbital periods ranging from ∼ 40 − 300 d and AstroSpectroSB1 solutions from Gaia DR3.These were selected from the broader sample of candidates selected by Shahaf et al. (2024).By jointly fitting the RVs from our spectra and astrometry, we validated the Gaia orbits and obtained tightened constraints on their stellar and orbital parameters.Our main findings are as follows: 1. Presence of a WD companion: Through our joint fitting, we directly constrained the flux ratio S = F 2 /F 1 of our objects.Given these values, we conclude that 26 objects must host dark secondaries to within 1-σ errors, in support of the majority indeed hosting WDs.For the remaining objects, while we cannot rule out the possibility of triples given the errors, they are also consistent with dark companions (Section 5.1).
2. Reliability of Gaia astrometric solutions: Comparing the astrometry-only and astrometry + RV solutions, we find that only 2 out of the 31 systems have significant discrepancies in the orbital parameters (Section 5.2), implying a high success rate for AstroSpectroSB1 solutions.
3. Comparison to literature PCEBs: These objects have P orb = 40−300 d.For a given WD mass, this is shorter than predicted for formation through stable MT but longer than the majority of pre-viously known PCEBs (Figure 7), making their mass transfer histories uncertain.
4. Eccentricity: We find eccentricities that range from 0.02 to 0.45, with a median at ∼ 0.07.This is more circular than MS + MS binaries at comparable periods but more eccentric than most MSP + WD binaries and significantly more eccentric than predicted by the theoretical eccentricity -period relation for stable MT products (Phinney 1992; Figure 8).
5. Formation through CEE : We ran single star MESA models to obtain realistic models of the donor star and trace the binding energy of its envelope.We used the α-formalism to calculate final orbital separations for mass transfer beginning at different points on the donor's evolution.We find that if mass transfer begins when the donor is on the TP-AGB, the wide orbits of our systems can form without the need to invoke recombination energy for sufficiently large values of α CE .However, by including recombination energy, such orbits can form over a wider range of parameters.We also find that if we only consider the outer region of the envelope (above the final orbital separation of the PCEB) to be ejected, it is significantly easier to form the wide orbits (Section 6.1).

Formation through Stable MT :
We also ran several MESA models to test whether MT can remain stable given the mass ratios of the companion and WD progenitor.Both single star models of the donor's response to mass loss as well as binary models of mass transfer from an RGB and AGB donor suggest that mass transfer will become unstable (Section 6.2).
Figures 17 and 18 shows the observed photometry and best-fit SEDs for all 31 objects.Where available, GALEX NUV/FUV points are plotted, but they were not included in the SED fitting to avoid potential contamination from hot WD.Table 5 summarizes the resulting best-fit parameters.
Figures 19 and 20 show RV curves for all 31 objects using parameters from 50 random draws of the posterior resulting from the fitting of just the Gaia solution.We also overplot our follow-up RVs -we see that these are generally in agreement with the predicted curves, suggesting the Gaia solutions are mostly reliable, though there are a few exceptions.The RV residuals from the best-fit model are plotted on the lower panels.We see that for most objects, the residuals are within ∼ 2.5 km s −1 .Figures 21 and 22 show RV curves resulting from the joint fitting of Gaia solutions and our follow-up RVs.The spread in the RV curves are noticeably smaller in this case, illustrating the tightened constraints on the parameters.We also see that the RV residuals are smaller, typically being less than ∼ 0.2 km s −1 .On Table 6, we summarize the best-fit parameters from the joint fitting for all objects.   2 0 1 6 .80 2 0 1 6 .8 5 2 0 1 6 .90 2 0 1 6 .9 5

Figure 1 .
Figure 1.Left: Distance and eccentricity against orbital period of the full sample of MS + WD candidates from Shahaf et al.(2024).The 31 objects presented in this paper are marked with star markers.We have also highlighted ∼ 40 more objects for which RV follow-up is ongoing with diamond markers.Right: The same objects plotted on top of a Gaia (extinction-corrected) CMD, in blue.

Figure 2 .
Figure 2. Results from SED fitting of two example objects, J0145-0208 (top) and J0547-2018 (bottom).Left: Corner plot showing the four fitted parameters.Center : Isochrone predictions of the radius versus temperature from 50 random draws of the posterior.The marker indicates the best-fit current values from the SED fitting.Right: Best-fit SED model plotted compared to the observed photometric points from WISE, 2MASS, Gaia, and GALEX along with the residuals.Note that GALEX points are plotted here, but were not used in the SED fitting.J0145-0208 does not show significant FUV excess, but J0547-2018 does.

Figure 4 .
Figure 4. Secondary mass M2 against the flux ratio S for several example objects.The black line is the expected relation obtained from the AMRF (Equation1) calculated using the best-fit orbital parameters derived in Section 4; black point marks the best-fit values of M2 and S. The blue solid line and orange dashed lines correspond to theoretical curves for single MS star and equal-mass inner binary companions.On the top row are "group 1" objects (as described in Section 5.1) which must host dark companions: the constraints on M2 and S are inconsistent with an inner MS+MS binary for any S. Bottom row shows "group 2" objects, for which the flux ratio constraint cannot rule out an inner MS+MS binary.

Figure 5 .Figure 6 .
Figure5.Left: Histogram showing the distribution of the best-fit flux ratio S for the 31 objects in our sample.The majority of objects have S ∼ 0, which implies the less luminous secondary contributes little light to the system.Right: Eccentricity e against the flux ratio S of our sample resulting from the joint Gaia + RVs fits.The markers and colors of the points represents the groups that each object belongs to, as described in Section 5.1 and shown in Figure4.All of the systems are consistent with a dark WD secondary.

Figure 7 .
Figure 7. Orbital period, P orb , against WD mass, M2.The colors of the points represent the masses of the luminous companion, M1.On the right axis, we show the minimum orbital separation, aperi, corresponding to the period assuming M1 = 1 M⊙ and M2 = 0.6 M⊙.We plot the "short-period PCEBs" discovered with SDSS by Rebassa-Mansergas et al. (2007) as well as those from Hernandez et al. (2021, 2022a,b) with circle markers."Wide PCEBs" which include IK Peg Wonnacott et al. (1993) and the systems from Yamaguchi et al. (2024b) are plotted with diamond markers.We also plot the self-lensing systems discovered with Kepler (Kruse & Agol 2014; Kawahara et al. 2018; Yamaguchi et al. 2024a) with triangle markers.The solid line marks the maximum mass that is reached by the giant progenitor at each WD mass.The dashed line is the track along which binaries undergoing stable MT are expected to evolve Rappaport et al. (1995).

Figure 9 .
Figure 9. Left: HR diagram of a MESA model of a 1.0 M⊙ star.Middle: Final separation, a f , against the initial separation at the onset of mass transfer, ai.The dashed line represents the case where all of the orbital energy loss goes into envelope ejection (αCE = 1), but where no recombination energy is included (αrec = 0).The three solid lines are cases for αCE = 0.3, 0.6 and 0.9, with the inclusion of 50% of the available recombination energy (αrec = 0.5).The horizontal dashed lines at 0.01, 0.15, and 0.7 AU represent the minimum separation below which the luminous primary would not fit in the orbit, the minimum separation of the "wide PCEBs", and the typical separation of the objects from this work.Black arrows indicate initial separations where the addition of recombination energy has resulted in the overall binding energy being positive, where the envelope is formally unbound.The upper and lower rows correspond to cases where the binding energy is calculated over the entire envelope vs. just above the final orbital separation of 0.7 AU, respectively.Right: Zoom-in of the region above 0.15 AU of the middle plot in the same row, isolating the regions with thicker line weight which can be reached by the donor in practice (i.e.where the radius is larger than at any previous point in the star's evolution).Wide orbits can form without the addition of any recombination energy if αCE is sufficiently large.This is true with both definitions of the envelope binding energy, though the range of initial separations over which they form significantly increased when just the outer envelope is considered.

Figure 10 .
Figure10.Same as Figure9, but for a 1.5 M⊙ model.For this model, we see that wide orbits like those of our systems (a f > 0.7 AU) do not form without recombination energy (orange dashed line) if the binding energy of the entire envelope is considered.Meanwhile, there is a small range of initial separations over which they can form if only the outer envelope is considered.

Figure 11 .
Figure 11.Same plots as Figure 9 but for a 3 M⊙ model.
Figure12.Envelope mass against separation at initial Roche lobe overflow for the 1, 1.5, and 3 M⊙ models.As in Figure9, regions where ai is larger than at any previous point have been plotted with a thicker line weight.The black arrow indicates the direction of evolution.We see that the range of initial separations needed to form the orbits of our systems in each of these cases (seeFigures 9, 10, and 11)   occur primarily during the thermal pulses.

Figure 14 .
Figure14.Period distribution of solar-type MS + MS binaries (e.g.Raghavan et al. 2010).In blue is the region corresponding to systems where the first interaction will occur when one star becomes a red giant.If such binaries undergo CEE, these are expected to form close PCEBs.In red are systems where the first interaction will occur when one star is on the AGB, which is expected to form wide PCEBs.These systems each make up ∼ 10% of the total population.
Figure 1 of Shahaf et al. (2024) (dotted line) and as shown in Figure 11 of El-Badry (2024), results in the condition M WD /M ⋆ > 0.6.A stable MT system with a WD mass of ∼ 0.4 M ⊙ would thus require M 1 ≲ 0.7 M ⊙ , and very few of our 31 objects have primary masses this low.

Figure 19 .
Figure19.RV curves resulting from the fitting of the Gaia solutions only for all 31 objects (continued to next figure).We plot curves corresponding to parameters from 50 random draws of the posterior to illustrate the uncertainity.We overplot our follow-up RVs, and plot their residuals from the best-fit model on the lower panels of each plot.

Figure 21 .Figure 22 .
Figure21.Same as Figure19, but for the joint fitting of both the Gaia solutions and follow-up RVs.We see that there is much less spread in the predicted RV curves, corresponding to improved constraints on the orbits.
Yamaguchi et al. (2024b)ram of literature MSP + WD binaries from the ATNF pulsar catalog(Manchester et al. 2005; compiled by Hui et al. 2018).The triangle and circle markers differentiate systems with WD masses above and below 0.45 M⊙ respectively, which marks the approximate boundary between He and CO/ONeMg WDs.We also plot the objects from this work with star markers and the wide PCEBs fromYamaguchi et al. (2024b)with diamond markers.The solid line is a relation expected for He WD systems formed from stable MT from Eggleton 1983 formula), and a f is the final separation at envelope ejection (i.e. the separation of the resul-

Table 2 .
All RVs collected for the 31 objects of this work.

Table 5 .
Best-fit parameters from SED fitting of the MS star.From left to right, the columns are: object name, extinction (with typical uncertainties of ∼ 0.02 mag), star mass, Equivalent Evolutionary Phase, parallax, initial metrallicity of the star, present-day metallicity, effective temperature, and radius.

Table 6 .
Best-fit parameters from the joint fitting of the astrometry and RVs.From left to right, the columns are: object name, period, eccentricity, inclination, position angle of the ascending node, argument of periapsis, periastron time, center-of-mass velocity, WD mass, MS star mass, flux ratio.