TATOOINE'S FUTURE: THE ECCENTRIC RESPONSE OF KEPLER'S CIRCUMBINARY PLANETS TO COMMON-ENVELOPE EVOLUTION OF THEIR HOST STARS

"Everything is in motion," Plato quotes Heraclitus (Cratylus, Paragraph 402, section a, line 8), "and nothing remains still." Despite the overwhelming timescales, from a human perspective, the natural life cycle of stars is a prime example of change on the cosmic stage. Following the laws of stellar astrophysics, over millions to billions of years, stars and stellar systems form, evolve, and ultimately die (e.g., Kippenhahn & Wiegert 1990 and references therein). And so do planetary systems, with their fate linked intimately to that of their stellar hosts.

The fate of planets orbiting single stars has been studied extensively, indicating that planets can survive their star's evolution if they avoid engulfment or evaporation during the red giant branch (RGB) and the asymptotic giant branch (AGB) stages (e.g., Livio & Soker 1984; Rasio & Livio 1996; Duncan & Lissauer 1998; Villaver & Livio 2007; also see Veras et al. 2016). There is also accumulating evidence of planetary or asteroidal debris surrounding white dwarfs and polluting their atmospheres (see recent reviews by Farihi 2016 and Veras 2016). A planet's survival depends both on the initial orbital separation and on its mass. For example, even a planet as massive as 15 M_Jup around a 1 ${M}_{\odot }$ main-sequence star can be destroyed inside the stellar envelope during the AGB phase (Villaver & Livio 2007). An unpleasant prospect for the future of our own planet is that it might not survive the evolving Sun (e.g., Schroder & Connon Smith 2008). Despite the complications, theoretical considerations indicate that planets can indeed survive the evolution of single stars (Villaver & Livio 2007), and there is mounting observational evidence supporting this (e.g., Reffert et al. 2015; Wittenmyer et al. 2016, and references therein⁷ ).

Single stars, however, do not have a monopoly over planetary systems. Nearly half of solar-type stars are members of binary and higher order stellar systems (Raghavan et al. 2010), and an increasing number of planets have been discovered in such systems. While the presence of a distant stellar companion will have little effect on a planet around one (evolving) member of a wide binary stellar system (≳100 au), a planet orbiting around both members of a close binary system (separation ∼10 au or less) will experience qualitatively and quantitatively different stellar evolution. Namely, where a single star loses mass and expands on timescales of millions to billions of years (Veras 2016) while on the main sequence (MS), the RGB, and the AGB—and its planets react accordingly—close binary stars can experience events of dramatic mass loss, orbital shrinkage on timescales of years, and common-envelope (CE; see Appendix for abbreviations and parameters) stages where the two stars share (and quickly expel) a common atmosphere (Paczynski 1976; Hilditch 2001). During the in-spiraling CE phase, the two stars can strongly interact with each other by transferring mass, smoothly coalesce or violently collide, or even explode as a supernova (SN). Both the amount of energy released during this stage and the timescale of the release are staggering: a close binary star can lose an entire solar mass over the course of just a few months (e.g., Livio & Soker 1984; Rasio & Livio 1996; Ricker & Taam 2008, 2012; Passy et al. 2012; Ivanova et al. 2013; Nandez et al. 2014, but also see Sandquist et al. 1998 and De Marco et al. 2009 for longer timescales); in the extreme case of an SN, the binary will be completely disrupted. Overall, the evolution of close binary stars is much richer than single stars, as the separation, eccentricity, and metallicity of the binary add additional complications to an already complex astrophysical process. Interestingly, to date there are nine close MS binary systems harboring confirmed circumbinary planets (CBPs; Doyle et al. 2011; Orosz et al. 2012a, 2012b; Welsh et al. 2012, 2015; Kostov et al. 2013, 2014; Schwamb et al. 2013). It is reasonable to assume that the reaction of these planets to the violent CE phase of their host binary stars will be no less dramatic.

Previous studies of the dynamical response of planets in evolving multiple stellar systems have focused on post-CE (PCE) CBP candidates (e.g., Portegies Zwart 2013; Mustill et al. 2013; Volschow et al. 2014), on circumprimary planets in evolving wide binary systems (Kratter & Perets 2012), and on the survival prospects of CB planets (Veras et al. 2011; Veras & Tout 2012). In particular, Portegies Zwart (2013) and Mustill et al. (2013) start with the final outcome of close binary evolution—PCE systems with candidate CBP companions (HU Aqu and NN Ser, respectively)—and reconstruct the initial configuration of their progenitors.

Inspired by recent discoveries of planets in multiple stellar systems, here we tackle the reverse problem: we start with the nine confirmed CBP Kepler systems with known initial binary configurations (using data from NASA's Kepler  mission), and we study their dynamical response as their stellar hosts evolve. We combine established, publicly available binary star evolution code (BSE; Hurley et al. 2002) with direct N-body integrations (using REBOUND, Rein & Liu 2012), allowing for stellar mass loss and orbital shrinkage on dynamical timescales. For simplicity, here we focus explicitly on dynamical interactions only (e.g., no stellar tides) and assume that the CBPs do not interact with the ejected stellar material. We expand on the latter limitation in Section 5.1.

We note that Mustill et al. (2013) assume an adiabatic mass-loss regime (Hadjidemetriou 1963) during the CE stage for the evolution of NN Ser. In this case, the orbital period of the CBP is much shorter than the CE mass-loss timescale (hereafter T_CE), and the planet's orbit expands at constant eccentricity.⁸ Volschow et al. (2014) study the two analytic extremes of mass-loss events for the case of NN Ser: the adiabatic regime and the instantaneous regime where ${T}_{\mathrm{CE}}$  is much shorter than the orbital period of the planets.⁹ The Kepler  CBP systems, however, occupy a unique parameter space where the planets' periods are comparable to the ${T}_{\mathrm{CE}}$  used by Portegies Zwart (2013) and suggested by Ricker & Taam (2008, 2012) for binary systems similar to Kepler  CBP hosts. Thus the adiabatic approximation may not be an adequate assumption when treating the CE phases of these systems. Specifically, the CE phase can occur on timescales as short as one year (Paczynski 1976; Passy et al. 2012), and during this year a binary can expel ∼2 ${M}_{\odot }$ worth of mass (Ricker & Taam 2008, 2012; Iaconi et al. 2016). Interestingly, such timescales are well within a single orbit of CBPs like Kepler-1647 with P_CBP = 3 years (Kostov et al. 2016). Such a dramatic restructuring will have a profound impact on the dynamics of the system—for example, a CBP can achieve very high eccentricity—and may even result in the ejection of the planet. Similar to Portegies Zwart (2013), here we tackle such complications by directly integrating the equations of motion during the CE phases for each Kepler  CBP system.

The transition from adiabatic to nonadiabatic regime can be characterized through a mass-loss index ψ. For a planet orbiting a single star, Veras & Tout (2012) define

$$\begin{eqnarray}{\rm{\Psi }} & \equiv & \displaystyle \frac{{\alpha }_{M}\,}{n\mu }={(2\pi )}^{-1}\left(\displaystyle \frac{{\alpha }_{M}}{1{M}_{\odot }\,{\mathrm{yr}}^{-1}}\right)\\ & & \times {\left(\displaystyle \frac{{a}_{p,0}}{1\mathrm{au}}\right)}^{\tfrac{3}{2}}{\left(\displaystyle \frac{\mu }{1{M}_{\odot }}\right)}^{-\tfrac{3}{2}},\end{eqnarray} \tag{ 1 }$$

 where ${\alpha }_{M}$  is the mass-loss rate, a_p,0 is the initial semimajor axis of the planet, and μ = M_star + M_p is the total mass of the system.¹⁰ If Ψ ≪ 1, then the evolution of the planet's orbit is in the adiabatic regime and its semimajor axis grows at constant eccentricity.

Alternatively, if Ψ ≫ 1, then the planet sees an "instantaneous" mass loss, and its orbital evolution is in a "runaway" regime. In this scenario, the fate of the planet depends on the ratio between the final and initial mass of the system, that is, β = μ_final/μ_init, and on the orbital phase of the planet at the onset of the mass-loss event. For example, if a highly eccentric planet is at pericenter at the beginning of the CE stage, it becomes unbound. Overall, ejection occurs when the system loses sufficient mass by the end of the mass-loss event. To first order, the critical mass ratio for ejection β_eject is given by Equation (39) of Veras et al. (2011), which we reproduce here for completeness:

$$\begin{eqnarray}&&{\beta }_{\mathrm{eject}}\equiv 0.5(1+{e}_{{\rm{p}},0})\end{eqnarray} \tag{ 2 }$$

where e_p,0 is the eccentricity of the planet at the beginning of the mass-loss event. A planet on a circular orbit or at pericenter is ejected from the system if β < β_eject.¹¹ A planet at apocenter remains bound even in the runaway regime, and its orbit expands to a_{runaway,apocenter} and circularizes or reaches an eccentricity of e_(post-circ). If a planet remains bound in the runaway regime, its semimajor axis is given by Equations (43) and (44) of Veras et al. (2011):

$$\begin{eqnarray}&&{a}_{\mathrm{runaway}}\equiv \displaystyle \frac{{a}_{{\rm{p}},0}(1\mp {e}_{{\rm{p}},0})}{2-\beta (1\pm {e}_{{\rm{p}},0})},\end{eqnarray} \tag{ 3 }$$

where the signs represent initial pericenter and apocenter, respectively.

According to the above prescription, the dynamical evolution of a CBP's orbit for Ψ ∼ 0.1–1 lies between these two regimes. For both Ψ ∼ 1 and Ψ ≫ 1, the evolution is not adiabatic, and e_CBP can vary (increase or decrease) as a_CBP  varies. Thus the planet may become unbound only if Ψ ≳ 1, with the caveat that the transition is not clear-cut and the Ψ ∼ 0.1–1 regime is more complex (Veras et al. 2011; Veras & Tout 2012; Portegies Zwart 2013; Veras 2016). As we show below, the fact that the central object in a CBP system is itself a binary instead of a single star adds yet another layer of complication to these theoretical considerations.

Based on Equation (1), and assuming a CE timescale of ${T}_{\mathrm{CE}}$  ∼ 10³ years, we would a priori expect that all Kepler CBPs remain bound to their hosts during their primary stars' CE stage (e.g., see Figure 3, Veras & Tout 2012). However, for the rapid CE timescale of ${T}_{\mathrm{CE}}$  ∼ 1 year suggested by recent results (e.g., Ricker & Taam 2008, 2012), some of the Kepler  CBP systems are within a factor of 5–10 of the transition to the nonadiabatic regime during the primary or secondary CE (i.e., Ψ ∼ 0.1–0.2), and within a factor of 2 of the transition during their respective secondary CE stage (i.e., Ψ ∼ 0.5). We might therefore expect significant variations in the respective CBPs' orbital eccentricities, as well as deviations from the final semimajor axes expected from adiabatic mass loss.

Because the transition between adiabatic and nonadiabatic regimes is difficult to study analytically, here we examine the dynamical evolution of all Kepler CBP systems numerically, using the adaptive-time-step, high-order integrator IAS15 (Rein & Spiegel 2015), which is available as part of the modular and open-source REBOUND package (Rein & Liu 2012). REBOUND is written in C99 and comes with an extensive Python interface.

This paper is organized as follows. In Section 2 we describe the algorithm we use for the evolution of the Kepler  CBP hosts. Section 3 details our implementation of an N-body code in exploring the dynamical reaction of said CBPs to the respective CE stages. We present the results for each system in Section 4, and we discuss our results and draw conclusions in Section 5.

The evolution of a gravitationally bound system of two stars depends primarily on its initial configuration. If the binary is wide enough, the stars will evolve in isolation, according to the prescription of single-star evolution theory. Alternatively, if the initial separation is sufficiently small, then the stars will influence each other's evolution. The evolution of such systems, described as close (or interacting) binary stars (Hilditch 2001), is more complex and includes a number of additional processes, such as tidal interactions and mass transfer. The nine Kepler CBP-harboring binary systems fall in the latter category, and we use the established, open-source binary star evolution code (Hurley et al. 2002) to study their evolution. Briefly, the code works as follows.

While the binary is in a detached state, the code evolves each star individually using a single-star evolutionary code (Hurley et al. 2000), which includes tidal and braking mechanisms and wind accretion. When the stars begin to interact, BSE uses a suite of binary-specific features, such as mass transfer and accretion, CE evolution, collisions and mergers, and angular momentum loss mechanisms. The evolution algorithm allows the specific processes to be turned on and off, or even modified based on custom requirements. By default, all stars are assumed to be initially on the zero-age main sequence (ZAMS), but any possible evolutionary state, such as corotation with the orbit, can be the starting point of the simulation. BSE uses an evolution time step small enough to prevent the stellar mass and radius from changing too much (not more than 1% in mass and 10% in radius), which allows identification of the time when, and if, a star first fills its Roche lobe.

At the onset of Roche-lobe overflow (RLOF), the two stars either come into contact and coalesce or initiate a CE stage. CE evolution is a complex mechanism that typically occurs when an evolved, giant star transfers mass to an MS star on a dynamical timescale.¹² When the giant star overfills the Roche lobe of both stars, its core and the MS star share a single envelope. This envelope rotates slower than the orbiting "cores" within due to expansion, and the resulting friction causes in-spiraling and energy transfer to the envelope, described by an efficiency parameter α_CE. The outcome of this phase is either ejection of the envelope (assumed to be isotropic) if neither core fills its Roche lobe, leaving a close white dwarf–MS pair (assumed to be in corotation with the orbit), or coalescence of the cores. A challenging process to study numerically, BSE recognizes the CE phase—beginning with RLOF—and forces the system through it on an instantaneous timescale, allowing the evolution to continue according to the appropriate PCE parameters. The code resolves the outcome of the CE phase based on the initial binding energy of the envelope and on the initial orbital energy of the two cores.¹³

On input, BSE requires a number of parameters. Namely, the CE efficiency parameter α_CE (in the range of 0.5–10), the masses and stellar types of both the primary and secondary stars (in the range of 0.1–100 ${M}_{\odot }$), the binary orbital period (in days), the eccentricity (0.0–1.0), the maximum time for the evolution (here we choose 15 Gyr, roughly the age of the universe, unless a particular Kepler CBP system stops evolving—or the stars coalesce—much sooner, at which point we stop the simulations), and the metallicity (0.0001–0.03) of the system. The code also checks whether tidal evolution is included, noted as "on" (hereafter TCP for tidal circularization path) or "off" (NTCP for no tidal circularization path). Given the still-uncertain mechanism of the CE stage, the parameters describing the efficiency of envelope ejection and the treatment of tidal decay—and their associated uncertainties—are likely the dominant source of error in our BSE results.

Depending on α_CE, a binary system may exit a CE stage as a merged star, as a very tight PCE binary star, or may even trigger a supernova explosion. For example, a higher α_CE  accounts for energy sources other than orbital energy, and the system then requires less energy to dissipate the envelope. The BSE code also allows an alternate CE model to be used, the de Kool CE evolution model (de Kool 1990), which first introduced the CE evolution binding energy factor λ.

Likewise, tides are integral to the evolution of close binary stars. The strength of tidal dissipation will affect the orbital separation and eccentricity when mass transfer begins and hence the future evolution under RLOF or a CE stage. BSE achieves orbital circularization mainly through tides, but also by accounting for more mass accretion at pericenter than apocenter of an eccentric orbit. Another effect is the exchange of orbital angular momentum with the component spins due to tidal synchronization, thus causing the evolution to proceed in a similar way to closer binary systems if tides were ignored. In addition, primary spin–orbit corotation can be forced within the code to avoid unstable RLOF, which gives two possible outcomes for each system. To account for the effect of tides on binary evolution, we use both the tides "on" (TCP) or "off" (NTCP) options in BSE. As a result, the binary starts the first RLOF at either e_bin,RLOF = 0 and a_bin,RLOF < a_bin,0 (for TCP), or at e_bin,RLOF = e_bin,0 and a_bin,RLOF = a_bin,0 (for NTCP).

Additionally, while the binary periods of the Kepler  CBP systems are known to very high precision, their stellar masses and metallicities have nonnegligible uncertainties.¹⁴ For example, the primary masses for Kepler-38 and Kepler-64 have 1σ errors of 0.05 M_Sun and 0.1 M_Sun, respectively, and the 1σ errors on the metallicities of half of the systems are as large as, or even larger than, the measured values. To better characterize the impact of these uncertainties on the evolution of the Kepler  CBP systems, we obtain BSE results for both the best-fit mass and metallicity values and also for their corresponding 3σ range.¹⁵

The entire parameter space we explored, shown in Table 1, includes a range of CE parameters as well and consists of 22,680 BSE simulations. Our default configurations are denoted in boldface in the table; for the rest of the input parameters, we use the default values from model A of Hurley et al. (2002), described as the most favorable and effective. Model A uses all default wind-related parameters and a Reimers mass-loss coefficient of η = 0.5.

As BSE evolves the stars, their types change, and each type is checked throughout the code to ensure the correct evolutionary tracks are used for both the individual stars and the system as a whole. These numerical values, as well as what they represent, are listed in Table 2. At each significant event, such as at the beginning and end of the RLOF and of the CE phase, when the stars coalesce or collide, the code also produces a "type" label. These labels, along with their meaning, are shown in Table 3.

An example of the BSE output for Kepler-38, α_CE  = 0.5, NTCP, and TCP is shown in Table 4. The BSE results for all Kepler  systems can be found in an online supplement to Table 4.

Table 4.  BSE Evolution of Kepler-38 for α_CE  = 0.5

Timea	M1	M2	Stell. Type	Stell. Type	abin	ebin	R1/RRoche	R1/RRoche	Evol. Stage
(Gyr)	(${M}_{\odot }$)	(${M}_{\odot }$)			(R⊙)
Kepler-38b
0.0	0.95	0.25	1	0	31.6	0.10	0.05	0.03	INITIAL
12.36	0.95	0.25	2	0	31.6	0.10	0.10	0.03	KW CHNGE
13.01	0.95	0.25	3	0	31.61	0.10	0.14	0.03	KW CHNGE
13.72	0.94	0.25	3	0	31.86	0.10	1.00	0.03	BEG RCHE
13.72	0.75	0.25	3	15	0.51	0.00	1.00	0.03	COMENV
13.76	0.47	0.00	10	15	0.00	−1.00	0.00	−1.00	KW CHNGE
15.00	0.47	0.00	10	15	0.00	−1.00	0.00	−1.00	MAX TIME
Kepler-38c
0.0	0.95	0.25	1	0	31.59	0.10	0.05	0.03	INITIAL
12.36	0.95	0.25	2	0	31.60	0.10	0.10	0.03	KW CHNGE
13.01	0.95	0.25	3	0	31.45	0.09	0.14	0.03	KW CHNGE
13.70	0.94	0.25	3	0	24.58	0.00	1.00	0.04	BEG RCHE
13.70	0.91	0.25	3	15	0.37	0.00	1.00	0.04	COMENV
13.75	0.66	0.00	4	15	0.00	−1.00	0.00	−1.00	KW CHNGE
13.88	0.63	0.00	5	15	0.00	−1.00	0.00	−1.00	KW CHNGE
13.89	0.55	0.00	6	15	0.00	−1.00	0.00	−1.00	KW CHNGE
13.89	0.52	0.00	11	15	0.00	−1.00	0.00	−1.00	KW CHNGE
15.00	0.52	0.00	11	15	0.00	−1.00	0.00	−1.00	MAX TIME

Notes. The BSE results for our entire set of simulations for all Kepler systems are published in their entirety in an online, machine-readable format.

^aTime in the online supplement is in Myr. ^bNTCP. ^cTCP.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

In the following section we describe how we use the initial and BSE-generated final binary masses and orbital separations for each Kepler CBP system to explore the dynamical response of the planets during the respective CE stages.

While the concept of CE evolution for binary stars was proposed 40 years ago (Paczynski 1976), detailed understanding of this important stage is a challenging task that requires computationally intensive hydrodynamical calculations. The outcome of this complex phase has been studied both numerically and analytically (see Table 1 of Iaconi et al. 2016 for a list of references). For example, hydrodynamical simulations by Ricker & Taam (2012) for a 1.05 ${M}_{\odot }$+0.36 ${M}_{\odot }$ binary star with an initial orbital period of 44 days (similar to Kepler CBP hosts) show that in the later stages of the CE the system loses mass at a rate of 2 ${M}_{\odot }$ yr⁻¹, the binary orbital separation decreases by a factor of 7 after ∼56 days, and ∼90% of the outflow is contained within 30 degrees of the binary's orbital plane.

To numerically account for the dynamical impact on Kepler CBPs from such drastic orbital reconfigurations and tremendous, short-timescale mass-loss phases, we use REBOUNDx, a library for incorporating additional effects beyond point-mass gravity in REBOUND simulations (https://github.com/dtamayo/reboundx). In particular, we shrink the binary orbits with the REBOUNDx modify_orbits_forces routine, which adds to the equations of motion a simple drag force for opposite particles' velocity vectors (Papaloizou & Larwood 2000). When orbit-averaged over a Keplerian orbit, this approach results in an exponential damping of the semimajor axis where the e-folding timescale is set through the particle parameter τ_a. Hydrodynamical simulations by Passy et al. (2012) and Iaconi et al. (2016) exhibit similar orbital decay during the CE stage (e.g., see their respective Figures 5 and 14), validating our approach.

Given the complexities and uncertainties of the mass-loss mechanism during the CE phases, we adopt a simple analytical parameterization similar to the method of Portegies Zwart (2013), who tackled the problem by using a constant mass-loss rate. Here we choose to instead use the same functional forms for the CE-driven changes in M_star as we do for a_bin (utilizing the REBOUNDx modify_mass routine), that is, an exponential mass loss for particles on their assigned e-folding timescales τ_M. Namely, for a given CE stage where one of the stars transitions from an initial mass M₀ to a final mass M_f (both provided by BSE), the mass loss occurs over a time ${T}_{\mathrm{CE}}$  = (M_f − M₀)/${\alpha }_{M}$. We note that an exponential mass loss has been explored by Adams et al. (2013) and Adams & Bloch (2013) as well, but for the case of single-star systems.

To ensure a smooth mass evolution (and avoid possible numerical artifacts) that asymptotes at M_f, we linearly evolved REBOUNDx's exponential mass-loss rate so that it vanished at time ${T}_{\mathrm{CE}}$. In particular, in terms of the e-folding timescale τ_M,

$$\begin{eqnarray}&&\dot{M}=-\displaystyle \frac{M}{{\tau }_{M}},{\tau }_{M}=\displaystyle \frac{{\tau }_{M,0}}{(1-{\rm{t}}/{T}_{\mathrm{CE}}\,)}.\end{eqnarray} \tag{ 4 }$$

 This differential equation can be solved analytically, yielding

$$\begin{eqnarray}&&M(t)={M}_{f}{\left(\displaystyle \frac{{M}_{0}}{{M}_{f}}\right)}^{{(1-{\rm{t}}/{T}_{\mathrm{CE}})}^{2}}\end{eqnarray} \tag{ 5 }$$

as long as τ_M,0 in Equation (4) is chosen as

$$\begin{eqnarray}&&{\tau }_{M,0}=\displaystyle \frac{{T}_{\mathrm{CE}}\,}{2\mathrm{ln}\left(\tfrac{{M}_{0}}{{M}_{f}}\right)}.\end{eqnarray} \tag{ 6 }$$

As the star's mass decreases from M₀ to M_f during the respective CE stage, the binary orbit must shrink from a_bin,0 to a_bin,PCE (as provided by BSE) as well, and both must occur over a time ${T}_{\mathrm{CE}}$. We apply the same prescription for evolving a_bin as that for the mass loss (i.e., Equations (4)–(6) with τ_a instead of τ_M), but we caution that while this would evolve a_bin from a_bin,0 to a_bin,PCE in isolation, the mass loss is simultaneously increasing a_bin. To correct for this and achieve the desired a_bin,PCE, we manually adjust the corresponding τ_a,0 parameter.

While tides are not directly accounted for in our dynamical simulations as the particles involved are treated as point masses, we capture their effects in a simple parameterized way. Specifically, we impose the same exponential decay on the evolution of e_bin as we do for the binary's semimajor axis shrinking and stellar mass loss (i.e., Equations (4)–(6) but with τ_e). As a result, e_bin undergoes dampened oscillations during the CE phase, at the end of which it settles within a few percent of zero, as required by BSE.

As mentioned in the previous section, while the BSE code yields initial and final values for the stellar masses and a_bin at the beginning and end of the CE phases, it does not resolve the extremely rapid mass-loss phase. We bridge this gap by exploring two regimes of mass loss, namely ${\alpha }_{1}$ = 1.0 M_⊙ yr⁻¹ and ${\alpha }_{0.1}$ = 0.1 M_⊙ yr⁻¹, denoted by their respective subscripts. The former is based on the results of Ricker & Taam (2008, 2012) and Portegies Zwart (2013) for systems similar to the Kepler  CBP hosts. Values of ${\alpha }_{0.1}$ (and smaller) typically guarantee that the orbital evolution of these CBPs is in the adiabatic regime (see Equation (1)).

To test the applicability of our numerical simulations to the stated problem, we examine the evolution of a mock Kepler-38-like system where the binary is replaced with a single star of the same total initial and final mass according to the BSE prescription for Kepler-38 (see Table 4) and using the same initial orbit of the planet. Here we test two mass-loss regimes: runaway (for ${\alpha }_{M}$  = 100.0 ${M}_{\odot }$ yr⁻¹) and adiabatic (for ${\alpha }_{M}$  = 0.1 ${M}_{\odot }$ yr⁻¹ ≡ ${\alpha }_{0.1}$). The results are shown in Figure 1, where the three panels represent (from upper to lower, respectively) the evolution of the stellar mass, of the planet's semimajor axis, and of the planet's eccentricity as a function of time (such that mass loss starts at t = 0 and ends at t = T_mass-loss). According to the analytic prescription, the planet should remain bound even in the runaway mass-loss regime as the mass lost is smaller than the critical mass limit for ejection, as indicated by the dotted line in the upper panel of the figure.¹⁶ The horizontal dashed lines in the middle and upper panels represent the theoretical adiabatic (green) and runaway (red) approximations. As seen from these two panels, the numerical simulations are fully consistent with the theoretical approximations, demonstrating the validity and applicability of our simulations.

Zoom In Zoom Out Reset image size

In order to interpret the results from numerical integrations, it is useful to define a critical mass-loss rate, α_crit, corresponding to Ψ = 1 (see Equation (1)). CBPs in systems with α_M ≫ α_crit  will experience "instantaneous" mass loss, while the orbits of those in systems with α_M ≪ α_crit  should evolve adiabatically. A CBP in a system where α_M becomes comparable to or larger than α_crit  at any point during the binary evolution risks becoming unbound. Secondary CE stages will have two different α_crit  values as the respective CBPs will have two different values for a_CBP  and e_CBP (but the same μ) at the end of the preceding CE stages (see Equation (1)): one for ${\alpha }_{1}$ and another for ${\alpha }_{0.1}$. We denote these critical values as α_crit (α₁) and α_crit (α_0.1). In the subsections below, for each CBP Kepler  system, we compare α_crit  to α₁ and α_0.1 for each CE phase. Where the adopted ${\alpha }_{1}$ is within a factor of 10 of the respective α_crit, we pay special attention to the CBP orbit and test its evolution for 0.1 ${\alpha }_{1}$, 0.2 ${\alpha }_{1}$, and 0.5 ${\alpha }_{1}$ as well. Unless indicated otherwise, the final fate of the CBP for these cases is similar to the case for ${\alpha }_{1}$.

The description of a three-body system containing an evolving binary star is inherently multidimensional, and minor changes in one parameter can cascade into major changes in the rest. In addition to the importance of α_CE, tides, and α_M, another key issue—as shown in the next section—is the phase offset between the binary star and the CBP at the onset of the CE phase, Δθ₀ ≡ θ_CBP,0 − θ_bin,0 − ω_bin (i.e., the difference between their true anomalies, taking into account the binary's argument of the pericenter).

Thus for each system we test four initial binary phases corresponding to the binary orbit turning points, that is, eastern and western elongations (EE and WE), and superior and inferior conjunctions (SC and IC), and 50 initial CBP phases (between 0 and 1, with a step of 0.02) for a total of 200 initial conditions. We explore each of these initial conditions for α_CE = (0.5, 1.0, 3.0, 5.0, 10.0), for TCP/NTCP, and for ${\alpha }_{1}$ and ${\alpha }_{0.1}$.

Finally, we note that the post-CE semimajor axis of a CBP (a_CBP,PCE) that was initially on an eccentric orbit (i.e., e_CBP,0 > 0) depends on its orbital phase at the onset of the CE stage. For those Kepler  systems that undergo a secondary CE phase with mass loss, we start the respective integrations at the modes of a_CBP,PCE and e_CBP,PCE from the preceding CE stage.

Here we describe the results of our simulations for those Kepler  CBP systems that experience at least one CE phase. Unless specified otherwise, the evolution of each system is for the default CE parameters listed in Table 1. The initial orbital parameters for each system are listed in Table 5. We note that our N-body integrations only cover the CE phases.

Table 5.  Kepler System Default Initial Parameters

Name	Kepler ID	M1	M2	Pbin	ebin	Z	PCBP	eCBP
		(${M}_{\odot }$)	(${M}_{\odot }$)	(days)			(days)
Kepler-16	12644769	0.69	0.20	41.08	0.16	0.010	228.78	0.01
Kepler-34	8572936	1.05	1.0208	27.79	0.52	0.016	288.82	0.18
Kepler-35	9837578	0.89	0.81	20.73	0.14	0.009	131.46	0.04
Kepler-38	6762829	0.95	0.25	18.79	0.1	0.015	105.59	0.03
Kepler-47a	10020423	0.96	0.34	7.45	0.02	0.011	49.5	0.03 (P1)
							187.4	0.02 (P2)
							303.16	0.04 (P3)
Kepler-64	4862625	1.53	0.41	20	0.22	0.031b	138.5	0.1
Kepler-413	12351927	0.82	0.54	10.12	0.04	0.012	66.26	0.12
Kepler-453	9632895	0.94	0.2	27.34	0.05	0.023	240.5	0.04
Kepler-1647	5473556	1.22	0.97	11.259	0.15	0.014	1107.59	0.06

Notes.

^aJ. Orosz (2016, private communication). ^bBSE allows maximum metallicity of Z = 0.03.

Download table as:  ASCIITypeset image

Below we outline the major binary evolution stages (as produced by BSE), along with the dynamical responses of the CBPs; the main results are presented in Tables 6–10. All reported times are in Gyr since ZAMS. For simplicity, we assume that the current Kepler  CBP systems are at ZAMS and that all CBPs start on initially circular, coplanar orbits.

Table 7.  PCE CBP Semimajor Axis and Eccentricity for ${\alpha }_{M}$  = 1.0 M_Sun yr⁻¹ and for Primary CE

αCE	Tides		aCBP,PCE (au)			eCBP,PCE		Notes
		Mode	Median	Range	Mode	Median	Range
Kepler-38								[aCBP,0 = 0.46]
0.5	NTCP	${0.6}_{-0.01}^{+0.2}$	0.8	0.58–1.7	${0.15}_{-0.07}^{+0.12}$	0.26	0.05–0.67	⋯
0.5a	TCP	${0.5}_{-0.01}^{+0.4}$	0.6	0.49–0.9	${0.4}_{-0.35}^{+0.01}$	0.2	0.03–0.42	⋯
1/3/5/10	NTCP	${1.3}_{-0.2}^{+0.1}$	1.3	0.9–3.8	${0.38}_{-0.01}^{+0.16}$	0.42	0.02–0.82	⋯
1/3/5/10	TCP	${1.4}_{-0.2}^{+0.01}$	1.4	1.1–1.7	${0.46}_{-0.08}^{+0.01}$	0.42	0.28–0.56	⋯
Kepler-47
0.5, P1b	NTCP	${0.31}_{-0.01}^{+0.09}$	0.34	0.3–0.44	${0.12}_{-0.03}^{+0.04}$	0.14	0.06–0.27	[aCBP,0 = 0.29]
0.5, P2b	NTCP	${0.66}_{-0.01}^{+0.42}$	0.88	0.61–1.67	${0.14}_{-0.05}^{+0.35}$	0.22	0.06–0.56	[aCBP,0 = 0.70]
0.5, P3b	NTCP	${0.9}_{-0.01}^{+0.61}$	1.2	0.8–2.87	${0.17}_{-0.05}^{+0.29}$	0.23	0.09–0.66	[aCBP,0 = 0.96]
0.5a, P1b	TCP	${0.32}_{-0.01}^{+0.07}$	0.34	0.31–0.39	${0.1}_{-0.02}^{+0.08}$	0.1	0.07–0.18	⋯
0.5a, P2b	TCP	${0.7}_{-0.01}^{+0.5}$	0.87	0.67–1.24	${0.08}_{-0.01}^{+0.32}$	0.21	0.06–0.42	⋯
0.5a, P3b	TCP	${0.93}_{-0.01}^{+0.91}$	1.18	0.88–1.89	${0.11}_{-0.01}^{+0.35}$	0.22	0.09–0.48	⋯
1.0, P1b	NTCP	${0.49}_{-0.02}^{+0.01}$	0.49	0.47–0.5	${0.08}_{-0.04}^{+0.05}$	0.08	0.04–0.13	⋯
1.0, P2b	NTCP	${1.33}_{-0.17}^{+0.17}$	1.41	1.12–1.97	${0.42}_{-0.14}^{+0.11}$	0.4	0.22–0.58	⋯
1.0, P3b	NTCP	${1.65}_{-0.01}^{+0.53}$	2.13	1.52–4.15	${0.5}_{-0.14}^{+0.18}$	0.51	0.29–0.75	⋯
1.0a, P1	TCP	${0.38}_{-0.01}^{+0.03}$	0.39	0.37–0.41	${0.11}_{-0.1}^{+0.01}$	0.07	0.01–0.12	⋯
1.0a, P2	TCP	${0.91}_{-0.01}^{+0.29}$	1.03	0.89–1.22	${0.36}_{-0.21}^{+0.01}$	0.26	0.14–0.38	⋯
1.0a, P3	TCP	${1.23}_{-0.01}^{+0.47}$	1.46	1.2–1.83	${0.44}_{-0.26}^{+0.01}$	0.32	0.17–0.45	⋯
3/5/10, P1d	NTCP	${0.67}_{-0.01}^{+0.01}$	0.67	0.66–0.68	${0.11}_{-0.01}^{+0.03}$	0.11	0.04–0.17	⋯
3/5/10, P2d	NTCP	${3.23}_{-0.01}^{+0.57}$	3.4	2.23–7.95	${0.71}_{-0.06}^{+0.1}$	0.73	0.57–0.89	⋯
3/5/10, P3d	NTCP	${7.64}_{-2.71}^{+7.21}$	9.97	4.48–22.52	${0.92}_{-0.02}^{+0.02}$	0.89	0.74–0.95c	46% ejected
3/5/10, P1d	TCP	${0.69}_{-0.01}^{+0.01}$	0.69	0.68–0.7	${0.1}_{-0.01}^{+0.03}$	0.11	0.08–0.14	⋯
3/5/10, P2d	TCP	${3.68}_{-0.23}^{+0.91}$	4.18	3.34–5.62	${0.75}_{-0.01}^{+0.07}$	0.78	0.72–0.83	⋯
3/5/10, P3d	TCP	${15.69}_{-0.01}^{+6.16}$	19.7	15.3–23.01	${0.93}_{-0.01}^{+0.02}$	0.94	0.93–0.95c	81% ejected

Notes. The modes are shown with their respective 68% range. We caution that most of the probability distribution functions (PDFs) are notably nonnormal, often with two dominant peaks at the range limits.

^aSinusoidal distributions in phase versus final a_CBP/e_CBP (see, e.g., Figure 5), typically with double-peaked PDFs. ^bP1/P2/P3 have 40%/15%/75% probability to become dynamically unstable within 1 Myr after the CE phase. ^cGiven the large parameter space and the associated uncertainties, we consider e_CBP > 0.95 as ejection. ^dP1/P2/P3 have 60%/15%/85% probability to become dynamically unstable within 1 Myr after the CE phase.

Download table as:  ASCIITypeset image

Table 8.  Same as Table 7

αCE	Tides		aCBP,PCE (au)			eCBP,PCE		Notes
		Mode	Median	Range	Mode	Median	Range
Kepler-64								[aCBP,0 = 0.65]
0.5	NTCP	${0.8}_{-0.2}^{+0.6}$	0.9	0.6–4.8	${0.16}_{-0.}^{+0.24}$	0.35	0.03–0.85	⋯
0.5a	TCP	${0.6}_{-0.0}^{+0.6}$	0.8	0.6–1.2	${0.19}_{-0.07}^{+0.23}$	0.19	0.11–0.43	⋯
1.0	NTCP	${1.4}_{-0.0}^{+1.1}$	2.0	1.3–7.5	${0.62}_{-0.31}^{+0.15}$	0.53	0.12–0.8	⋯
1.0a	TCP	${0.9}_{-0.0}^{+0.1}$	0.9	0.9–1.0	${0.2}_{-0.18}^{+0.0}$	0.13	0.01–0.21	⋯
3/5/10	NTCP	${1.6}_{-0.1}^{+2.6}$	3.1	1.4–19.6	${0.8}_{-0.26}^{+0.11}$	0.66	0.04–0.95b	4% ejected
3/5/10	TCP	${2.4}_{-0.1}^{+0.}$	2.4	2.2–2.6	${0.49}_{-0.09}^{+0.01}$	0.44	0.38–0.5	⋯
Kepler-1647c								[aCBP,0 = 2.71]
0.5	NTCP	${1.6}_{-0.0}^{+2.3}$	2.4	1.4–39.5	${0.72}_{-0.43}^{+0.06}$	0.69	0.3–0.95	47% ejected
0.5a	TCP	${2.3}_{-0.0}^{+10.6}$	3.6	2.0–13.2	${0.26}_{-0.01}^{+0.5}$	0.33	0.23–0.79	⋯
1.0	NTCP	${1.9}_{-0.0}^{+2.1}$	2.6	1.7–37.5	${0.53}_{-0.11}^{+0.11}$	0.57	0.4–0.95	47% ejected
1.0a	TCP	${3.9}_{-0.0}^{+5.6}$	5.5	3.8–9.7	${0.29}_{-0.0}^{+0.4}$	0.5	0.27–0.71	⋯
3/5/10	NTCP	${2.3}_{-0.3}^{+1.5}$	3.1	1.9–47.1	${0.4}_{-0.12}^{+0.12}$	0.44	0.16–0.95	54% ejected
3/5/10	TCP	${9.6}_{-1.6}^{+4.6}$	11.5	7.6–23.2	${0.7}_{-0.01}^{+0.14}$	0.76	0.63–0.88	⋯
Kepler-1647d
5.0a	NTCP	${11.3}_{-1.1}^{+5.3}$	15.3	10.2–53.3	${0.76}_{-0.03}^{+0.03}$	0.79	0.68–0.95	86% ejected
10.0a	NTCP	${8.2}_{-0.6}^{+6.4}$	10.5	7.6–36.9	${0.64}_{-0.07}^{+0.07}$	0.71	0.57–0.95	83% ejected
10.0	TCP	${14.5}_{-1.0}^{+11.1}$	18.2	13.5–96.7	${0.52}_{-0.4}^{+0.28}$	0.5	0.01–0.95	85%ejected
Kepler-35e								[aCBP,0 = 0.60]
0.5	NTCP	${0.87}_{-0.01}^{+0.23}$	1.04	0.75–10.13	${0.44}_{-0.18}^{+0.27}$	0.41	0.03–0.93	4% ejected
0.5a	TCP	${0.87}_{-0.01}^{+0.47}$	1.01	0.84–1.37	${0.43}_{-0.37}^{+0.01}$	0.25	0.04–0.45	⋯
1/3/5/10	NTCP	${0.84}_{-0.01}^{+0.16}$	1.03	0.77–6.93	${0.43}_{-0.26}^{+0.01}$	0.36	0.04–0.91	⋯
1/3/5/10	TCP	${1.0}_{-0.04}^{+0.03}$	0.98	0.86–1.24	${0.18}_{-0.04}^{+0.04}$	0.19	0.01–0.39	⋯

Notes.

^aSinusoidal-like distributions in a_CBP,PCE/e_CBP,PCE as a function of Δθ₀; typically with double-peaked PDFs. ^bGiven the large parameter space and the associated uncertainties, we consider $e\gt 0.95$ as ejection. ^cPrimary RLOF and CE. ^dSecondary RLOF and CE. ^eExperiences CE for Z − 1σ_Z.

Download table as:  ASCIITypeset image

Table 10.  Same as Table 9

αCE	Tides	aCBP,PCE (au)			eCBP,PCE			Notes
		Mode	Median	Range	Mode	Median	Range
Kepler-64								[aCBP,0 = 0.65]
0.5	NTCP	0.82	0.82	⋯	0.25	0.25	0.22–0.27	⋯
0.5${}^{{\rm{f}}}$	TCP	0.72	0.72	⋯	0.02	0.02	⋯	⋯
1.0	NTCP	1.58	1.58	⋯	0.31	0.33	0.31–0.35	⋯
1.0${}^{{\rm{f}}}$	TCP	0.91	0.91	⋯	0.01	0.01	⋯	⋯
3/5/10	NTCP	1.83	1.83	⋯	0.3	0.35	0.31–0.39	⋯
3/5/10	TCP	1.9	1.9	⋯	0.02	0.02	⋯	⋯
Kepler-1647b								[aCBP,0 = 2.71]
0.5	NTCP	${2.52}_{-0.01}^{+1.53}$	3.55	2.33–40.61	${0.7}_{-0.08}^{+0.02}$	0.71	0.62–0.95a	27% ejected
0.5${}^{{\rm{f}}}$	TCP	${3.1}_{-0.01}^{+0.19}$	3.19	3.08–3.3	${0.11}_{-0.08}^{+0.01}$	0.07	0.02–0.12	⋯
1.0	NTCP	${3.91}_{-0.01}^{+0.67}$	4.45	3.83–20.27	${0.72}_{-0.01}^{+0.22}$	0.81	0.71–0.95	27% ejected
1.0${}^{{\rm{f}}}$	TCP	4.1	4.1	⋯	0.07	0.07	⋯	⋯
3/5/10	NTCP	${4.7}_{-0.26}^{+0.63}$	4.93	4.29–19.14	${0.7}_{-0.01}^{+0.19}$	0.76	0.57–0.95	4% ejected
3/5/10	TCP	4.93	4.93	⋯	0.08	0.08	⋯	⋯
Kepler-1647c
5.0e	NTCP	${7.9}_{-0.7}^{+6.3}$	9.9	7.2–55.9	${0.48}_{-0.1}^{+0.9}$	0.55	0.2–0.95	82% ejected
10.0e	NTCP	${7.7}_{-0.8}^{+4.0}$	9.2	6.9–31.1	${0.32}_{-0.08}^{+0.25}$	0.54	0.24–0.95	81% ejected
10.0	TCP	⋯	⋯	⋯	⋯	⋯	⋯	100% ejected
Kepler-35d								[aCBP,0 = 0.60]
0.5	NTCP	0.94	0.94	⋯	0.37	0.37	0.33–0.41	⋯
0.5	TCP	0.95	0.95	⋯	0.03	0.03	⋯	⋯
1/3/5/10	NTCP	0.94	0.94	⋯	0.27	0.27	0.18–0.39	⋯
1/3/5/10	TCP	0.95	0.95	⋯	0.02	0.02	⋯	⋯

Notes.

^aGiven the large parameter space and the associated uncertainties, we consider e > 0.95 as ejection. ^bPrimary RLOF and CE. ^cSecondary RLOF and CE. ^dExperiences CE for Z − 1σ_Z. ^eNonnormal PDF with multiple peaks of comparable strength, with mode not well defined.

Download table as:  ASCIITypeset image

Overall, five of the nine Kepler CBP systems experience at least one CE phase for their default metallicities (see Table 1) over the duration of the BSE simulations (15 Gyr): Kepler-34, -38, -47, -64, and -1647. Kepler-35 experiences a CE for Z − 1σ_Z, and Kepler-453 for Z − 2σ_Z.

4.1. Kepler-34

4.2. Kepler-38

The binary consists of 0.95 ${M}_{\odot }$ and 0.25 ${M}_{\odot }$ stars on a slightly eccentric (e_bin,0 = 0.1), 19 day orbit that changes little in the first ∼13 Gyr. The CBP has an initial orbit of 0.46 au. The main evolution stages of the binary star and their effect on the CBP are as follows.

4.2.1. The Binary

The primary star fills its Roche lobe at t ∼ 13.7 Gyr, and the binary enters a primary CE stage. During this stage, Ψ(${\alpha }_{1}$) ≈ 0.04 and α_crit  ≈ 26 ${M}_{\odot }$ yr⁻¹, much larger than ${\alpha }_{1}$ and ${\alpha }_{0.1}$. As a result of the CE, the binary merges as a first giant branch star and loses ∼25–40% of its mass for α_CE = 0.5; here β > β_eject (see Equation (3)), and the planet should remain bound in a runaway regime (i.e., Ψ ≫ 1) even at pericenter (see Equations (37)–(44), Veras et al. 2011). The system continues to slowly lose mass by the end of the BSE simulations (15 Gyr), so a_CBP,PCE adiabatically expands by up to a factor of 2.

For α_CE = 1.0, 3.0, 5.0, and 10.0, the binary orbit shrinks by a factor of 4–30 and its mass decreases by ∼60%; here β < β_eject, that is, the planet should be ejected in a runaway regime. For α_CE  = 1 and TCP, the system experiences a secondary RLOF at t ∼ 13.74 Gyr and coalesces into a first giant branch star without mass loss (thus with no changes in a_CBP,PCE or e_CBP,PCE). The binary does not lose additional mass by the end of the BSE simulations.

Overall, by 15 Gyr the binary evolves into either a merged helium or carbon/oxygen white dwarf (HeWD, COWD) or a very close HeWD–MS star binary.

The orbital reconfiguration of the binary system during the primary CE stage is shown on the left panels of Figure 2 for α_CE = 0.5 (upper panels) and α_CE = 1.0, 3.0, 5.0, and 10.0 (lower panels). The observer is looking along the y axis, from below the figure. The maximum (red dashed line; binary initially at EE, CBP initially at phase 0.73) and minimum (red solid line; binary at EE, CBP at phase 0.08) orbital expansion of the CBP is shown on the right panels of the figure for ${\alpha }_{1}$ (red), as well as the maximum expansion for ${\alpha }_{0.1}$ (green solid line). To test the dynamical stability of the CBP at maximum expansion (dashed red line), we integrated the system for a thousand planetary orbits. The orbit precesses but does not exhibit chaotic behavior for the duration of the integration. We leave examination of the long-term behavior and stability of the system for future work.

Zoom In Zoom Out Reset image size

4.2.2. The CBP

In Figure 3 we show the α_CE  = 0.5 CE evolution of M₁, of a_bin (blue lines, upper panel), and of a_CBP  (middle panels) and e_CBP (lower panel) for ${\alpha }_{1}$ (red) and ${\alpha }_{0.1}$ (green); on Figure 4 we show the corresponding evolution for α_CE  = 1, 3, 5, and 10. Unlike the case of a planet orbiting a single star, where the evolution of the planet's orbit follows a single path for each ${\alpha }_{M}$ (see Figure 1), the dynamical response of a CBP to a binary undergoing a CE stage depends on the initial configuration of the system, that is, the phase difference Δθ₀ ≡ θ_CBP,0 − θ_bin,0 − ω_bin, where θ_CBP,0and θ_bin,0 are the initial true anomaly of the CBP and the binary, respectively, and ω_bin is the binary's argument of pericenter. This is reminiscent of the importance of the initial true anomaly of a planet orbiting a single star on an eccentric orbit (see Veras et al. 2011). The dependence is shown in the second panel of Figure 3, where each curve represents the evolution of a_CBP  for Δθ₀ varying between 0 and 1 (for clarity, the curves here represent every tenth phase where in the simulations the phase difference varies with a step of 0.02).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Since our simulations encompass a discrete set of initial conditions for α_CE, tides (NTCP or TCP for minimum to maximum strength, respectively), and Δθ₀, the true evolution of a_CBP  and e_CBP would be continuous. We represent this by the hatched regions on the lower panels of Figures 3 and 4, where the different hatches correspond to TCP or NTCP. For clarity, we show individual a_CBP  evolutionary tracks only in Figure 3; these are replaced with hatched regions in all subsequent figures. We note that the distribution of a_CBP,PCE and e_CBP,PCE  (i.e., a_CBP  and e_CBP at the end of the CE phase) is neither uniform nor necessarily normal, as seen from the second panel on Figure 3, but instead depends on Δθ₀. This dependence is demonstrated on Figure 5 for α_CE  = 0.5 and on Figure 6 for α_CE  = 1, 3, 5, and 10, and the corresponding probability distribution functions (PDFs) are shown in Figures 7 and 8. The PDFs can have a prominent peak with a long one-sided tail (e.g., upper left panel on the first figure), or be double-peaked with the peaks near the edges (e.g., lower left panel on the first figure); this is a natural consequence of the vaguely sinusoidal dependencies seen in Figures 5 and 6. Given the prominent diversity of the PDFs, the distributions are better represented by their modes than their medians. Thus for completeness we list the mode, median, and min/max range for each a_CBP,PCE and e_CBP,PCE  in Tables 7 through 10, but we cite the 68%-range upper and lower bounds on the modes only because the medians may not be appropriate in some cases. Unless otherwise noted, we refer to the modes of a_CBP,PCE and e_CBP,PCE  as the default results.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It is important to note that the osculating orbital elements of a CBP do not describe a closed ellipse during the CE phase, and its orbital motion is not Keplerian. Instead, the orbit spirals outward with time and, as seen from Figures 3 and 4, while the binary is losing mass, a_CBP  and e_CBP oscillate together with the oscillations in e_bin (see inset figure, lower panel). Similar behavior was observed in numerical simulations of planets around single stars experiencing mass loss (Veras et al. 2011; Adams et al. 2013), where the planet's orbital elements oscillate on the planet's period due to the different effects of mass loss at pericenter and apocenter. This is reproduced in our simulations for the binary star itself; the orbital elements of the secondary star oscillate on the binary period as seen from the inset figure in the lower panel of Figure 3. The CBP, initially on a circular orbit, quickly gains eccentricity and responds to the binary's oscillations.

Overall, the orbital expansion of the Kepler-38 CBP for ${\alpha }_{0.1}$ is fully consistent with the adiabatic approximation, and the planet gains slight eccentricity (see Table 9). The outcome for the case of ${\alpha }_{1}$ is diverse: a_CBP,PCE ranges from below the corresponding adiabatic approximation up to 3.8 au, and e_CBP,PCE  can reach up to 0.8 (see Table 7).

4.3. Kepler-47

The binary consists of 0.96 ${M}_{\odot }$ and 0.34 ${M}_{\odot }$ stars on a nearly circular (e_bin,0 = 0.02), 7.45 day orbit that changes little in the first ∼12 Gyr. The three CBPs have initial orbits of 0.3, 0.72, and 0.99 au and masses of 2 M_⊕, 19 M_⊕, and 3 M_⊕ (J. Orosz 2016, private communication), from inner to outer, respectively. For completeness, we integrate the system for the appropriate CBP masses to take into account planet–planet interactions. The main evolutionary stages are as follows.

4.3.1. The Binary

The primary star experiences an RLOF at t ∼ 12 Gyr, and the binary enters a CE stage. During this stage, Ψ(${\alpha }_{1}$) ≈ 0.02/0.06/0.1 and α_crit  ≈ 60/16/10 ${M}_{\odot }$ yr⁻¹ for planets 1, 2, and 3, respectively. For ${\alpha }_{1}$ the time of CE mass loss (${T}_{\mathrm{CE}}$ ) is comparable to the orbital periods of planets 2 and 3, that is, T_CE/P_CBP ∼ 0.5–2, indicating that these two CBPs evolve in the transition regime.

As a result of the primary CE, for α_CE  = 0.5/1 the binary merges as a first giant branch star and loses ∼15–40% of its mass; β < β_eject, and the CBPs should remain bound even in a runaway regime. From the end of the CE to the end of the BSE simulations, the system slowly loses 50% of its mass, and a_CBP,PCE of those planets that remain long-term stable after the CE (discussed below) expand adiabatically by a factor of 2.

For α_CE  = 3/5/10, the binary evolves into a HeWD–MS star pair and loses 56% of its mass, and a_bin decreases by a factor of 4–10; here β > β_eject, and the CBPs should be ejected in a runaway regime. The system does not experience further mass loss by the end of the BSE simulations.

4.3.2. The CBPs

In Figures 9 and 10 we show the evolution of a_CBP  and e_CBP for the three CBPs (and for ${\alpha }_{1}$ and ${\alpha }_{0.1}$) caused by the primary RLOF and CE. For simplicity here and for the rest of the Kepler  systems, we do not show the evolution of M₁ and a_bin since it is qualitatively very similar to the case of Kepler-38. We note that, unlike the corresponding figures for Kepler-38, here the colors represent the three planets (magenta, green, and red for planets 1, 2, and 3, respectively), and the different panels represent the evolution of a_CBP  and e_CBP for ${\alpha }_{1}$ (upper and middle panels) and for ${\alpha }_{0.1}$ (lower panels).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As seen from the lower panels in Figures 9 and 10, the orbital expansion of all three of Kepler-47's CBPs is fully consistent with the adiabatic approximation for ${\alpha }_{0.1}$, and the planets gain slight eccentricity (not shown in the figures, but listed in Table 9). The orbital evolution of planets 2 and 3 for the primary RLOF is much richer for ${\alpha }_{1}$, as shown on the upper and middle panels of Figures 9 and 10 (see also Table 7). The planets gain high eccentricities and for α_CE  = 3/5/10 may even become ejected during the CE stage, where, given the complexity, uncertainties, and large parameter space of the evolution of the binary and of the CBP, here and throughout the paper we define ejection as e_CBP > 0.95. Specifically, planet 3 becomes unbound in 46% of the NTCP simulations and in 81% of the TCP simulations.

We note that, during the primary CE stage, for the case of ${\alpha }_{1}$, the orbit of the inner CBP grows adiabatically, while those of the outer two do not (see upper panel of Figure 10). This is because the orbital period of the inner planet is much shorter than the duration of ${T}_{\mathrm{CE}}$, while those of the outer two are comparable. Thus a multiplanet CB system can experience two regimes of planetary orbital evolution during the same CE stage.

The dependence of a_CBP,PCE and e_CBP,PCE of planet 2 on Δθ₀ is shown in Figures 11 and 12. The planet's orbit evolves adiabatically for ${\alpha }_{0.1}$ (upper right panels), nonadiabatically for ${\alpha }_{1}$, and remains bound during the CE in both cases, regardless of Δθ₀.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Even if an ejection does not occur during the primary CE phase, a_CBP,PCE and e_CBP,PCE of all three CBPs for ${\alpha }_{1}$ are such that the planets will continue to interact with each other after the end of the CE in such a way that further ejections or collisions are possible. To study the long-term evolution of the system, we integrated the equations of motion for all five bodies (two stars and three planets) for 1 Myr, using the planets' mode(a_CBP,PCE) and mode(e_CBP,PCE), and randomizing their initial arguments of pericenter and true anomalies. Any planet that achieves e_CBP > 0.95 is marked as ejected and removed from the simulations, and collision is defined as a planet–planet separation smaller than d_min = R_P,i + R_P,j.

Overall, planet 2 dominates the 1 Myr dynamical evolution and has the highest probability of remaining bound to the system. For α_CE  = 0.5/1 the modes of a_CBP,1/2/3 have 40%/15%/75% ejection probability in 1 Myr. Where planet 2 remains bound, its orbital elements remain mostly within the respective 68% range of mode(a_CBP,PCE,2) and mode(e_CBP,PCE,2).¹⁷ Here the modes of a_CBP,1/2/3 have 60%/15%/85% ejection probability in 1 Myr, and the orbital elements of planet 2 again remain mostly within the 68% range of mode(a_CBP,PCE,2) and mode(e_CBP,PCE,2).

4.4. Kepler-64

The binary consists of 1.53 ${M}_{\odot }$ and 0.41 ${M}_{\odot }$ stars on a slightly eccentric (e_bin,0 = 0.22), 20 day orbit that changes little in the first ∼3 Gyr; the CBP has an initial orbit of 0.65 au. The system evolves as follows.

4.4.1. The Binary

The primary star fills its Roche lobe at t ∼ 3.02 Gyr, and the binary enters a primary CE stage. During this stage, Ψ(${\alpha }_{1}$) ≈ 0.03 and α_crit  ≈ 32 ${M}_{\odot }$ yr⁻¹, much larger than ${\alpha }_{1}$ and ${\alpha }_{0.1}$. However, for ${\alpha }_{1}$ the time of mass loss (T_CE) is comparable to the orbital period of the CBP, that is, T_CE/P_CBP ∼ 1.1–3.4, depending on α_CE. Thus while the CBP is theoretically in the adiabatic regime, the comparable timescales result in a nonadiabatic dynamical response.

The binary merges as a first giant branch star and loses ∼20–58% of its mass for α_CE  = 0.5/1. In this regime, β > β_eject for α_CE  = 0.5, and for α_CE  = 1, TCP and the planet should remain bound in a runaway regime even at pericenter; β < β_eject for α_CE  = 1.0 and NTCP, indicating a potential ejection of the CBP in a runaway regime. For α_CE  = 0.5 and 1, the system continues to slowly lose mass by the end of the BSE simulations (15 Gyr), so a_CBP,PCE expands adiabatically by up to a factor of 3.

For α_CE = 3.0, 5.0, and 10.0 the binary shrinks by a factor of 5–25 and its mass decreases by ∼65%; here, β < β_eject, and the CBP should be ejected in a runaway regime. For α_CE  = 3 and TCP, the system experiences a secondary RLOF at t ∼ 6.2 Gyr and coalesces without mass loss (no changes in a_CBP  or e_CBP). Except for α_CE  = 3 and TCP, the binary does not lose mass after the primary CE.

The binary reaches 15 Gyr as either a merged COWD or a very close HeWD–MS star binary.

4.4.2. The CBP

In Figure 13 we show the evolution of a_CBP  and e_CBP for ${\alpha }_{1}$ (red) and ${\alpha }_{0.1}$ (green), respectively, for α_CE  = 0.5/1 (upper two panels), which result in a binary merger, and for α_CE  = 3/5/10 (lower two panels). The dependence of a_CBP,PCE and e_CBP,PCE  on the initial phase difference between the binary and the CBP (Δθ₀) is shown in Figure 14 for α_CE  = 0.5/1.0 and in Figure 15 for α_CE  = 3/5/10.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Similar to Kepler-38, the orbital expansion of Kepler-64's CBP is fully consistent with the adiabatic approximation for ${\alpha }_{0.1}$, and again the planet gains slight eccentricity (Table 10). The outcome for the case of ${\alpha }_{1}$ is quite diverse. Depending on Δθ₀, the planet can gain very high eccentricity and even become ejected (e_CBP > 0.95). This is shown as missing points in Figures 14 and 15. The Kepler-64 CBP can reach such eccentricities in ∼4% of our simulations for ${\alpha }_{1}$ and α_CE  = 3/5/10. Interestingly, for α_CE  = 0.5 and TCP, a_CBP  can decrease during the CE, and mode(a_CBP,PCE) = ${0.6}_{-0.0}^{+0.6}\,[\mathrm{au}]$ is smaller than the initial orbit a_CBP,0 = 0.65 au, though the mode's 68% range is significant. On the other side of the spectrum, a_CBP  can reach up to ∼20 au for α_CE  = 3/5/10 and NTCP.

4.5. Kepler-1647

The initial binary consists of 1.22 ${M}_{\odot }$ and 0.97 ${M}_{\odot }$ stars on a slightly eccentric (e_bin,0 = 0.15), 11 day orbit that changes little in the first ∼5 Gyr. The CBP has an initial orbit of 2.7 au. The main stages of the evolution of the system are as follows.

4.5.1. The Binary

The primary star fills its Roche lobe at t ∼ 5.4 Gyr, and the binary enters a primary CE stage. During this stage, Ψ(${\alpha }_{1}$) ≈ 0.2 and α_crit  ≈ 4.6 ${M}_{\odot }$ yr⁻¹, comparable to ${\alpha }_{1}$. Additionally, for NTCP and ${\alpha }_{0.1}$, the time of mass loss (T_CE) is comparable to the orbital period of the CBP, that is, T_CE/P_CBP ∼ 1.3/3.0 for α_CE  = 0.5/1.0, respectively. Thus the orbit of the CBP evolves in the transition regime.

Compared to the other Kepler  CBP systems, the binary star of Kepler-1647 experiences the richest evolution. For α_CE  = 0.5/1 the binary merges as a first giant branch star and loses ∼15–40% of its mass. In this regime, β > β_eject, and the planet should remain bound in a runaway regime even at pericenter. The system continues to slowly lose mass by the end of the BSE simulations (15 Gyr), so a_CBP,PCE further expands (adiabatically) by up to a factor of 3.

For α_CE = 3/5/10, the binary shrinks by a factor of 2–7 into a HeWD–MS star binary, and its mass decreases by ∼45%, close to the critical mass loss for runaway ejection (0.5). Here, β < β_eject and the CBP should remain bound even in a runaway regime.

The binary experiences a secondary RLOF and CE for α_CE  = 3/5/10 (see Table 6 for details), and its subsequent evolution can follow several paths. By the end of the BSE simulations, the binary (1) merges without mass loss (i.e., no changes in a_CBP,PCE and e_CBP,PCE ) into a first giant branch for α_CE  = 3/5 and TCP, which by 15 Gyr evolves into a COWD; and (2) evolves into a PCE WD–WD binary (with mass loss, thus a_CBP,PCE and e_CBP,PCE  change) for α_CE  = 5 (and for NTCP) and for α_CE  = 10 (for both TCP and NTCP), which can experience a third and final RLOF, triggering an SN explosion (for α_CE  = 5/10 and TCP).

4.5.2. The CBP

The corresponding evolution of a_CBP  and e_CBP for the primary RLOF and for ${\alpha }_{1}$ and ${\alpha }_{0.1}$ is shown in Figure 16. Unlike the previous systems, the ${\alpha }_{0.1}$ case for Kepler-1647 represents an adiabatic evolution only for the TCP. For ${\alpha }_{0.1}$ and NTCP, where ${T}_{\mathrm{CE}}$  and P_CBP are comparable for all α_CE, the CBP's orbit evolves in the transition regime (as it does for ${\alpha }_{1}$).

Zoom In Zoom Out Reset image size

Specifically, for ${\alpha }_{0.1}$ and NTCP, the CBP reaches sufficiently high eccentricities (e_CBP > 0.95) in ∼5–25% of the simulations to be ejected from the system; for ${\alpha }_{1}$ the planet is ejected in ∼45–55% of the simulations. In addition, for α_CE  = 0.5/1 and NTCP, and for α_CE  = 3/5/10 and NTCP, a_CBP  can decrease during the primary CE phase, and the corresponding mode of a_CBP,PCE (1.6–2.3 au, depending on α_CE, with a wide 68% range) is smaller than the initial semimajor axis of the CBP (2.71 au). Some of the simulations produce the opposite result, namely orbital expansion by up to a factor of ∼20 such that a_CBP,PCE reaches ∼50 au. The corresponding dependence of a_CBP,PCE and e_CBP,PCE  on Δθ₀ is shown in Figures 17 and 18, where the missing Δθ₀ phase coverage in all panels represents planet ejection.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Using the respective a_bin,PCE and e_bin,PCE and the modes of a_CBP,PCE and e_CBP,PCE  from the end of the preceding primary RLOF as initial conditions, we further explore the evolution of the CBP for the secondary RLOF and CE stage for the three cases where the binary evolves into a WD–WD pair with mass loss: (1) α_CE  = 5, NTCP; (2) α_CE  = 10, TCP; and (3) α_CE  = 10, NTCP. The results are as follows.

During this secondary CE phase, for (1) and (2) Ψ(${\alpha }_{1}$) ≈ 0.4 and α_crit  ≈ 2.4 ${M}_{\odot }$ yr⁻¹, comparable to ${\alpha }_{1}$, indicating that the CBP evolves in the transition regime. For (3) Ψ(${\alpha }_{1}$) ≈ 3.6, α_crit  ≈ 0.3 ${M}_{\odot }$ yr⁻¹, and the evolution of the planet's orbit is in the runaway regime. The time of mass loss (T_CE) is close to the orbital period of the CBP, that is, T_CE/P_CBP ∼ 0.3–2.5. The binary shrinks by a factor of ∼5–10 into a WD–WD pair and loses ∼60–70% of its mass. As a result, β < β_eject and the CBP should become unbound in a runaway regime.

The evolution of a_CBP  and e_CBP during the secondary CE phase is shown in Figure 19, where the upper panel represents case (1) (α_CE  = 5, NTCP), the middle panel case (2) (α_CE  = 10, NTCP), and the lower panel case (3) (α_CE  = 10, TCP). The CBP is ejected in more than ∼80% of the simulations. For ${\alpha }_{1}$, the CBP's a_CBP  expands to a mode value of ∼8–15 au (depending on the scenario described above), with a maximum of up to 50–100 au; the mode of e_CBP is in the range of 0.5–0.75, reaching a maximum of 0.95. For ${\alpha }_{0.1}$, the CBP's orbit can expand up to ∼30–50 au.

Zoom In Zoom Out Reset image size

By the end of the BSE simulations, the CBP can remain bound only for case (2) as the other two scenarios trigger an SN shortly after the secondary CE.

4.6. Kepler-35

Several key results emerge from our study:

(1) Five of the nine Kepler  CBP hosts experience at least one RLOF and a CE phase for their default parameters (masses, metallicities) by the end of our BSE simulations (15 Gyr); Kepler-35 experiences a primary RLOF for Z − 1σ_Z. Depending on the treatment of the CE stage, the binaries either merge after a primary or secondary CE (typically for α_CE  = 0.5/1) or evolve into very short-period WD–MS pairs (for α_CE  = 1/3/5/10)¹⁸ ; Kepler-1647 evolves into a WD–WD pair after a secondary RLOF and CE for α_CE  = 5/10. Two systems trigger a double-degenerate supernova explosion: Kepler-34 for primary RLOF, and Kepler-1647 for a third RLOF and α_CE  = 5, TCP, and α_CE  = 10, NTCP. The binary systems lose a tremendous amount of mass—up to ∼60–70%—during the CE stages and can shrink to sub-R_⊙ orbital separations.

(2) From the dynamical perspective presented here, Kepler CBPs predominantly remain bound to their host binaries after the respective CE phases even for mass-loss rates as high as 1.0 ${M}_{\odot }$ yr⁻¹. There are only four scenarios where the respective CBP has a 100% probability of becoming unbound: (a) for Kepler-34 where the SN explosion occurs with the first CE; (b) for the secondary CE phase of Kepler-1647 for ${\alpha }_{0.1}$, α_CE  = 10, TCP; (c) for the SN caused by the third RLOF of Kepler-1647, α_CE  = 5, NTCP; and (d) SN caused by the third RLOF of Kepler-1647, α_CE  = 10, TCP. In all other scenarios (106 total, not accounting for the initial phase differences Δθ₀), the CBPs remain bound in the majority of cases (only Kepler-1647 has nonnegligible ejection probability, again highly dependent on Δθ₀).

For mass-loss rates of 0.1 ${M}_{\odot }$ yr⁻¹, the orbits of the CBPs evolve adiabatically—well reproduced by REBOUNDx—except for Kepler-1647. The orbital expansion of some of the Kepler  CBPs is consistent with the adiabatic approximation even for mass-loss rates of 1.0 ${M}_{\odot }$ yr⁻¹ (e.g., inner planet of Kepler-47, Kepler-64); in other cases, the mode of a_CBP,PCE is smaller than the adiabatic approximation (e.g., Kepler-1647).

According to the analytical prescription, some systems should retain their CBPs even in a runaway mass-loss regime (Ψ ≫ 1). Our code fully reproduces the runaway approximation for single-star systems.

(3) The transition regime (Ψ ∼ 0.1–1)—where the evolution of some of Kepler  CBPs falls for ${\alpha }_{1}$—is complex, and the final eccentricities and semimajor axes of the planets depend both on the treatment of the CE stage (i.e., α_CE, treatment of tides) and on the initial configuration of the system (i.e., Δθ₀). We find that the orbits of Kepler  CBPs can expand by more than an order of magnitude over the course of a single CBP year, reaching a_CBP,PCE of tens of astronomical units (during the secondary CE of Kepler-1647, a_CBP  can expand from ∼10 au to ∼100 au; see Table 8).

Multiplanet CBP systems add yet another level of complexity because they can experience both regimes simultaneously. For example, where the orbit of the inner Kepler-47 CBP expands adiabatically, the middle and outer CBP orbits grow nonadiabatically—all during the primary CE phase.

Overall, the CE-induced orbital evolution of CBPs is a dynamically rich and complex process in which a planet can migrate adiabatically during one CE phase and nonadiabatically during another. As a result, a CBP cannot experience the same CE twice.

(4) If CE mass-loss rates are indeed high (e.g., ∼1.0 ${M}_{\odot }$ yr⁻¹)—as suggested by theoretical work—and the orbits of Kepler-like CBPs evolve (and survive) according to our simulations, we should expect to detect potentially highly eccentric planets orbiting PCE systems on very large orbits. Alternatively, if Ψ ≪ 1, then we should find PCE CBPs on low-eccentricity orbits. Interestingly, recent observational efforts have produced a number of PCE eclipse-time variations (ETV) CBP candidates on wide, notably eccentric orbits (e.g., Zorotovic & Schreiber 2013 and references therein), suggesting a possible connection with PCE Kepler-like CBPs.

In Figure 20 we show a_CBP,PCE versus e_CBP,PCE of Kepler's CBPs¹⁹ and compare them to those of the PCE ETV CBP candidates.²⁰ We caution that the comparison is not direct as the CE stages of the different Kepler systems do not occur all at the same time, and the PCE ETV CBPs do not have identical ages. Instead, given the handful of known targets—Kepler's nine compared to estimated millions of similar CBPs (Welsh et al. 2012) and a dozen PCE ETV CBPs—the lower two panels on the figure portray a potential distribution of an underlying PCE CBP population with a range of ages and evolutionary stages.

Zoom In Zoom Out Reset image size

Overall, the PCE orbital configurations of Kepler's CBPs are qualitatively consistent with those of the observed population of the currently known PCE ETV CBP candidates.²¹ Thus our results assist in both interpreting the nature of these candidates and in guiding future observational efforts to discover new systems. Such discoveries will provide a deeper understanding of the evolution of planets in a binary star system and also much-needed observational constraints on the stellar astrophysics of the complex CE phase, as "CE is one of the most important unsolved problems in stellar evolution," according to Ivanova et al. (2013), "and is arguably the most significant and least-well-constrained major process in binary evolution" (but also see Taam & Sandquist 2000; Taam & Ricker 2010 and Webbink 2008 for alternative reviews).

(5) A CB body gravitationally perturbs its host binary star, and if the latter is eclipsing, these perturbations can be manifested as deviations from linear ephemeris in the measured stellar eclipse times. The cause of these variations can be either dynamical, where the third body's influence changes the orbital elements of the binary star, or a light travel time effect, where the tertiary object and the binary revolve around a common center of mass. ETVs are indeed a powerful and highly productive method of discovering and studying stellar triple and higher order systems (e.g., Borkovits et al. 2016; Orosz 2015, and references therein), and they have recently been used to study CBPs as well. In particular, the former effect has been measured for several of the Kepler CBP systems and used to constrain the respective planetary masses (Doyle et al. 2011; Welsh et al. 2012; Kostov et al. 2016). The latter effect has been suggested as the cause for measured ETVs for a number of PCE EB systems with proposed CB companions (as mentioned above; also see Volschow et al. (2016) and references therein).

Assuming that the orbits of the Kepler  CBPs studied here evolve only dynamically and their masses remain constant (i.e., no changes due to, for example, interactions with a CE-triggered CB disk), it is informative to evaluate the respective planets' ETV-based detectability after the CE stages of the relevant systems. To calculate the expected amplitudes of the two effects discussed above, A_dyn and A_LTTE, we use the formalism of Borkovits et al. (2012, 2015).

The respective amplitudes and periods for Kepler-35 (M_p = 0.13 M_Jup) and Kepler-1647 (M_p = 1.5 M_Jup)—the only two systems with well-constrained CBP masses—are listed in Table 11 (for the respective range of binary and CBP parameters listed in Tables 8) and compared to the detected signals for the two PCE CBP candidates NN Ser and HW Vir. As seen from Table 11, detection of such ETVs is feasible in terms of both the amplitudes and periods of the expected signal, which are qualitatively similar to the ETV signals of the two PCE CBP candidates.

Table 11.  Expected Dynamical and Light Travel Time ETV Amplitudes and Periods after the Respective CE Stages for Kepler CBPs with Known Masses (Kepler-1647 and Kepler-35)

αM	Adyn (min/max)	ALTTE(min/max)	P(min/max)	Notes
(${M}_{\odot }$ yr−1)	(s)	(s)	(year)
Kepler-1647ba
1.0	0.2/13.4b	68/1682	2.4/292	⋯
0.1	0.1/13	154/680	8.1/75	⋯
Kepler-1647bc
1.0	0./0.2	736/9680	31/1490	⋯
0.1	0./0.2	668/5420	27/624	⋯
Kepler-35ba
1.0	0.1/56.5	4.3/36.9	0.7/17	⋯
0.1	0.2/44.2	5.0/5.1	0.9/0.9	⋯
NN Ser bc
	∼5/30d	∼5/30	∼7/15	⋯
HW Vir bc
	∼55/550d	∼50/550	∼13/55	⋯

Notes. For comparison, we also list the detected ETV amplitudes and periods for the CBP candidates in NN Ser and HW Vir.

^aAfter primary CE. ^bMin/max for the range of parameters from Table 8. ^cAfter secondary CE. ^dShowing the combined effect of inner/outer CBPs (instead of min/max).

Download table as:  ASCIITypeset image

Another option for the detection of a PCE CBP is to directly image the planet. The benefits for such detection are twofold: (a) the contrast ratio between the binary and the CBP decreases after each CE phase as the primary or secondary evolves from a luminous MS star to a much fainter WD; and (b) the angular separation between the binary and the CBP increases²² (e.g., Hardy et al. 2015). There is potentially a third benefit, where the CBP can accrete mass after the CE stage (e.g., Zorotovic & Schreiber 2013; Bear & Soker 2014) and becomes brighter, thus further decreasing the contrast ratio.

Recently, Hardy et al. (2015) used VLT/SPHERE to observe the V471 Tau system—an eclipsing binary star with measured ETVs—with the goal to directly image a substellar-mass CB candidate suspected to be the cause of the ETVs. Given the known age of the system and the masses of the binary and the CB candidate, direct detection of such a tertiary body should have been well within the capabilities of the instrument. Hardy et al. (2015), however, report a null detection, casting doubt on the interpretation of the measured ETV signal (but see Vaccaro et al. 2015 for an alternative explanation). Nevertheless, based on our results for the dynamical survivability of CBPs around evolving binary stars, and on estimates of the occurrence rates of planets orbiting short-period MS binary stars (e.g., Welsh et al. 2012), we encourage the continuation of such direct imaging efforts, aimed at both PCE, very short period binaries and PCE, WD-coalesced binaries where the contrast ratio is even more favorable.

(6) The closer a planet's orbit is to a single star during the MS stage, the worse its prospects are for survival when the star ascends the RGB and AGB. For example, Mercury, with an orbit of 0.46 au, will be destroyed as the Sun evolves off the MS (e.g., Schroder & Connon Smith 2008). In contrast, our results suggest that, despite their dramatic evolution, close binary stars may in fact be better for the survival of planets compared to single stars. For example, the Kepler-38 CBP has an initial orbital separation similar to Mercury's, and yet it survives the evolution of a binary star with a total mass of ∼1 ${M}_{\odot }$. A CBP with an even smaller orbit—as small as 0.3–0.35 au—can survive the evolution of a 1.3–1.4 ${M}_{\odot }$ (total mass) central binary star, as is the case for the inner planet in Kepler-47. A planet with the same orbit around a single star of the same mass will have no chance for survival (Villaver & Livio 2007). Thus, despite being a disruptive stage for the binary star itself, the CE phase may in fact promote planet survivability.

(7) With the exception of Kepler-1647, the Kepler CBPs currently orbit their binary star hosts within a factor of 2 of the critical limit for dynamical stability (Holman & Wiegert 1999, HW hereafter), that is, ${a}_{\mathrm{CBP},0}(1-{e}_{\mathrm{CBP},0})/{a}_{\mathrm{crit},0}\sim 2$ (the factor for Kepler-1647 is ∼7). However, as the binaries shrink and lose mass during the CE and the planets migrate to larger, eccentric orbits, this ratio will change. To calculate a_crit,PCE—the PCE critical limits (for those scenarios where the binaries do not coalesce, i.e., α_CE  = 3/5/10)—we use Equation (3) from HW using a binary mass ratio μ = M_A,PCE/(M_A,PCE + M_B,PCE).²³ As the CE phases result in a variety of a_bin, a_CBP,PCE, and e_CBP,PCE, here we quote only the most constraining critical limits, where the respective separation between a_bin and a_CBP,PCE (1 − e_CBP,PCE ) is at a minimum. For ${\alpha }_{1}$, the limits for Kepler-38, Kepler-47, Kepler-64, Kepler-1647 (NTCP), and Kepler-1647 (TCP) in terms of a_CBP,PCE (1 − e_CBP,PCE)/a_crit,PCE are 3.3, 7.9, 5.3, 15.9, 14.7, and 29.5, respectively, indicating that the CBPs remain in a dynamically stable regime.

While it is beyond the scope of this study, we note that tidal evolution of the binary orbit prior to the CE phase may affect the CBP orbits. Specifically, as most of the planets are currently orbiting close to their host binaries, a pre-CE decay in a_bin means that strong mean motion resonances (MMR) may sweep over the planet, leading to eccentricity excitation and potential destabilization. As an example, a_bin of Kepler-38 decreases from 31.6 R_⊙ to 24.6 R_⊙ prior to the CE phase for TCP (see Table 4), indicating that the CBP's orbit crosses 6:1, 7:1, and 8:1 MMR. The MMR crossings are (a) 7:1, 8:1, and 9:1 for Kepler-47 (planet 1); (b) 7:1 through 11:1 for Kepler-64; and (c) 99:1 through 116:1 for Kepler-1647.

(8) "...It is interesting to consider the situation at a much earlier time in the past when the primary was near the zero-age main sequence...," note Orosz et al. (2012a) for the case of Kepler-38b: "The primary's luminosity would have a factor ≈3 smaller...." While Kepler-38b is currently closer to its binary host than the inner edge of the habitable zone (HZ) for the system, four of the Kepler  CBPs (Kepler-16, Kepler-47c, Kepler-453, Kepler-1647) reside in the HZ. As their host binaries evolve, however, the location of the HZ will change.

Unlike the case for single stars where both the HZ and a (surviving) planet's orbit will expand during the RGB and AGB stages (Villaver & Livio 2007), as a close binary star evolves through a CE stage the HZ can shrink while its CBP migrates outward. The shrinking is caused by the luminosity drop of the central binary as one of its stars rapidly transitions from the MS to a compact object. Thus a CBP residing in the HZ during the pre-CE binary will leave the zone after the CE phase (e.g., Kepler-1647), whereas a CBP that is initially interior to the HZ may migrate to a PCE orbit coinciding with the PCE HZ.

Consequently, it would be equally interesting to consider the configuration of a CBP system at a much later time in the binary evolution. For example, the current orbital configuration of Kepler-35b is such that the planet's orbit is internal to the HZ (Welsh et al. 2012). However, as its host binary evolves through the CE phase, for α_CE  = 1/3/5/10 the primary star evolves into a HeWD (while the secondary remains on the MS as a 0.81 ${M}_{\odot }$ star), and the CBP migrates to mode(a_CBP,PCE) = 0.84–1.0 au, mode(e_CBP,PCE) = 0.0–0.4 (see Tables 8). This migration places the planet near the conservative HZ of 0.98–1.77 au (for 1 M_⊕; Kopparapu et al. 2014), using the luminosity, 0.4 L_Sun, and temperature of the secondary star, 5200 K (Welsh et al. 2012), and ignoring the flux contribution of the WD.

5.1. Limitations

5.2. Conclusions

We thank the anonymous referee for the insightful comments that helped us improve this paper. The authors are grateful to Jarrod Hurley and Marten van Kerkwijk for valuable discussions. VBK gratefully acknowledges support by an appointment to the NASA Postdoctoral Program at the Goddard Space Flight Center. DT was supported by a postdoctoral fellowship from the Centre for Planetary Sciences at the University of Toronto at Scarborough and is grateful for additional support from the Jeffrey L. Bishop Fellowship. This work was supported in part by NSERC grants to RJ. We acknowledge conversations with Daniel Fabrycky, Nader Haghighipour, Kaitlin Kratter, Boyana Lilian, Jerome Orosz, and William Welsh.

• BSE: Binary star evolution code
• CE: Common envelope
• CBP: Circumbinary planet
• EB: Eclipsing binary
• NTCP: No tidal circularization path (tides "OFF" in BSE)
• PCE: Postcommon envelope
• RLOF: Roche-lobe overflow
• TCP: Tidal circularization path (tides "ON" in BSE)

• ${\alpha }_{M}$ : Common-envelope mass-loss rate
• ${\alpha }_{1}$: 1 ${M}_{\odot }$ yr⁻¹ mass-loss rate
• ${\alpha }_{0.1}$: 0.1 ${M}_{\odot }$ yr⁻¹ mass-loss rate
• α_crit: Critical common-envelope mass-loss rate
• α_CE: Common-envelope efficiency parameter
• a_CBP,0: Initial semimajor axis of the circumbinary planets
• a_CBP,PCE: Semimajor axis of the circumbinary planets at the end of the common envelope
• β: Ratio between initial and final mass of the system
• β_eject: Runaway ejection ratio between initial and final mass of the system
• e_CBP,0: Initial eccentricity of the circumbinary planets
• e_CBP,PCE: Eccentricity of the circumbinary planets at the end of the common-envelope phase
• ψ: Common-envelope mass-loss index
• ${T}_{\mathrm{CE}}$: Common-envelope mass-loss timescale