Planets surviving common envelope evolution should be rare, massive, and orbit very close to their host star

Introduction

The fate of exoplanets is tied to the evolution of their environments, including the life cycle of their host stars. It has been established that exoplanets are present at different stages of the evolution of the host star, including its death. Observational evidence for planetary survival in the late stages of stellar evolution includes planets around pulsars (e.g. A. Wolszczan & D. A. Frail [65]), accreting planetary material onto polluted white dwarfs (e.g. B. Zuckerman et al. [72]; D. Koester et al. [29]; M. Hollands et al. [22]), and some intact planet candidates around white dwarfs (e.g. K. L. Luhman et al. [34]; T. Gänsicke et al. [18]; K. Zhang et al. [69]; J. W. Blackman et al. [2]). Three of these intact candidates are in close orbits around their star (WD 1856+534 b, A. Vanderburg et al. [61]; WD 0310-688 b, M. A. Limbach et al. [32]; and WD J0914+1914 B. T. Gänsicke et al. [18]). WD 1856b was recently confirmed to be a planet with a mass between $1 - 6\,M_{\mathrm{Jup}}$ (M. A. Limbach et al. [33]) and an orbital period of 1.4 days (A. Vanderburg et al. [61]).

The main challenge of explaining the presence of planets in close orbits around stellar remnants is that the planet has to survive the final stages of the star’s post-main-sequence evolution, specifically the Red Giant branch (RGB) and Asymptotic Giant Branch (AGB). Planets with larger orbital separations should be able to avoid engulfment and are expected to survive the giant stages of their host star. In contrast, the present-day orbits of close-in planets around white dwarfs, like WD 1856b, lie well within the maximum radius reached by the host star during its giant phases. This implies that either the planets survived engulfment or migrated inward after the white dwarf formed.

That challenge has made WD 1856b the subject of considerable speculation regarding its formation history (C. E. O’Connor et al. [45]; A. P. Stephan et al. [59]; D. J. Muñoz & C. Petrovich [38]; F. Lagos et al. [30]; L. Chamandy et al. [6]; A. Merlov et al. [37]). Two

[email protected] primary formation pathways have been proposed: common envelope (CE) evolution and high-eccentricity (HE) migration. HE migration involves an initial excitation of the planet’s orbit to a high eccentricity, potentially via Kozai-Lidov oscillations, planet-planet scattering, or stellar flybys, followed by tidal circularization to a close orbit. Although high eccentricity migration post-RGB or -AGB phase could, in principle, produce the observed configuration of systems hosting intact giant planets such as WD 1856b (as well as perturb smaller bodies, e.g. (Y. Li et al. 2025)), there is still not enough observational evidence to tell whether this mechanism is an active formation channel [45].

In the CE scenario, the planet initially orbits a progenitor star that undergoes expansion during its red giant phase, leading to engulfment and subsequent orbital decay due to dynamical friction within the envelope. Once orbiting within the giant star’s envelope, the planet must deposit enough energy into the gas to eject the star’s outer layers and end migration before the planet is destroyed. Since early works and observations on CE evolution [47, 4 and others that followed], this mechanism has become widely accepted as the explanation for the presence of close-in stellar companions to white dwarfs. However, it is less clear whether the process can extend to planetary companions. Because planets have much lower masses than stars, they release less orbital energy during inspiral, making it harder to unbind and eject the envelope. However, previous studies have concluded that CE evolution for planetary-mass companions is borderline energetically feasible [30, 67], especially if additional energy sources, like recombination energy and orbital energy from other planetary companions can contribute [6]. This could make it possible for CE to form close-in planets around white dwarfs.

While a number of studies have investigated the orbital evolution and consequent survivability of planets during CE evolution, these works have focused on individual systems or on the mechanisms of survival. In this paper, we focus on predicting the observable population of planets that could emerge if CE evolution is responsible for producing close-in planets around white dwarfs. Other studies [42, 35] have performed similar analyses, but they mostly focused on planets that survived beyond the outer reaches of the envelope. Here, we focus on the planets that get engulfed, see which of them can survive, and predict the properties of the post CE population.

In this work, we estimate the number of planets that could result from CE evolution and predict the post-CE planet population distributions of semi-major axis and planet mass. We simulate a population of giant-planet hosts drawn from observed exoplanet mass and orbital semi-major axis distributions. We vary certain simulation parameters, like the tidal $Q$ value in our tidal prescription, the efficiency of envelope ejection $\alpha$ in the CE prescription, and the host star mass. We find that planets surviving CE evolution should be rare, massive, and orbit near their tidal disruption radii (the Roche limit). Although WD 1856b is the only confirmed example of a closely-orbiting giant planet around a white dwarf, its rarity may be reflect the incompleteness of WD surveys rather than true scarcity. Therefore, the characteristics of our resulting CE population stand in contrast to those of WD 1856b, suggesting it may have formed by a different method.

This paper is organized as follows. In § 2, we describe our methods for simulating the population of giant planet hosts. We also discuss the tidal evolution and CE prescription that we use to evolve the systems, as well as the termination conditions for our simulations. § 3 details the trends in the surviving populations for each simulation case study. In § 4, we discuss the occurrence rates of the surviving planets and compare them to the known and candidate white dwarf planets to determine how well their observed properties match the predictions of the CE scenario. § 5 summarizes our main findings.

Methods

To predict the population-level characteristics of planets that migrated through CE evolution, we ran simulations that evolve single-planet systems along with their host stars. Our approach consists of generating a population of planets drawn from observationally motivated distributions and coupling their orbital evolution to the host star’s evolution from the main-sequence to the white-dwarf stage. We conduct several simulation “runs,” each with a different combination of assumptions to test the robustness of our predictions. We visually summarize the steps of our full evolutionary sequence in Figure 1 and provide an in-depth description of each element of our model in the following subsections.

Initial System Parameters

For each run, we randomly generated 10 million planets and evolved the systems with different tidal prescriptions, stellar masses, or CE efficiency parameters. The planet’s initial parameters were randomly drawn from known giant planet distributions. The mass $M_{\mathrm{p}}$ and semi-major axis $a$ were sampled from the California Legacy Survey (CLS) distributions [14]. We sampled the initial masses from 0.1 to 13 $M_{\mathrm{Jup}}$ and semi-major axis from 0.03 to 30 AU, but chose to optimize our computation time by only running our full evolutionary sequence for planets likely to survive CE. In particular, we only ran simulations for planets with mass $> 2 M_{\mathrm{Jup}}$ and $a$ between 1 and 15 AU. Preliminary tests confirmed that there were no CE survivors outside of these bounds for most of our scenarios. These initial tests also showed that planets below 2 $M_{\mathrm{Jup}}$ and 1 AU never survived; therefore, we immediately mark these planets as destroyed. For planets with $a > 15$ AU, initial tests showed that they never go through CE; therefore, we mark such planets as non-CE survivors. We note that for the cases with $\alpha= 20$, we relaxed the inner cutoff on $a$ and evolved planets with $a < 15$ AU, because initial tests showed survivors down to smaller separations. The full initial planet mass and semi-major axis distribution are shown in Figure 2. The eccentricity values $e$ were drawn from the distribution found by D. M. Kipping [28]. We calculated the planet radius $R_{\mathrm{p}}$ from the J. Chen & D. Kipping [7] mass-radius relationship. Each planet was assigned to a host star with a main-sequence mass chosen to be either 1 or 2 $M_{\odot}$. We chose 2 $M_{\odot}$ for our realistic and illustrative cases because stars of this mass are the most common progenitors of white dwarfs.

MESA Models

We computed stellar evolution models for initial masses of 1.0 and 2.0 $M_{\odot}$ at Solar metallicity ($Z = 0.02$) using Modules for Experiments in Stellar Astrophysics (MESA, version r24.08.1, [48, 49, 50, 51, 52, 27]). Each model was evolved from the pre-main-sequence through the thermally pulsating asymptotic giant branch (TP-AGB) and subsequent white dwarf cooling. We adopted the GS98 solar composition [17] with Type 2 OPAL opacities. The energy equation was solved with the dedt option during core hydrogen burning, switching to the eps_grav formulation for all subsequent phases to improve numerical stability through the helium core flash and AGB thermal pulses. The gold2 solver tolerances were enabled throughout to maintain tight energy conservation.

Mass loss was treated with the D. Reimers [56] prescription on the red giant branch and the T. Blöcker (1995) prescription on the asymptotic giant branch, with the transition governed by a central hydrogen mass fraction threshold of $10^{-4}$. Angular momentum transport was included via rotationally induced instabilities, with the Tayler–Spruit dynamo, Eddington–Sweet circulation, Goldreich–Schubert–Fricke instability, and secular shear all enabled. During the AGB phase, we applied exponential convective overshoot [21] at all convective boundaries with $f_{\mathrm{ov}} = 0.014$ and $f_{0} = 0.004$, and employed additional mesh refinement in and around overshoot regions to resolve thermal pulses. The AGB phase was terminated when the envelope mass fell below $10^{-2} M_\odot$.

Orbital Evolution

Once we initiated the population, we evolved the planet’s semi-major axis $a$ and eccentricity $e$ with a 4th-order Runge-Kutta integration method. We over-sampled the MESA time array by a factor of 4 to obtain shorter timesteps and improve the accuracy of the time-evolution calculations. The change in semi-major axis, $da/dt$, had two components: one due to tides and one due to stellar mass loss. The eccentricity only changed due to tides:

The change in $a$ due to stellar wind mass loss is given by:

where $M_{\star}$ is the stellar mass at the corresponding timestep and $\dot{M}_{\star}$ is the stellar mass loss rate computed directly by MESA rather than calculated manually, which minimizes numerical errors. We also expected that tides raised on both the planet’s and star’s surface would play an important role in the planet’s orbital evolution. To model the tidal evolution of orbits, we used the equilibrium tide formalism of P. Goldreich & S. Soter [16] to second order in eccentricity, written as in J. Nordhaus et al. (2010):

where $R_p$ and $R_{\star}$ are the radii of the planet and the star at a given time, $Q'_p$ and $Q'_{\star}$ are the modified tidal dissipation factors of the planet and star, respectively. This low-eccentricity approximation describes well the evolution expected for the majority of our sample. For the small portion of our sample that are initialized with larger eccentricities, we may be slightly underestimating tidal dissipation, but the computed orbital evolution will be qualitatively accurate across all tested values of eccentricity. The constants $K_1$ and $K_2$ are given by:

We ran simulations with both a constant and a variable star and planet tidal quality factor $Q'_{\star,p}$. These modified tidal factors are given by $3Q/2k_2$, where $Q$ is the specific tidal dissipation function, and $k_2$ is the Love number. The stellar static tidal parameter $Q'_{\star}$ was taken to be $10^5$ or $10^6$, which are within the range of accepted values in the literature and are based on observational constraints ([46, 25]). The planet static tidal parameter $Q'_p$ was taken to be $10^{5.5}$. The variable $Q'$s were calculated at each timestep, where the subscripts $\star$, p stand for star or planet following D. Fabrycky & S. Tremaine [12]; B. M. Hansen [19]:

where the apsidal motion constant $k_{\star,\mathrm{p}} = 0.5\frac{Q_{\mathrm{E}_{\star,\mathrm{p}}}}{1-Q_{\mathrm{E}_{\star,\mathrm{p}}}}$ is related to the P. P. Eggleton et al. [11] coefficient $Q_{\mathrm{E}}$, for which we used 0.028 for stars and 0.30 for gas giants. The viscous timescale $t_{V_{\star,\mathrm{p}}}$ was taken to be 1.5 years for stars and 0.01 years for planets (S. Naoz et al. [39]; C. Petrovich [53]). Figure 3 shows the time evolution of variable $Q'_{\star}$ and $Q'_{\mathrm{p}}$ in one of our simulations as the planet and star evolve.

For the stages of evolution where the host star has a significant convective envelope, the viscous time scale had to be modified when the eddy turnover timescale (F. A. Rasio et al. [55]) became greater than half the orbital period of the planet (i.e., when $\tau_{\mathrm{d}} > P/2$) (J. P. Zahn et al. [68]; E. Villaver & M. Livio [62]). We calculated this effective viscous timescale $t_{V\mathrm{eff}} = t_V/f$, where:

$R_{\mathrm{env}}$ is the radius at the bottom of the convective envelope, and we defined the core boundary within the radial profile of the star where it has a hydrogen mass fraction less than 0.1, as was assumed in J. D. M. Dewi & T. M. Tauris [9]. The term $L_{\star}$ is the luminosity of the star. We used $f = 1$ otherwise.

Common Envelope Evolution

The planet can be engulfed by its host star if it is orbiting close enough during stages where the star increases in radius, such as the RGB or AGB. This event initiates a phase called Common Envelope evolution. During this process, the primary star fills its Roche lobe and wraps the companion planet into a common gas envelope. As the planet moves through this envelope, it experiences drag forces and begins to lose angular momentum and orbital energy. This usually causes a decrease in the planet’s orbital separation. In the end, either enough orbital energy is injected into the envelope to eject it, and the planet remains in orbit, or the orbit continues to shrink until the planet is at or below its host star’s convective zone. Then there can be two outcomes, either the planet is dissolved into the convective zone or it merges into the star’s core (R. Yarza et al. [67]). In our simulations, we do not distinguish between these two destruction scenarios, we designate them both as “accreted”.

The system will undergo common envelope evolution if $R_{\star} > R_{\mathrm{Roche}}$. This means that mass transfer begins when the star radius $R_{\star}$ reaches or exceeds its Roche radius $R_{\mathrm{Roche}}$ defined by the P. P. Eggleton [10] approximation:

where $q = M_{\star}/M_{\mathrm{p}}$. After CE, the stripped core with mass $M_{\mathrm{c}} = M_{\mathrm{WD}}$ will be left. We define a new variable, $a_{\mathrm{engulf}}$, to describe the semi-major axis of the planet at the instant when Roche lobe overflow begins. For the CEE approach, we used the energy formalism by R. F. Webbink [64] combined with the numerical parameter $\alpha\lambda$ studied by M. de Kool et al. [8], where $\alpha$ is related to the CE efficiency and $\lambda$ to the central concentration of the body:

where $M_{\mathrm{c}}$ is the core mass of the host star, $M_{\mathrm{env}}$ is the mass of the stellar envelope, $M_{\star}$ is the total mass of the star, and $M_{\mathrm{p}}$ is the planet mass. We note that $\alpha$ is more specifically defined as the ratio of the energy needed to unbind the envelope to the orbital energy released during inspiral. In the simulations, this efficiency parameter was fixed at 0.41, 1, 2, 2.5, or 20 (we motivate such high values of $\alpha$ below).

We calculated $\lambda$ at each timestep by integrating the binding energy of the envelope:

where $r_{\mathrm{L}} = R_{\mathrm{Roche}}/a_{\mathrm{engulf}}$ is the dimensionless Roche lobe radius of the star, so we can use $a_{\mathrm{engulf}}r_{\mathrm{L}} = R_{\mathrm{Roche}}$. We defined $M_{\mathrm{env}}$ as the total mass above the star’s core-envelope boundary, which is the location where the hydrogen mass fraction first exceeds 0.1 as defined in § 2.3. Here we calculate the binding energy $E_{\mathrm{env}}$ in a similar manner to Equation 5 of J. D. M. Dewi & T. M. Tauris [9], which included contributions of gravitational and internal thermal energy as:

where $k$ is Boltzmann’s constant, $\mu$ is the mean molecular weight $\simeq0.615$, $m_{\mathrm{H}}$ is the atomic mass of hydrogen, and $\rho(r)$ and $T(r)$ are the density and temperature radial profiles of the star, respectively.

Note that we did not include ionization (recombination) energy in our calculations. This contribution can be comparable in magnitude to the thermal energy when the host star reaches larger radii and, if included, would drive the total envelope energy toward zero or even positive values if the ionization energy becomes sufficiently large:

where $E_{\mathrm{grav}}$ is the gravitational potential energy with magnitude $V$ and $E_{\mathrm{th}}$ is the thermal energy defined as with $E_{\mathrm{th}} = -E_{\mathrm{grav}}/2 = V/2$ from the virial theorem. At large stellar radii, the ionization energy can approach or even exceed the thermal energy ($kT \sim$ binding energy per particle), which we express as $E_{\mathrm{ion}} \sim E_{\mathrm{th}}$. This would imply that no additional energy is required for envelope ejection. In this limit, the envelope becomes effectively unbound, and the parameter $\lambda$ would diverge to very large values (e.g. $\lambda\gtrsim100$). Given the uncertainty in how efficiently ionization energy contributes to unbinding the envelope, we decided to exclude it from our calculations of the envelope binding energy. Previous work has also deemed it unnecessary to include recombination energy in population synthesis studies ([70]). However, to illustrate the potential impact of this additional energy, we included models assuming an envelope ejection efficiency parameter $\alpha$ greater than 1 (up to 20) in our simulations. Figure 4 shows the time evolution of $\lambda$ for a 1 and 2$M_{\odot}$ host star.

The final semi-major axis $a_{\mathrm{f}}$ can then be solved for by substituting Equation (9) into Equation (8):

However, because we already pre-calculated $\lambda$ at Roche overflow in our simulations, we substituted $\lambda$ explicitly into Equation (8) after solving for $a_{\mathrm{f}}$ to track the orbit evolution:

Termination conditions

Not all planets survive during the CE phase of the simulations. We therefore defined planetary “termination” as the point at which a planet can no longer maintain a stable orbit or a temperature suitable for survival. Once the host star has filled its Roche lobe ($R_{\star} > R_{\mathrm{Roche}}$), we calculated $a_{\mathrm{f}}$ for CE evolution. Then we assessed the potential destruction of the planet under two conditions. First, we considered whether a planet is tidally disrupted inside the star, and is either dissolved in the convective zone or completely accreted onto the stellar core. We check if this takes place by testing whether $a_{\mathrm{f}}$ is less than the planet’s tidal disruption limit $a_{\mathrm{disrupt}}$. We note that the tidal disruption radius for the planet is also known as the Roche limit, but to avoid confusion with the Roche radius of the star discussed in § 2.4, we will refer to this distance as the tidal disruption radius. For Jupiter-like planets, we can solve for $a = a_{\mathrm{disrupt}}$ from Equation (7) after setting $R_{\mathrm{Roche}} = R_{\mathrm{p}}$ and $q' = M_{\mathrm{p}}/M_{\mathrm{WD}}$:

The second condition is when the planet ends up with a very short orbit, enough to be under the evaporation radius of the star $a_f < a_{\mathrm{evap}}$ and become thermally destroyed. This is where the planet can experience significant irradiation from the stellar core during or after the envelope has been ejected. This distance was computed as in N. Soker [58]:

where $M_p$ is the planet mass and $M_{\mathrm{Jup}}$ is Jupiter mass. Together, these criteria provide physical constraints for planetary survival after CE evolution.

Results

Possible Outcomes

Our simulations yield three distinct evolutionary outcomes for the planets in our modeled systems. The first possibility is that the planet never enters CE evolution and its orbit remains wide and only evolving in response to stellar mass loss. This outcome is typically associated with planets that begin at large separations, well beyond the maximum radius of the host star during its AGB phase. In such cases, the planet avoids engulfment altogether, remaining gravitationally bound to the remnant white dwarf. The second scenario is where the planet undergoes CE evolution and is destroyed because the orbital energy released during inspiral is insufficient to eject the envelope material and avoid either tidal disruption or evaporation. The third outcome is that the planet undergoes CE evolution and survives. We find that these survivor planets are generally massive and usually located slightly beyond their tidal disruption radii, where tidal disruption by the star and the planet’s self-gravity are balanced. Although this is a rare outcome in our simulations, it is feasible for planets to survive CE evolution under specific assumptions (especially assuming high values of envelope ejection efficiency) and finely-tuned initial orbital conditions. Figure 5 illustrates the orbital evolution of systems representative of the three outcomes using an illustrative set of assumptions, which is detailed in Section 3.3.

Realistic Case: $2M_{\odot}$, variable $Q'_{*,p}$, $\alpha= 0.41$

We first performed a batch of simulations for the selected subset of the $10^7$ initial systems (as explained in § 2.1) with what we consider to be the most realistic set of assumptions we could choose. We chose a stellar mass of $2M_{\odot}$ because this closely matches the most common progenitor for white dwarf stars in our galaxy. We used a tidal prescription that allows the tidal quality factor to change with the evolution of the host star,

and we used a common envelope efficiency value that is in the upper range of plausible values of $\alpha$ by observations of the population of post-common-envelope white dwarf/brown dwarf binary stars ([70]). Out of $10^7$ initial systems in this instance of parameters, there were no systems where the planets survived CE evolution. This implies that if our assumptions indeed hold, planets surviving common envelope evolution are likely very rare, with an occurrence rate $\lesssim10^{-7}$.

Illustrative Case: $2M_{\odot}$, static $Q'_{*,p}$, $\alpha= 2$

Because our most realistic set of simulation inputs and assumptions did not result in any surviving planetary systems, we also performed simulations with slight permutations to those realistic assumptions in order to illustrate the properties of the common envelope planet population. For this illustrative case, we continued to assume the planet orbited a $2M_{\odot}$ star, and used a slightly simpler static tidal quality factor prescription with $Q'_{*}=10^5$ and $Q'_{p}=10^{5.5}$. Most importantly, we assumed a higher envelope ejection efficiency parameter $\alpha$ of 2.0 which implies some modest amount of additional energy has been added into the system, perhaps from recombination reactions. These parameters do result in a diversity of post-CE outcomes; a portion of the initial planetary distribution of our illustrative-case simulations is shown in Figure 7, with color coding corresponding to the planets’ post-CE fates. After the cuts in parameter space made for computational efficiency (as described in § 2.1), we were left with 1,461,550 systems to run through our model. About 0.074% of our simulated systems result in CE survivors, 18.27% of planets are destroyed via evaporation, 1.7% are destroyed via tidal disruption, and $\sim80\%$ of systems never go through CE.

Our simulations resulted in 1084 surviving planets, which is a small fraction ($\sim10^{-4}$) of the total 10 million initial systems. Even in this case with the additional energy to help the planet unbind the stellar envelope, CE survival is still a rare outcome. Importantly, because we have a sample of surviving planets, we can now examine the population characteristics of post-CE planets. We find that the surviving planets are usually massive ($\gtrsim10M_{\mathrm{Jup}}$) and orbiting very close to their star ($\lesssim0.006$ AU, $\sim1.3R_{\odot}$). The final distribution of the planets surviving CE in our illustrative case for a $2M_{\odot}$ is shown in Figure 6.

We further investigated the distribution of the final semi-major axes for surviving planets by comparing them to the orbital radii at which they would either be tidally disrupted ($a_{\mathrm{disrupt}}$, also called the tidal disruption radius or the Roche limit) or the location within the star that they would have likely been evaporated ($a_{\mathrm{evap}}$, which we call the evaporation radius). We show the ratio of the planets’ final semi-major axes to these limits in Figure 8. We find that surviving planets pile up around as close an orbital separation as they can survive, either around the evaporation radius or the tidal disruption limit. All of the planets in this illustrative case had semi-major axes $a_f < 1.4\,a_{\mathrm{disrupt}}$ and $a_f < 1.6\,a_{\mathrm{evap}}$.

Varying Parameters

In addition to running our simulations with these two specific sets of assumptions, we also systematically explored how the variations in key physical parameters (host star mass, tidal prescription, and common envelope efficiency) influence the final population of planetary systems that emerge from CE evolution. We ran a grid of simulations, each evolving the subset out of $10^7$ initial systems, over all possible combinations of two stellar mass values, three tidal prescriptions/quality factors, and five values of common envelope efficiency. We summarize the number of surviving planets and their typical orbital parameters across our parameter space in Table 1.

Parameters	1M$_\odot$					2M$_\odot$
$Q'_{\star}$	$Q'_{\mathrm{p}}$	$\alpha$	$N_{\mathrm{surv}}$	$a_{\mathrm{f\,med}}$ (AU)	$M_{\mathrm{p\,med}}$ (M$_{\mathrm{Jup}}$)	Occ. Rate	$N_{\mathrm{tr}}$	$N_{\mathrm{surv}}$	$a_{\mathrm{f\,med}}$ (AU)	$M_{\mathrm{p\,med}}$ (M$_{\mathrm{Jup}}$)	Occ. Rate	$N_{\mathrm{tr}}$
variable	variable	0.41	0	-	-	-	0	-	-	-	-
		1	0	-	-	-	142	0.0054	11.17	$9.9 \times10^{-6}$	10
		2	0	-	-	-	578	0.0075	10.93	$4 \times10^{-5}$	53
		2.5	737	0.0042	12.63	$2.5 \times10^{-5}$	85	1046	0.0062	11.17	$7.3 \times10^{-5}$	101
		20	229907	0.0098	9.54	$7.7 \times10^{-3}$	14949	113407	0.0078	10.27	$7.9 \times10^{-3}$	8734
$10^{5}$	$10^{5.5}$	0.41	0	-	-	-	0	-	-	-	-
		1	0	-	-	-	0	-	-	-	-
		2	0	-	-	-	1084	0.0049	12.03	$7.6 \times10^{-5}$	147
		2.5	458	0.0042	12.62	$1.6 \times10^{-5}$	79	3845	0.0048	11.98	$2.7 \times10^{-4}$	465
		20	147215	0.011	9.15	$4.9 \times10^{-3}$	8473	155628	0.013	9.10	$1.1 \times10^{-2}$	8991
$10^{6}$	$10^{5.5}$	0.41	0	-	-	-	0	-	-	-	-
		1	0	-	-	-	0	-	-	-	-
		2	0	-	-	-	989	0.0048	12.15	$6.9 \times10^{-5}$	126
		2.5	598	0.0041	12.67	$2 \times10^{-5}$	86	4072	0.0048	12.11	$2.9 \times10^{-4}$	533
		20	135472	0.012	8.91	$4.6 \times10^{-3}$	7297	197933	0.013	8.78	$1.3 \times10^{-2}$	10343

Progenitor Mass

We first investigated how the mass of the host star influences the resulting planet population. In addition to performing simulations with a $2M_{\odot}$ progenitor (like in our illustrative and realistic cases), we also conducted simulations for a lower host mass of $1M_{\odot}$ to explore how a less massive envelope and different evolutionary timescale affect the outcome. The comparison between the resulting planetary populations produced in these two cases is shown in Figure 9 and Table 1.

We find that overall in our simulations, there are roughly the same number of CE survivors around a host star of $1M_{\odot}$ or $2M_{\odot}$. In a case by case basis (that is, across different permutations of common envelope efficiency and tidal prescriptions), it is often but not always true that there are fewer CE survivors for a $1M_{\odot}$ host than a $2M_{\odot}$ host. This result is slightly counterintuitive because a $1M_{\odot}$ star has a lower-mass envelope, which would (all else being equal) be easier to eject, allowing more planets to survive. However, even though the mass of the envelope is lower for a $1M_{\odot}$ star, it turns out that $\lambda$ is lower, implying the envelope is more centrally concentrated and more strongly bound to the star, making it harder to eject. The mean $\lambda$ at onset for CE survivors around a $1M_{\odot}$ host is 0.18, while for the survivors around a $2M_{\odot}$ this is 0.97. The larger values of $\lambda$ mean that less energy from the planet is required to expel the envelope, therefore allowing a larger number of CE survivors.

Overall, regardless of the mass of the host star, we found that the population of post-CE planets are qualitatively similar to our illustrative case. In both host-star scenarios, only the most massive planets survive, and they end up on very small final semi-major axes.

In scenarios where a planet does not undergo CE, stellar evolution will lead to a $0.57M_{\odot}$ WD for a $2M_{\odot}$ progenitor and $0.54M_{\odot}$ WD for a $1M_{\odot}$ progenitor. However, when a planet and star undergo CE, the resulting WD mass spans a range of values set by the evolutionary phase at which CE occurs. Core mass grows continuously as the star ascends the giant branch. Therefore, earlier CE episodes (such as during the RGB) will result in lower-mass white dwarfs while later ones will leave slightly more massive remnants. This effect produces a spread in final masses of WDs hosting CE survivors, rather than the single mass value expected for systems that never undergo CE. Figure 10 shows the final WD mass distributions for the systems that survive CE around a $1M_{\odot}$ and $2M_{\odot}$ progenitor with different values of $\alpha$. For the $1M_{\odot}$ host, white dwarf masses are consistent with CE occurring during the RGB phase. The $\alpha= 20$ case produces a broader distribution that includes lower masses ($\sim0.3\,M_{\odot}$), while the $\alpha= 2.5$ case has WD masses clustered closer to the upper end of the range ($\sim0.45\,M_{\odot}$). For the $2M_{\odot}$ host, there can be two distinct WD mass distributions that can arise. The $\alpha= 20$ case produces survivors from both the RGB

Tidal Prescriptions

($\sim0.4\,M_{\odot}$) and AGB ($\sim0.54\,M_{\odot}$) stages. In contrast, the $\alpha=2.5$ case only produces AGB survivors clustered toward the upper end of the mass range.

Tidal Prescriptions

We next tested how different tidal prescriptions affect the results of our simulations. We compared two tidal prescriptions to assess the role of tidal physics in shaping the post-CE planetary populations. One prescription assumed constant tidal quality factors for both the host star $Q'_{\star}$ and the planet $Q'_{\mathrm{p}}$, while the other allowed them to vary. Variable tidal quality factors may be more realistic because tidal dissipation efficiency is tied to the orbital decay rate and can evolve meaningfully over the course of the stellar lifetime. In the static case, we tested two values for $Q'_{\star}$: $10^{5}$ and $10^{6}$ as explained in § 2.3. The results of these tests are summarized in Table 1, and a comparison of the populations for a particular host star mass ($2M_\odot$) and common envelope efficiency ($\alpha=2$) are shown in Figure 11.

We found that generally, using different tidal prescriptions and different static tidal quality factors does change the number of surviving post-CE planets, but the changes are usually small, and there are no clear rules or patterns to predict whether one prescription will produce more CE survivors than another. An exception to this trend is in the limit of low common envelope efficiency, and cases where the planet has relatively low amounts of orbital energy (planets with relatively low-masses and large $a_{\mathrm{f}}$). In these cases, it appears that a variable $Q'_{\star,\mathrm{p}}$ prescription increases the likelihood that low-energy planets survive CE evolution. This is visible in Figure 11, where the surviving planets in the variable $Q'_{\star,p}$ prescription have masses as low as 7 $M_{\mathrm{Jup}}$, while planets surviving in the static $Q'_{\star,p}$ case all have masses greater than 10 $M_{\mathrm{Jup}}$.

Overall, however, we continue to see the same qualitative picture of the post-CE planet population regardless of the tidal prescription: planets that survive CE are massive, and have small semi-major axes.

Efficiency Parameter

Finally, we examined how variations in the envelope ejection efficiency parameter, $\alpha$, influence the survival of planets after CE. This parameter encapsulates the efficiency with which orbital energy is converted into kinetic energy to eject the envelope. In our simulations, we kept $\alpha$ constant in each run and varied the test value from 0.41 to 20 to capture conservative and optimistic scenarios. We include $\alpha= 0.41$ as the upper limit within the range of 0.2–0.41 estimated from population-wide observations of post-CE WD binaries ([70, 71, 60, 5]). We note that while some studies incorporate some internal energy sources into their calculations of $\lambda$ (specifically M. Zorotovic et al. [71] and J. Camacho et al. [5]), which makes direct comparisons of inferred $\alpha$ challenging, the value we adopt of $\alpha= 0.41$ inferred by M. Zorotovic & M. R. Schreiber [70] assumes no additional energy sources in their calculation of $\lambda$.

The results of these experiments are summarized in Table 1, and a comparison of the population for a particular tidal prescription and stellar mass is shown in Figure 12. We find that modifying the common envelope efficiency parameter $\alpha$ causes by far the largest and most dramatic changes to the post-CE planet population. Increasing $\alpha$ from 0.41 to 20 increases the number of surviving planets by more than 5 orders of magnitude.

While we inferred that post-CE planets are likely rare based on the results of our illustrative and most realistic cases, which yielded survival fractions of $\sim10^{-5}$ and $\sim10^{-7}$ respectively, if a considerable amount of additional energy is available to help unbind the envelope, the survival rates could be as high as $\sim10^{-2}$.

Like in our other tests, the population of surviving planets still skews towards massive planets and close-in orbits. However, this preference for massive planets is not surprising. We have derived an analytical expression (see Appendix) for the minimum planet mass threshold needed to survive CE: $\alpha\lambda M_{\mathrm{p}} \gtrsim20\,M_{\mathrm{Jup}}$. For the $\alpha=20$ case, it is possible for planets to survive with lower-masses and larger semi-major axes than in any of our other simulations.

While this high value of $\alpha$ is not needed to explain the population of stellar and brown dwarf companions in our galaxy [70], it remains unclear whether the same efficiencies apply in planet-star interactions, so it is not out of the question that a considerable amount of additional energy can be injected in planetary-mass CE events.

Discussion

Occurrence Rates of CE survivors

Our Monte Carlo approach allows us to predict the population of CE planets, rather than just the analytical parameters needed for a planet to survive CE. In many of our simulations, we noted that the fraction of planets that undergo and survive CE evolution is small. Here, we quantify this result by calculating the predicted occurrence rates of post-CE planets around white dwarfs.

Since every white dwarf originates from a main-sequence progenitor, we can calculate the predicted occurrence rate of common envelope survivors by multiplying the fraction of surviving planets in our simulations ($N_{\mathrm{surv}}/10^7$) by the occurrence rate of planets around main sequence stars. We source the occurrence rate of giant planets around main sequence stars (0.339 planets per star) from B. J. Fulton et al. [14], which was calculated for planets in a sample of stars with masses similar to the Sun. We also account for the dependence on host star mass [15], so our final distributions are scaled by:

where $M_\star$ is the host star mass of 1 or $2M_\odot$. We report the results of this calculation, the occurrence rate of post-CE planets per white dwarf, in Table 1.

The number of survivors primarily depends on the value of the common envelope efficiency $\alpha$ that we use. For the observationally-motivated value of $\alpha=0.41$ around a $2M_\odot$ WD progenitor, our simulations predict that we should find almost no CE-survivors; which indicates an upper limit of one survivor for every $10^7$ white dwarfs. For $\alpha=2$, we should find roughly one CE-survivor planet for every $10^4$ white dwarfs. In the scenario in which we tested the most additional energy injected, $\alpha=20$, we predict one planetary CE survivor out of every $10^2$ white dwarfs.

Comparison to Known Planets

One of the key motivations for this work was to shed light on how planets like WD 1856b and the candidate WD 0310b come to arrive in their present-day orbits.²

One way we can better understand the origin of these planets is to compare their properties (in particular, mass and orbital separation) to theoretical predictions of planets that undergo various formation pathways. Here, we review the properties of WD 1856b and WD 0310b and compare them to the outcomes of our CE simulations.

WD 1856b orbits its white dwarf host every 1.4 days, and has an orbital separation of about 0.02 AU (A. Vanderburg et al. [61]). The mass of WD 1856b is likely about 5 Jupiter masses ($M_{\rm p} = 5.2^{+0.7}_{-0.8}\,M_{\rm Jup}$), but could be as low as 0.84 $M_{\rm Jup}$ ([33] (S. Xu et al. 2021)). Meanwhile, if real, the candidate planet WD 0310b has a semi-major axis between 0.1 and 2 AU, and a mass between 0.2 and 3 $M_{\rm Jup}$.

Immediately, the masses and orbits of these two planets stand apart from the population of CE survivors predicted in this work. In our simulations, we found that typical survivors of common envelope evolution have masses $\gtrsim10\,M_{\rm Jup}$ and semi-major axes $\lesssim0.01$ AU.

Both WD 1856b and WD 0310b are less massive and have larger orbital separations than the typical common envelope survivor. The physical reason why common envelope survivors are dominated by closely-orbiting and massive planets is because these planets can release the most gravitational potential energy to unbind the stellar envelope. We therefore compared the mass, semi-major axis, and maximum amount of gravitational potential energy released ($U_{\max} = GM_{\rm c}M_{\rm p}/a_{\rm f}$) for the known planets and our simulation results in Figure 13.

The known closely-orbiting white dwarf planets are highly discrepant from the results of all of our simulations. None of our simulations produce planets with as large semi-major axes and low masses as WD 0310b, and only the most extreme $\alpha= 20$ case with the most additional energy used to eject the stellar envelope can match the parameters of WD 1856b. Even in this case WD 1856b occupies the lowest-energy tail of the distribution of common envelope survivors. If common envelope evolution were operating efficiently enough to form these two planets, it seems unlikely that two early discoveries would occupy such sparsely-populated regions of parameter space. Instead, it seems more likely that another formation mechanism is responsible for these objects. WD 1856b’s orbit, in particular, is highly consistent with the population of hot Jupiters around main sequence stars (see Figure 10 of S. Xu et al. 2021). We show this in the context of transiting hot Jupiters (planet mass between $0.25$-$13M_{\mathrm{Jup}}$ and orbital period $<10$ days) and our CE population around a $2M_{\odot}$ main sequence host run with different values of $\alpha$ in Figure 14.

We can also test whether the occurrence rate for CE survivors we predict is compatible with the presence of WD 1856b and WD 0310b in the solar neighborhood. Both systems are extremely close to the Sun in cosmic terms. WD 1856 is among the 250 closest white dwarfs within 40 pc of the Sun, while WD 0310 ranks among the 25 closest (M. W. O’Brien et al. [44]). We found that in most of our simulations, CE survivors are rare; our most realistic simulations predict occurrence rates less than $10^{-7}$. It is therefore very unlikely that two examples of such a rare process would be found so nearby.

Population of Transiting Survivors

The populations and occurrence rates we have predicted in this work do not account for the probability that the objects will be detected in observations. Because WD 1856b, the most confident example of a closely orbiting planet around a white dwarf, was discovered via transits, we simulate how the requirement that detected planets transit their white dwarfs affects our predicted population.

To correct for this bias, we can select systems based on their transit probability:

where $R_{\mathrm{p}}$ is a planet radius, $R_{\mathrm{WD}}$ is the radius of the white dwarf in our simulations ($0.015\,R_{\odot}$ for a $2M_{\odot}$ progenitor), and $a_{\mathrm{f}}$ is the final semi-major axis of each planet. For each planet, we draw a random number from a uniform distribution between 0 and 1. If the random draw is $\le P_{\mathrm{tr}}$, then we include the planet in the transiting sample. We summarize and compare the final distributions of transiting CE survivors with different $\alpha$ values in the bottom panels of Figure 13. We collect the number of transiting planets per white dwarf for each of our simulation runs and present them in Table 1.

We find that the outcome of this calculation is that the population of transiting CE survivor planets is even more strongly biased towards small semi-major axes, making the relatively wide orbits of WD 1856b and WD 0310b stand out even more. Another consequence, as expected, is that only a small fraction of planets transit their stars, so transiting CE survivors are very rare. We predict that the fraction of white dwarfs with transiting post-CE planets is $\sim10^{-5}$. This fraction is small enough that existing observations might not be expected to detect any such objects. For example, one of the main observatories that has been used to search for transiting planets around white dwarfs (and that discovered WD 1856b) is the Transiting Exoplanet Survey Satellite (TESS, G. R. Ricker et al. [57]), which has looked at an order of $\sim10^{4}$ WDs. If the fraction of white dwarfs with transiting CE survivors is really as low as our simulations predict, TESS has not looked at enough white dwarfs to likely detect such a system.

Fortunately, other datasets exist that provide an opportunity to detect CE survivors. In particular, there are ongoing efforts (e.g. Householder in prep.) to search through ground-based surveys like the Zwicky Transient Facility (ZTF, E. C. Bellm et al. [1]). Rich datasets like ZTF or the future Legacy Survey of Space and Time (LSST) conducted by the Vera C. Rubin Observatory (Ž. Ivezić et al. [24]), which have and will observe billions of objects, present an opportunity to detect very rare objects like post-CE planetary systems.

Future Work

This work is the first attempt to model a full post-CE planet population. However, there is more physical complexity that can be implemented into our approach to build a more complete picture of CE evolution and its outcomes for the formation of WD-planet systems. The first avenue for improvement is to consider alternative CE outcomes beyond the canonical ejection scenario. For example, scenarios like the one studied in L. Chamandy et al. [6] in which one planetary companion plunges into the stellar core and may deposit sufficient energy or angular momentum to help the other companion survive. These multi-body interactions could significantly affect the post-CE architecture and merit targeted simulations. Second, our treatment of the envelope’s internal energy has so far relied on simplified prescriptions. A more detailed approach, following methods such as those in F. Lagos et al. [30], may improve estimates of the binding energy and shift the predicted survival boundary. Third, we note that our tidal evolution modeling has relied on two prescriptions and a single planet in the system. Additional work is needed to test the sensitivity of the planet’s final orbit to different tidal formulations (J. Nordhaus & D. S. [41] (Spiegel 2013), including self-consistent dynamical-tide frameworks that couple stellar evolution to the full linear tidal response [40], which can shorten engulfment timescales by factors of several relative to equilibrium tides. Furthermore, additional planetary companions could be included in these simulations, potentially resulting in different system architectures. These alternatives could shift the range of parameter space in which outer planets remain bound and dynamically stable.

Lastly, our simulations have not considered accretion of the host star onto the planet and what that means for survival or detectability. Recent work by S. Jenkins et al. [26] simulated planet engulfment, concluding that the planet can accrete 3–48% of its initial mass during this stage at an accretion efficiency of 1.0. This can be a significant increase which would allow initially lower-mass planets to become more massive and survive CE evolution. The cited work also predicts the observational signatures that can flag such accretion events.

Future detections of wide-orbit planets around compact binaries via Gaia astrometry, direct imaging, or transit timing could directly test these survival scenarios and calibrate key ingredients of common-envelope evolution and tidal models.

Conclusion

In this work, we have simulated planetary systems via stellar evolution to understand the population of planets that survive common-envelope evolution. We explored a grid of simulation input parameters, varying assumptions such as host mass, tidal prescription, and envelope ejection efficiency. Our conclusions are as follows:

1. Planets that survive CE evolution are rare. We tested $1M_\odot$ and $2M_\odot$ host star masses and static or variable tidal prescriptions for 10 million systems of each variation. We find that most planets fail to unbind the envelope before being destroyed. Our simulations for the most realistic values of $\alpha=0.41$ or 1, based on population-wide estimates, result in a survival fraction of $\lesssim10^{-7}$. In systems where we include equal envelope binding energy to gravitational energy ($\alpha=2$), the survival fraction is $10^{-5}$. Finally, in the most optimistic case we considered ($\alpha=20$), the fraction of survivors is $10^{-3}$ to $10^{-2}$.

2. CE survivors are preferentially massive, typically $\gtrsim10M_{\mathrm{Jup}}$. High-mass planets are more effective at transferring energy and angular momentum during inspiral. Lower-mass planets rarely halt inspiral in time, suggesting that CE survival selects for the most massive end of the planet distribution.

3. CE survivors end up in close-in orbits (within 1.4 tidal disruption radii) around the resulting white dwarf. This is in contrast to planets that undergo high-eccentricity migration, which typically end up in orbits $>2a_{\mathrm{disrupt}}$ ([13]). The very close orbits of post-CE planets are therefore a distinct observational signature that

could help differentiate them from planets that migrate via other mechanisms. The properties of the known planets around white dwarfs WD 1856b and WD 0310b stand at odds with our simulated CE outcomes. Their orbital separation and masses lie outside the parameter space of our typical survivors. While CE evolution can explain the orbits of these planets (when additional energy sources are included), if this mechanism were operating we would expect to have also seen a larger and easier to detect population of massive planets with even smaller orbital separations. Since no such ultra-closely-orbiting massive planets around white dwarfs are yet known, our results suggest that the majority of close-in planets around white dwarfs arise from mechanisms other than CE evolution, such as high eccentricity migration.

Acknowledgments

We thank Hailiang Chen and Smadar Naoz for helpful conversations. This research has made use of NASA’s Astrophysics Data System Bibliographic Services. The authors acknowledge the MIT Office of Research Computing and Data for providing high performance computing resources that have contributed to the research results reported within this paper.

Software: numpy ([20]), scipy ([63]), pandas ([36]), matplotlib ([23]), MESA ([48, 49, 50, 51, 52, 27]).

Appendix

We can derive an approximate expression to estimate whether or not a planet of mass $M_p$ will survive given a scenario with specific $\alpha,\lambda$ values, under the assumption that survival is determined primarily by whether the planet survives tidal disruption (as is typically the case for planets more massive than about 11 $M_{\mathrm{Jup}}$ in our simulations), rather than photoevaporation. Start with equation (13) in the paper:

Because $M_{\mathrm{env}} \gg M_p$ we can neglect the “1”. This leads to

Now, we make the assumption that the planet survives if it avoids tidal disruption. From equation (4) of S. Rappaport et al. [54] we can express the minimum separation before Roche-lobe overflow of a planet by

where $M_c$ is the mass of the white-dwarf progenitor, and $\rho_p$ is the mean density of the planet in $g\ \mathrm{cc}^{-1}$. If we substitute this result into equation (A.2), and solve for $\alpha\lambda$ we find:

where we now call the mass of the giant $M_g$ (same as $M_\star$, in $M_\odot$), we have replaced $a_{\mathrm{engulf}}$ with the radius of the giant, $R_g$ in $R_\odot$, when it starts to overflow its RL, and we have approximated $r_L$ as unity. We can do this because the mass ratio $(M_g/M_p)$ is so large that $R_L \simeq a_i$.

The minimum value of $\alpha\lambda$ comes when $R_g$ is a maximum, which we can take to be $\sim$400 $R_\odot$ for a 1–3 $M_\odot$ star on the AGB. If we take $M_g$ to be 1.5 $M_\odot$ as representative, and $M_c \simeq0.6M_\odot$ as a typical WD mass. This comes to

where we have taken the mean density of gas giants to be near unity. This leads to a simple analytic approximation for a planet to successfully eject the common envelop of its host star:

This result is broadly consistent with our findings using more detailed calculations, and provides an analytic estimate of the planetary masses we can expect to survive the common envelope. Of course, the planet must also survive the heat bath of both the red giant interior and the nascent white dwarf if it is to last to the current epoch.