Exomoons in Systems with a Strong Perturber: Applications to α Cen AB

The secular forced eccentricity from the three-body problem can provide a simple analytical estimate on the region of stability around a planet orbiting one star of a binary star system. In such a system, the exomoon exists in a three-body hierarchy, where it orbits a planet that in turn orbits its host star. In addition, a secondary star orbits the center of mass of three-body system at such a distance so that the planetary orbit remains stable to produce a hierarchical four-body configuration.

Various aspects of the stability of hierarchical four-body configurations have been investigated in some detail (e.g., Milani & Nobili 1983; Liu & Gong 2019 and references therein.). Most approaches are based on Jacobi (quasi) integrals, c² E, where c is the angular momentum and E the total energy of the system (i.e., Marchal & Bozis 1982). Jacobi integrals are constants of motion in the circular restricted three-body problem related to zero-velocity curves. In order for a three-body configuration to be Hill stable, the value of c² E must be less than the critical value at the Lagrange point L1. If that is the case, the zero-velocity surface prevents a (massless) exomoon from escaping its orbit around the planet. Although this concept can be shown to yield rigorous and sufficient criteria, derived expressions for stability limits are often lengthy and cumbersome to use, especially in the general three- and four-body problems where the Jacobi integrals are not constant.

We consider the hierarchical four-body problem as two coupled three-body problems: 1) planet-moon-host star and 2) planet-moon barycenter-host star-stellar companion. In the first three-body hierarchy, the planetary semimajor axis a_p relative to the host star is much larger than the moon's semimajor axis a_sat relative to the host planet (i.e., a_p ≫ a_sat). Then, the radius of the Hill sphere around the planet depends on the planetary semimajor axis and masses through,

$$\begin{eqnarray}&&{R}_{{\rm{H}}}={a}_{{\rm{p}}}{\left(\displaystyle \frac{{M}_{{\rm{p}}}+{M}_{\mathrm{sat}}}{3{M}_{\star }}\right)}^{1/3},\end{eqnarray} \tag{ 1 }$$

where M_p, M_sat, and M_⋆ are the masses of the host planet, moon, and host star, respectively. The above definition corresponds to the Hill radius for a circular planetary orbit. In a more general formula a_p is replaced by r_p, which also includes the planetary eccentricity e_p. As a consequence, R_H depends on the position of the planet on its orbit with a minimum at

$$\begin{eqnarray}&&{R}_{{\rm{H}}}^{\prime} =(1-{e}_{{\rm{p}}}){R}_{{\rm{H}}}.\end{eqnarray} \tag{ 2 }$$

Since the planet, moon, and host star form a hierarchical triple, the semimajor axis of the moon around the planet remains approximately constant over time. To guarantee the Hill stability of a moon around the planet, we have to identify pathways for the moon to escape the Hill sphere around the planet, which typically occurs when the planet is at periastron.

Under an external perturbation, the Hill radius R_H' is modified as the planetary eccentricity evolves (Georgakarakos 2003), where a simple model for the planetary orbital eccentricity evolution in a hierarchical system dates back to Heppenheimer (1978). In essence, the eccentricity of the planetary orbit can be decomposed into forced () and free (η) components (Andrade-Ines & Eggl 2017). The former is determined by system parameters, such as the ratio of semimajor axes and the eccentricity of the stellar binary, and can be derived from initial conditions. In their simplest form the equations for the forced and free amplitudes of the eccentricity vector read as:

$$\begin{eqnarray}&&{\epsilon }_{{\rm{p}}}\approx \displaystyle \frac{5}{4}\displaystyle \frac{{a}_{{\rm{p}}}}{{a}_{\mathrm{bin}}}\displaystyle \frac{{e}_{\mathrm{bin}}}{1-{e}_{\mathrm{bin}}^{2}},\end{eqnarray} \tag{ 3 }$$

and

$$\begin{eqnarray}&&{\eta }_{{\rm{p}}}\approx | {e}_{{\rm{p}},{\rm{o}}}-{\epsilon }_{{\rm{p}}}| ,\end{eqnarray} \tag{ 4 }$$

where e_p,o and a_p are the initial eccentricity and semimajor axis of the planet p. The semimajor axis and eccentricity of the stellar binary are denoted by a_bin and e_bin, respectively. The maximum planetary eccentricity (e_p,max) can be estimated by adding the amplitudes of the forced and free components of the respective eccentricity vector,

$$\begin{eqnarray}&&{e}_{{\rm{p}},\max }\approx {\eta }_{{\rm{p}}}+{\epsilon }_{{\rm{p}}}.\end{eqnarray} \tag{ 5 }$$

In our four-body configuration, the influence of the secondary star induces variations in the planetary eccentricity as the the system evolves. As a result, we determine the smallest extent of the planet's Hill radius as

$$\begin{eqnarray}&&R{{\prime} }_{{\rm{H}},\min }={a}_{{\rm{p}}}(1-{e}_{{\rm{p}},\max }){\left(\displaystyle \frac{{M}_{{\rm{p}}}+{M}_{\mathrm{sat}}}{3{M}_{\star }}\right)}^{1/3}.\end{eqnarray} \tag{ 6 }$$

To become unbound, a moon's orbital energy h must exceed the gravitational potential U, where this can occur when ${r}_{\mathrm{sat}}\geqslant R{{\prime} }_{H,\min }/2$. We ignore the moon's eccentricity e_sat because the amplitude of its forced eccentricity is small, _sat ∝ a_sat/a_p ≪ 1. Thus, we can construct a simple stability criterion for moons around planets that orbit one star of a binary star system, namely

$$\begin{eqnarray}&&{a}_{\mathrm{sat}}^{\mathrm{crit}}=\displaystyle \frac{1}{2}R{{\prime} }_{H,\min }\approx \displaystyle \frac{1}{2}\left(1-{\epsilon }_{{\rm{p}}}-{\eta }_{{\rm{p}}}\right){R}_{{\rm{H}}},\end{eqnarray} \tag{ 7 }$$

where _p and η_p are the forced and free eccentricities of the planet (see Equations (3) and (4)) and R_H is the Hill radius for a circular planetary orbit. Moons on stable orbits would have semimajor axes ${a}_{\mathrm{sat}}\lt {a}_{\mathrm{sat}}^{\mathrm{crit}}$. Note that in Equation (7) ${a}_{\mathrm{sat}}^{\mathrm{crit}}$ is given in units of distance. For the remained of this article we will measure ${a}_{\mathrm{sat}}^{\mathrm{crit}}$ in units of R_H .

For moderate planetary eccentricity (e_p ≳ 0.3) N:1 mean-motion resonances (MMRs) between the planetary and satellite orbits add an additional component to the eccentricity evolution of the satellite. For example, the 8:1 MMR occurs for a_sat ∼ 0.36 R_H and the 9:1 MMR occurs for ${a}_{\mathrm{sat}}\sim 0.\bar{3}$ R_H following ${a}_{\mathrm{sat}}^{N}={(3/{N}^{2})}^{1/3}$ R_H (Kipping 2009a; Rosario-Franco et al. 2020). The libration width of MMRs grows with the perturber's eccentricity (i.e., e_p ), which can allow for the overlap of MMRs to destabilize a significant range in a_sat (Murray & Dermott 1999; Mudryk & Wu 2006; Morais & Giuppone 2012). From Mardling 2013, we calculate the libration width Δσ_N of the N:1 MMRs using the following:

$$\begin{eqnarray}&&\begin{array}{l}{\rm{\Delta }}{\sigma }_{N}=\displaystyle \frac{6{{ \mathcal H }}_{22}^{1/2}}{{\left(2\pi \right)}^{1/4}}\sqrt{\displaystyle \frac{{M}_{\star }}{{M}_{T}}+{N}^{2/3}\displaystyle \frac{{M}_{p}{M}_{\mathrm{sat}}}{{M}_{T}^{2/3}{\left({M}_{p}+{M}_{\mathrm{sat}}\right)}^{4/3}}}\\ \displaystyle \frac{\sqrt{{e}_{\mathrm{sat}}}}{{e}_{p}}\sqrt{1-\displaystyle \frac{13}{24}{e}_{\mathrm{sat}}^{2}}{\left(1-{e}_{p}^{2}\right)}^{3/8}{N}^{3/4}{{\rm{e}}}^{-N\xi ({e}_{p})/2},\end{array}\end{eqnarray} \tag{ 8 }$$

where ${{ \mathcal H }}_{22}=0.71$ is a scale factor from the spherical harmonic expansion of a pendulum model for resonance, M_⋆ is the mass of the host star, M_T is the total mass of the star-planet-satellite system, and $\xi ({e}_{p})={\cosh }^{-1}(1/{e}_{p})-\sqrt{1-{e}_{p}^{2}}$. Note that Δσ_N = 0 for e_sat = 0, which implies that circular satellite orbits are always stable in this model, and this is not the case, as noted by Mardling (2013). When calculating Δσ_N for circular orbits, we use e_sat = 10⁻³ or e_p = 10⁻³ to avoid the complication noted by Mardling or the singularity for e_p → 0, respectively.

Using an Earth-mass moon on a circular orbit around a Jupiter-mass planet, we produce Figure 1 to illustrate the locations of the N:1 MMRs for such exomoons. The host star is assumed to equal α Cen A in terms of its mass (1.133 M_⊙) and luminosity (1.519 L_⊙). The white curves in Figure 1 are determined through ${a}_{\mathrm{sat}}^{N}\pm {\rm{\Delta }}{\sigma }_{N}$ and the color code denotes the order N of the MMR. The black regions represent configurations between the MMRs and the light green (solid) curve in Figure 1 traces the first intersection points between adjacent MMRs and can be calculated numerically using a root-finding function for ${a}_{\mathrm{sat}}^{N}-{\rm{\Delta }}{\sigma }_{N}={a}_{\mathrm{sat}}^{N+1}+{\rm{\Delta }}{\sigma }_{N+1}$. An initial condition to the right of the light green curve becomes unstable due to the overlap in the MMRs and this process can be slow due to the eccentricity growth timescale of the satellite. However, in our four-body configuration the planetary eccentricity undergoes oscillations due to secular forcing by the secondary stellar companion. As a result, the light green curve can shift leftward and more efficiently destabilize prograde satellites. Retrograde satellites are less affected and can survive for a longer timescale at a much larger satellite separation (Henon 1970), but chaotic diffusion within the MMR overlap can lead to instabilities on longer timescales.

Zoom In Zoom Out Reset image size

N-body simulations provide a direct test of exomoon stability and we use the orbital integration package Rebound (Rein & Liu 2012; Rein & Spiegel 2015) because it is a versatile tool for the evolution of complex N-body systems. We use the WHFAST and IAS15 integration schemes within Rebound to evaluate the stability of exomoons (i.e., Rosario-Franco et al. 2020) on 10⁵ or 10⁶ yr timescales using an initial timestep that is 2.5% of the exomoon's initial orbital period. Since our initial timestep is very small compared to the period of stellar binary, the WHFAST integration scheme maintains sufficient accuracy. We verified this using the IAS15 adaptive integration scheme at large planetary eccentricity. The 10⁵ yr timescale is sufficient for the most cases because the typical orbital timescale for an exomoon is much shorter (∼0.1 yr) and the secular timescale for perturbations from the stellar companion is ∼10⁴ yr (e.g., α Centauri AB; Quarles & Lissauer 2018; Quarles et al. 2018a, 2019). Perturbations from the stellar companion on retrograde-orbiting exomoons are weaker, especially in the habitable zone of each star, and require the longer 10⁶ yr integrations. Our definition of stability is for an exomoon that begins on a circular, coplanar orbit to survive for the full integration time. Unstable conditions are those that cause the exomoon separation r_sat to exceed the host planet's Hill radius (r_sat > R_H ) or crash into the host planet (r_sat < R_p ).

Our simulations use the α Centauri AB (α Cen AB) binary to identify the extent of exomoon stability because the stellar masses, separation, and eccentricity are well known, which is important for our criteria in Section 2. Pourbaix & Boffin (2016) provide values for α Cen AB based on analysis of archival data from the High Accuracy Radial Velocity Planet Searcher at the ESO La Silla 3.6 m telescope. Pourbaix & Boffin find the masses of star A and B are 1.133 M_⊙ and 0.972 M_⊙, respectively, while the binary orbital semimajor axis and eccentricity are 23.7 au (17.66 arcsec) and 0.524, respectively. These values update the previous analysis Pourbaix et al. (2002) and find good agreement with estimates from asteroseismology (Lundkvist et al. 2014). In our simulations, the binary companion is coplanar, begins at periastron, and is nodally aligned (see Table 1). Our choice for the binary initial phase affects the stability of planetary orbits sufficiently near the stability limit (Quarles et al. 2018b, 2020a), where the initial separation of the host planets in this work are in a regime that is long-term stable considering a wide range of initial phases for both the binary and planetary orbit (Quarles et al. 2020a).

Since the existence of planets in α Cen AB is currently unknown (much less their atmospheric composition), we choose an initial condition close to the inner edge of the conservative habitable zone (Eggl et al. 2020). The habitable zone (HZ) is commonly defined as the region where liquid water could potentially exist on the surface of a rocky planet and depends, among other things, on the properties of the incident radiation and composition of the planet's atmosphere (Kasting et al. 1993; Kopparapu et al. 2014; Eggl et al. 2020). We place the planet at 1.232 au around star A and 0.707 au around star B. Terrestrial planets on circular orbits at these semimajor axes receive one solar constant's worth of radiative flux annually. For a planet orbiting α Cen A, the planet's semimajor axis a_p is sampled at 1.232 au, 2 au, 2.25 au, or 2.5 au as star A has a greater chance to host planets as far out as ∼2–2.5 au (Quarles & Lissauer 2016). For a planet orbiting α Cen B, the planet's semimajor axis begins at 0.707 au or at 2 au, respectively. In both cases, the planetary orbit starts coplanar (i_p = 0°) with the binary orbit, with zero obliquity ψ_p = 0°, and is nodally aligned (ω_p = Ω_p = 0°). Most of the observed circumstellar planets in binary systems are Jupiter-like giants (Schwarz et al. 2016). For the planetary mass M_p and radius R_p , we use values identical to Jupiter's in our orbital stability simulations. Prior studies of the satellites in the solar system (Estrada & Mosqueira 2006; Canup & Ward 2006) developed formation pathways akin to the terrestrial planets, but for a circumplanetary disk. Canup & Ward (2006) suggested that in situ formation of moons from circumplanetary disks around giant planets are limited to ∼10⁻⁴ M_p due to the estimated flux of infall material from the circumstellar disk, but this prediction awaits confirmation from exoplanet observations. We use an Earth-mass exomoon M_sat in our simulations due to the potential of tidal capture through scattering events and the flux of smaller infall material could be larger than was for the solar system due to the formation environment of the binary stars. However, we switch to an Earth-Moon analog for the planet-moon system when considering the possible effects of tidal evolution (see Section 3.4). This is motivated by observational constraints for α Cen AB (Zhao et al. 2018) and the comparisons that we make with the solar system.

For the exomoons, we explore prograde and retrograde orbits because both are ubiquitous within the solar system (Jewitt & Haghighipour 2007), where there could be a higher incidence of scattering during the early stages of planet formation (Quintana et al. 2002, 2007; Haghighipour & Raymond 2007). Since we are investigating long-lived exomoons, we expect that the tidal dissipation in the host planet to efficiently align the exomoon's orbit with the planetary equator (i.e., coplanar with the planetary and binary orbits) and circularize the exomoon's orbit. As a result, our numerical model begins the exomoons on circular (e_sat = 0) orbits that are nodally aligned ω_sat = Ω_sat = 0° with the planet and binary orbits. We vary the planet-satellite separation a_sat as a fraction of the planet's Hill radius starting at 0.25 R_H and up to 0.7 R_H in 0.01 R_H steps, where the 0.25 R_H is just interior to planet-satellite MMRs and 0.7 R_H is the theoretical limit for stability for retrograde small bodies in the circular restricted three-body problem (Quarles et al. 2020a).

Following Rosario-Franco et al. (2020), we evaluate 20 initial phases f_sat that are randomly drawn between 0^°–180^° from a uniform distribution for each exomoon semimajor axis a_sat. When constructing our initial conditions, we use the above-prescribed orbital elements and convert pairwise (e.g., [A, (p,sat),B]) to the Cartesian positions and velocities. The center of mass between the planet and satellite typically resides within the planet and thus our results are indistinguishable from a Jacobian construction of initial conditions. Table 2 summarizes the parameter ranges that we explore for the orbital stability simulations.

Table 2. Parameter Ranges Used in Orbital Stability Simulations

parameter	unit	range/value
Mp	M⊕	318
ap	au	$\sqrt{{L}_{\star }}$, 2, 2.25, 2.5
ep		0.00–0.60
Msat	M⊕	1
asat	RH	0.25–0.70
isat	°	0, 180
fsat	°	0–180

Download table as:  ASCIITypeset image

Previous studies (Henon 1970; Innanen 1979; Hamilton & Burns 1991; Domingos et al. 2006) show that the extent of stable orbits for retrograde-orbiting moons is substantially larger than for a prograde moon. Domingos et al. (2006) explored how the stability limit for retrograde satellites varied with the planet and satellite eccentricity. However, the simulations by Domingos et al. evaluated only a single initial phase (f_sat = 0°) to formulate the stability boundary. As a result, the empirically determined retrograde stability limit was actually the upper critical orbit (Dvorak 1986). This distinction to their study was discussed recently for the case of prograde-orbiting satellites (Rosario-Franco et al. 2020). Moreover, Quarles et al. (2020a) numerically expanded the Jacobi constant stability criterion (Eberle et al. 2008) to retrograde orbits and determined that the retrograde stability limit was a_c ≈ 0.7 R_H for the μ = 10⁻³ case studied by Domingos et al. (2006).

We perform simulations of retrograde Earth-mass exomoons hosted by the Jupiter-mass planet that orbits 1 au from its Solar-mass host star (see Section 3.1). The exomoon begins on a circular orbit, where the host planet's eccentricity is varied from 0–0.6 in 0.01 steps. We test initial exomoon semimajor axis values starting within 0.25–0.70 R_H, where the outer limit is the determined from the Jacobi constant criterion (Eberle et al. 2008; Quarles et al. 2020a). The planet and satellite orbits begin nodally aligned, where we evaluate 20 random values for the initial satellite orbital phase f_sat. Figure 2 illustrates the relation between the initial satellite separation a_sat and the planetary eccentricity e_p using a color code for the fraction of initial phases f_stab that are stable for the full integration timescale of 10⁵ yr. Since these simulations lack a stellar companion the longer 10⁶ timescale described in Section 3.1 is unnecessary, but will be used for retrograde systems in Section 3.3.

Zoom In Zoom Out Reset image size

The white cells in Figure 2 mark when all 20 trials either collide with or are stripped from the host planet. The black cells mark when all 20 trials are stable, where the colored cells show the transition. The black (dashed) curve marks the prediction from the fitting formula by Domingos et al. (2006), where the red (solid) curve denotes our fitted stability boundary using a linear function in the planetary eccentricity in units of the Hill radius R_H,

$$\begin{eqnarray}&&{a}_{{\rm{s}}{\rm{a}}{\rm{t}}}^{{\rm{c}}{\rm{r}}{\rm{i}}{\rm{t}}}={C}_{1}(1-{C}_{2}{e}_{p}),\end{eqnarray} \tag{ 9 }$$

and the gray curves show the uncertainty in our estimate. Equation (9) has a similar structure as Equation (7) where the planetary eccentricity is used directly (instead of the maximum eccentricity) because the planetary eccentricity is nearly constant in time for these Sun-Jupiter-satellite systems with negligible tidal interactions. Figure 2 demonstrates that the boundary from Domingos et al. (2006) and our boundary represents the upper and lower critical orbits (Dvorak 1986), respectively. Our stability boundary recovers the limit set by the Jacobi constant criterion (Eberle et al. 2008; Quarles et al. 2020a) for circular planetary orbits. The red (solid) curve is fit by the coefficients: C₁ = 0.668 ± 0.006R_H and C₂ = 1.236 ± 0.019. Using our outer stability limit coefficient C₁ is important because it significantly affects the limiting satellite mass (M${}_{\mathrm{sat}}\propto {C}_{1}^{13/2}$) through tidal migration (Barnes & O'Brien 2002).

Planets orbiting either component in a stellar binary experience strong perturbations, but the perturbation strength scales with the pericenter distance of the secondary star and the initial semimajor axis of the planet (e.g., David et al. 2003; Quarles et al. 2018b). We perform numerical simulations (see Section 3.1) of prograde and retrograde-orbiting satellites using the α Cen AB binary and consider host planets whose orbits begin at the inner edge of their respective HZs (1.232 au or 0.707 au) or closer to their respective stability limits (a_p ∼ 2–2.5 au). α Cen AB is representative of moderately wide binary systems (Raghavan et al. 2010; Moe & Di Stefano 2017) in terms of its eccentricity and binary orbital period, which makes it important as a case study. The distribution of binary eccentricity and orbital period is quite broad and we explore the possible impacts on our results in Section 4.2.

Starting with planets and moons in the respective HZ of their host star in α Cen AB, we numerically determine the outer stability boundary for satellites (Figure 3). We also plot two gray curves that represent the upper or lower critical orbit using the formalism derived in Section 2 in units of the planetary Hill radius relative to the host star (Table 2),

$$\begin{eqnarray}&&{a}_{{\rm{s}}{\rm{a}}{\rm{t}}}^{{\rm{c}}{\rm{r}}{\rm{i}}{\rm{t}}}={C}_{{\rm{s}}{\rm{a}}{\rm{t}}}(1-{\epsilon }_{{\rm{p}}}-{\eta }_{{\rm{p}}}){R}_{{\rm{H}}},\end{eqnarray} \tag{ 10 }$$

where C_sat is a scale factor for R_H in Equation (7). The coefficient C_sat corresponds to values determined for exomoons in single star systems (Domingos et al. 2006; Rosario-Franco et al. 2020) for prograde and retrograde orbits (see Table 3). For the solid gray curve representing the upper critical orbit, C_sat is 0.5 (prograde) or 0.95 (retrograde). For the dashed gray curve representing the lower critical orbit, C_sat is 0.4 (prograde) or 0.65 (retrograde). Due to the absolute value imposed in Equation (4), the gray curves have a peak at the forced eccentricity, _p. Similar to Figure 2, the color code in Figure 3 illustrates the transition between stable and unstable orbits. In each panel, the solid curve encapsulates nearly all the potentially stable cells, while the dashed curve marks the more conservative inner boundary. The curves fit reasonably well because the eccentricity growth of the satellite, which is the main driver for instabilities, grows secularly in this regime. Moderately eccentric planets (e_p ≥ 0.3) undergo eccentricity oscillations and MMR overlap causes significant eccentricity growth for the satellite. The light green (solid) curves in Figure 3 (lower panels) mark the boundary where MMR overlap begins in the host star-planet-satellite system. Those curves shift left and right as the system evolves.

Zoom In Zoom Out Reset image size

There are a handful of stellar binaries with a separation of ∼20 au that are known to host Jovian exoplanets near the stability limit (e.g., HD 196885AB, γ Cephei AB; Hatzes et al. 2003; Chauvin et al. 2007, 2011; Torres 2007; Correia et al. 2008; Fischer et al. 2009) with extensive dynamical studies focusing on the stability of the planet (Thebault 2011; Giuppone et al. 2011, 2012; Satyal et al. 2013, 2014). We perform numerical simulations of similar conditions for a putative Jupiter-like exoplanet orbiting either star in α Cen AB at 2, 2.25, or 2.5 au from the host star to investigate the stability of moons under more significant perturbations from the stellar companion. Zhao et al. (2018) may have already excluded Jupiter-sized planets at these distances using archival radial velocity measurements, but our results are scaled by the planetary Hill radius and would still be applicable for lower-mass planet-satellite pairs as long as M_sat ≪ M_p.

Figure 4 illustrates our results for prograde-orbiting satellites and is similar to Figure 3 (top row). Satellites can be stable at larger separations when the host planet orbits its star at the forced eccentricity _p (see Section 2) and is apsidally aligned (ω_p = ω_bin) with the binary orbit (Andrade-Ines & Eggl 2017; Quarles et al. 2018b). Apsidal alignment is important because it limits the degree of precession for the planetary orbit and minimizes the planet's eccentricity oscillations, where the maximum planetary eccentricity e_p can reach ∼0.2 when starting from a circular orbit (Quarles et al. 2018b). For an eccentric planetary orbit in the HZ of α Cen AB, e_p can also vary by ∼0.2 over a secular cycle (Quarles et al. 2018b), but the magnitude depends on the proximity to the forced eccentricity. At the maximum of the planetary eccentricity oscillation, the truncation of the planetary Hill radius is more significant and significant overlap with MMRs can occur, where an ejection of the satellite becomes more probable.

Zoom In Zoom Out Reset image size

Figure 4 demonstrates this effect through: 1) the gray curves marking the upper and lower critical satellite orbits and 2) the light green curves denoting the boundary for overlapping MMRs. The dashed (light green) curves represent a shift of the MMR based curves (solid) to smaller initial planetary eccentricities by about ∼0.2, which is representative of the maximum oscillation of the planetary eccentricity that permits planetary stability from previous studies (Quarles & Lissauer 2016; Quarles et al. 2018b). Figure 5 illustrates the planetary eccentricity oscillations (Figures 5(a)–(c)) for three different initial values of e_p, where the remaining initial conditions are drawn from Figure 4(c) and a_sat begins at 0.3 R_H. Figures 5(a) and 5(c) show a higher variation due to their relative distance from the planetary forced eccentricity (_p ∼ 0.08 when a_p = 2.25 au). Figure 5(d)–(f) demonstrate the evolution of the satellite's apocenter Q_sat (black dots) in response to the planet's eccentricity variation, where the gray (dashed) curve marks the lower critical orbit boundary, light green (solid) curve denotes the boundary for MMR overlap when e_sat = 0, and the magenta (dashed) shows the shift of the MMR overlap boundary at the maximum satellite eccentricity. The boundary for MMR overlap depends on e_p and e_sat (see Equation (8)), where increases in e_sat lead to shifts in the MMR overlap boundary to lower e_p values and increases in e_p (relative to _p) correlates with a higher likelihood of reaching the new boundary. Instability can then ensue once the satellite enters the MMR overlap region, but this process is chaotic (Mudryk & Wu 2006). Figure 6 illustrates similar calculations as Figure 4, but for retrograde-orbiting satellites. The overall area is greater due to the enhancement to stability from retrograde orbits and the possible truncation due to planetary eccentricity oscillations is less severe. Even though more of the parameter space is stable for retrograde satellites, the maximum satellite eccentricity is also larger and the possible existence of multiple satellites would constrain the parameter space further (Giuppone et al. 2013).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In our solar system, the lifetime of moons is significantly constrained by tidal interactions (Barnes & O'Brien 2002; Sucerquia et al. 2019; Lainey et al. 2020). While constraints have been placed on the long-term tidal evolution of exomoons in single stellar systems (Sasaki et al. 2012; Sasaki & Barnes 2014), this has yet to be established for exomoons in stellar binary systems. To obtain a full picture of orbital stability of exomoons in stellar binary systems, it is necessary to consider the contribution of planetary and stellar tides. In this section, we apply a secular constant time lag tidal model and evaluate the migration lifetimes of exomoons within the HZs of each stellar binary component in α Cen AB.

Analyses of observational surveys (Zhao et al. 2018) have largely excluded Jupiter-mass exoplanets from the HZs of α Cen AB, where the detection thresholds are 53 M_⊕ and 8.4 M_⊕ for star A and B, respectively. Additionally, the low tidal time lag in Jupiter-like planets (k_2,JΔt_J ∼ 10⁻² s or Q_J ∼ 10⁶; see Leconte et al. 2010) greatly reduces the outward migration of satellites and prevents us from placing a strong constraint using tidal migration on possible exomoons, although recent results for Saturn suggest substantially faster migration rates could be possible (Lainey et al. 2020). More complicated models for the tidal dissipation in gas dominated exoplanets could be employed (e.g., Guenel et al. 2014; Alvarado-Montes et al. 2017), but such considerations are beyond the scope of this work. As a consequence we use an Earth-Moon analog in our tidal model instead, due to its potential to place meaningful constraints on exomoons in α Cen AB, and because the Earth-Moon system has tidal parameters that are better understood. Similar to Section 3.3, the planetary orbits are sampled from the inner edge of their HZ, but values for exomoon's semimajor axis begin at 8.64 R_⊕ (three times the Roche radius of an Earth-Moon analog).

We implement a constant time lag secular model (Leconte et al. 2010; Hut 1981) and evaluate the tidal evolution up to 10 Gyr. This model calculates the secular changes to the semimajor axes (a_p and a_sat), eccentricities (e_p and e_sat), and mean motion (n_p and n_sat) averaged over an orbit. The model is scaled by the tidal Love number k₂ and the time lag Δ t, where the latter is proportional to (nQ)⁻¹ in the constant phase lag (i.e., constant Q) tidal models (Leconte et al. 2010; Piro 2018). A reasonable time lag for the Earth-Moon system is ∼100 s, and we assume Δt_p = 100 s, unless stated otherwise. The initial rotation period of the host planet is 5 hr, which is consistent with expectations from terrestrial planet formation (Kokubo & Ida 2007). The external perturbations from the secondary star on the exomoon are weak compared to the planetary gravitational interactions and a tidal model considering the interactions between the host star-planet-moon are sufficient. We use a new module called tides_constant_time_lag from reboundx (Baronett et al. 2021) to test this assumption, which is based upon the well-established code mercury-T (Bolmont et al. 2015). Our implementation of the tides_constant_time_lag module evaluates the instantaneous torques at each timestep using the ias15 integrator and the N-body simulations are evolved for 10 Myr, which is an appropriate timescale to see the effects on the satellite orbit due to tides. At the time of this writing, the tides_constant_time_lag module does not evolve the changes to primary body's spin due to outward tidal migration of a secondary body. Therefore, we update the spin rate of the primary body periodically using interpolated values of the host planet's spin rate from the afore mentioned secular model. Figure 7 illustrates the results from reboundx (solid black) and the corresponding secular model (dashed red), where the two approaches are in agreement to a high degree. In both models, the exomoon's eccentricity grows at similar rates and orbital energy exchanges between the planet and moon become possible. Since the N-body method is more computationally expensive (12 CPU-days), we use the secular model for the remainder of this section.

Zoom In Zoom Out Reset image size

The outward migration of a satellite varies with both the assumed time lag and with the mass of the satellite. We evaluate the outward tidal migration using a broad range in the time lag (10–600 s) and the satellite mass (0.001–0.05 M_p) over a 10¹⁰ yr timescale (i.e., main sequence lifetime of a G dwarf). Figure 8 shows variation in outcomes due to these parameters for an Earth-Moon analog orbiting either stellar host in α Cen AB. The color code delineates the assumed parameter in each panel. Figures 8(a) and (b) evolve a Moon-like satellite (0.0123 M_⊕). Satellites are stable for the full range of time lags considered for α Cen A. But, exomoons can escape after 1 Gyr when their host planet orbits α Cen B for a time lag ≳100 s. For α Cen B, the inner edge of the HZ is closer to the host star and the tidal despinning of the host planet is faster, which can accelerate the outward migration. Figures 8(c) and (d) consider a range of satellite masses, where the time lag is held fixed at 100 s. There is a smaller range of outcomes when the satellite mass is varied, where larger planet-satellite mass fractions (≳0.02 M_sat/M_p) shorten the despinning timescale for the host planet. Once the host planet is despun, the migration changes direction and the satellite begins to fall toward the planet. In this case, the stellar tide on the host planet is weak enough so that the infall rate is slow as is shown by the flattening of the 0.05 M_⊕ curves in Figures 8(c) and (d). If the host planet's orbit was ∼0.5× smaller, then a putative satellite would collide with the host planet within 10 Gyr (e.g., Sasaki et al. 2012; Sasaki & Barnes 2014).

Zoom In Zoom Out Reset image size

Previous studies have historically used tidal migration to place upper limits on satellite masses using a constant phase lag (i.e., constant Q_p) tidal model (Goldreich & Soter 1966; Barnes & O'Brien 2002; Domingos et al. 2006). However, these efforts are usually limited to cases where the satellite mass is much smaller than the planet (i.e., M_sat/M_p ≪ 1) so that the satellite mass is neglected in the formulation of the satellite's mean motion n_sat. Additionally, the final expression is determined in terms of a fraction of the planets Hill radius f (Barnes & O'Brien 2002, see their Equation 8). In the constant phase lag model, the satellite's semimajor axis changes via the following differential equation (Murray & Dermott 1999):

$$\begin{eqnarray}&&{\dot{a}}_{\mathrm{sat}}=\displaystyle \frac{3{k}_{2,{\rm{p}}}}{{Q}_{{\rm{p}}}}\sqrt{\displaystyle \frac{{{GM}}_{{\rm{p}}}}{{a}_{\mathrm{sat}}}}\displaystyle \frac{{M}_{\mathrm{sat}}}{{M}_{{\rm{p}}}}\sqrt{1+\displaystyle \frac{{M}_{\mathrm{sat}}}{{M}_{{\rm{p}}}}}{\left(\displaystyle \frac{{R}_{{\rm{p}}}}{{a}_{\mathrm{sat}}}\right)}^{5},\end{eqnarray} \tag{ 11 }$$

where G is the gravitational constant, R_p is the planetary radius, k_{2, p} is the planetary Love number, and Q_p is the tidal quality factor. When the initial planetary rotation period is sufficiently short (Piro 2018), the satellite migrates past the stability limit (Equation (10)) on a timescale T that is determined through the integral of Equation (11) resulting in:

$$\begin{eqnarray}\begin{array}{rcl}T & = & \displaystyle \frac{2}{13}\left({a}_{{\rm{f}}}^{13/2}-{a}_{{\rm{o}}}^{13/2}\right)\displaystyle \frac{{Q}_{{\rm{p}}}}{3{k}_{2,{\rm{p}}}{R}_{{\rm{p}}}^{5}}\\ & & \times \ {\left(\sqrt{{{GM}}_{{\rm{p}}}\left(1+\displaystyle \frac{{M}_{\mathrm{sat}}}{{M}_{{\rm{p}}}}\right)}\displaystyle \frac{{M}_{\mathrm{sat}}}{{M}_{{\rm{p}}}}\right)}^{-1}\end{array}\end{eqnarray} \tag{ 12 }$$

where a_f and a_o are the final and initial semimajor axes of the satellite, respectively. For moons that form close to the host planet, the term proportional to a_o can be neglected because a_f ≫ a_o. To determine the limiting mass ratio, we make the substitution x_m = M_sat/M_p and set a_f equal to fR_H, where f = C_sat(1 − _p − η_p) from Equation (10). Rewriting R_H in terms of x_m as ${R}_{{\rm{H}}}={a}_{{\rm{p}}}{[(1+{x}_{{\rm{m}}}){M}_{{\rm{p}}}/(3{M}_{\star })]}^{1/3}$ and inserting a_f into Equation (12) yields the following expression after straight forward algebraic manipulation:

$$\begin{eqnarray}&&\displaystyle \frac{{x}_{{\rm{m}}}}{{\left(1+{x}_{{\rm{m}}}\right)}^{5/3}}\leqslant \displaystyle \frac{2}{13}{\left[{{fa}}_{{\rm{p}}}{\left(\displaystyle \frac{{M}_{{\rm{p}}}}{3{M}_{\star }}\right)}^{1/3}\right]}^{13/2}\displaystyle \frac{{Q}_{{\rm{p}}}}{3{k}_{2,{\rm{p}}}{{TR}}_{{\rm{p}}}^{5}\sqrt{{{GM}}_{{\rm{p}}}}}.\end{eqnarray} \tag{ 13 }$$

Equation (13) represents the upper limit in the mass ratio x_m between the satellite and the planet for a given configuration after a set time T. For the range of parameters considered in this work Equation (13) has only one real root which is most easily found numerically using a root-finding algorithm, such as the root_scalar function in numpy (Harris et al. 2020). Equation (13) reduces to Equation (8) from Barnes & O'Brien (2002) in the low-mass ratio limit (i.e., x_m ≪ 1). As pointed out by Piro (2018), Equation (13) will hold as long as sufficient spin angular momentum in the host planet is available for the outward migration. Once the planetary rotation synchronizes with the satellite's mean motion, the outward migration stops and inward migration toward the Roche limit can occur.

We test the validity of this Equation (13) using an Earth-mass host planet in α Cen AB over 10 Gyr in a constant time lag model (Δt_p = 100 s) following the system setup in Figure 8, where the outer stability limit ${a}_{\mathrm{sat}}^{\mathrm{crit}}$ is truncated by the planetary eccentricity e_p. In addition, we evaluate a comparable constant Q tidal model for a more direct comparison. Figure 9 shows that Equation (13) (black curve) fits the boundary between stable (colored cells) and unstable (gray cells) for planet-satellite mass ratios ≲10⁻², where the assumption for only outward migration remains valid. Above this curve, there are solutions where the infall timescale (after synchronization) is longer than 10 Gyr. The color code denotes how far the satellite migrates outward over 10 Gyr relative to the respective critical semimajor axis and the ⊕ symbol denotes parameters for an Earth-Moon analog at the planet's forced eccentricity. The constant time lag models (Figures 9(a) and (b)) allow for a smoother transition when the migration direction changes (outward to inward migration) as compared to the respective constant Q models (Figures 9(c) and (d)). We find that Equation (13) is more applicable for planets orbiting α Cen B because its HZ is closer to the host star and when the planet-satellite mass ratio is ≲2–3%. Different assumptions on the dissipation parameter (Δt_p or Q_p) can affect the outcome, but those also run into limits for plausibility since we are considering Earth-like values estimated from the terrestrial planets in the solar system (Quarles et al. 2020b).

Zoom In Zoom Out Reset image size

The presence of exomoons for single star systems can be deduced using transit-timing variations (TTVs) and transit duration variations of the planet as it crosses its host star (Kipping 2009a, 2009b), where these measurements also depend on the planet-satellite mass ratio m_sat/m_p and separation a_sat. These factors affect the maximum TTV (rms) amplitude through the displacement of the host planet from the center of mass (Sartoretti & Schneider 1999), where the larger planet-satellite separations increase the TTV magnitude for a given M_sat/M_p. We examine a hypothetical case of an Earth analog near the inner edge of the conservative HZ of each star in α Cen (see Figure 3) and how the presence of the stellar companion affects the maximum TTV amplitude for the putative transiting planet due to a prograde-orbiting satellite. The Earth analog is assumed to be at its forced eccentricity to set the outer stability limit.

The binary companion primarily affects potential exomoons through the truncation of the outer stability limit. Planetary tides could also limit the allowed parameters, but the tidal dissipation could be plausibly adjusted to allow for long-lived exomoons and thus cannot place a meaningful constraint for Earth analogs in α Cen AB. However, satellites of Earth analogs in other binary systems with less luminous secondary stars can be meaningfully constrained (Sasaki et al. 2012; Quarles et al. 2020b). Figure 10 shows the parameters that can produce 1, 2, 5, 10, 20, and 40 minutes TTVs (color coded) induced by prograde or retrograde satellites. An analog of the current Earth-Moon system in terms of mass ratio and separation is denoted by the ⊕ symbol, which would produce a ≳2 minutes TTV. Much larger TTVs (≲40 minutes) are possible for retrograde-orbiting satellites because the outer stability limit is farther from the host planet. An instrument similar to the Kepler Space Telescope is ideal for observing circumstellar planets and measuring potential TTVs (Ford et al. 2011). The Transiting Exoplanet Survey Satellite (Ricker et al. 2016) has already detected a transiting planet within a stellar triple consisting of M-stars (LTT 1445ABC; Winters et al. 2019), however the planet is likely already tidally locked which may limit its chances to host long-lived moons (Sasaki et al. 2012).

Zoom In Zoom Out Reset image size

Although an Earth analog in α Cen AB is a good case study, results based on the latter may not translate to other populations of binary star systems. In this section we, therefore, examine the frequency of Earth analogs hosting moons in binary systems more generally through a Monte Carlo experiment. ⁷ In our previous work (Quarles et al. 2019, 2020a), we used empirically derived probability density functions (PDFs) from surveys of binary stars with a solar-type primary (Raghavan et al. 2010; Moe & Di Stefano 2017). The PDF for the binary period distribution is a log-normal distribution for the binary period P in days, ${p}_{\mathrm{log}{\rm{P}}}\propto {e}^{-{(\mathrm{log}{\rm{P}}-\xi )}^{2}/2{\sigma }^{2}}$ where ξ = 5.03, σ = 2.28, and 4 ≤ log P ≤ 7. For the mass ratio q = M_B/M_A, we use a broken power-law PDF ${p}_{q}\propto {q}^{{\gamma }_{n}}$, where γ₁ = 0.3 for q ≤ 0.3 and γ₂ = − 0.5 for q > 0.3. Following Moe & Di Stefano (2017), we add to the PDF an excess twin fraction of 0.1 for q ≥ 0.95 to account for the observed stellar twins. An additional power-law PDF ${p}_{e}\propto {e}_{\mathrm{bin}}^{0.4}$ is included to account for the binary eccentricity. Coefficients for these PDFs are numerically determined so that the total probability is equal to unity, ∫p_x dp = 1. In contrast to our previous work, we ignore the uncertainty in power law exponents as they do not substantially alter the results.

To generate random binary system parameters, we numerically determine the cumulative distribution function for each of the PDFs described above, use a draw from a uniform distribution ${ \mathcal U }(0,1)$, and then numerically determine the random variate by parallel array matching. The planetary semimajor axis a_p is calculated at 1 Earth Flux (S_⊕ = 1L_⊙/(1 au)²) and can vary depending on which star is hosting the planet. As the PDF for the binary mass ratio varies, we adjust the secondary star's luminosity using the common power-law approximation of the mass–luminosity relation for main sequence stars. We also assume that the planet begins at its forced eccentricity (see Equation (3)). Our final results use the semimajor axis ratio a_p/a_bin, where the planetary semimajor axis is normalized by the binary semimajor axis a_bin. We calculate the maximum satellite mass (${M}_{\mathrm{sat}}^{\max }$) assuming an Earth-Moon analog for tides with a 10 Gyr lifetime (see Section 3.4) and normalize by the planetary mass M_p = M_⊕. We choose 50 random values for the satellite mass and semimajor axis using uniform distributions within the ranges $-3\,\leqslant \mathrm{log}{M}_{\mathrm{sat}}\leqslant \mathrm{log}[{M}_{\mathrm{sat}}^{\max }]$ and $0.05{R}_{{\rm{H}}}\leqslant {a}_{\mathrm{sat}}\leqslant {a}_{\mathrm{sat}}^{\mathrm{crit}}$. We then use the interpolation map technique (Quarles et al. 2020a) to verify that the generated host planet orbit is stable and our formalism from Section 2 to determine the outer stability boundary for the exomoon ${a}_{\mathrm{sat}}^{\mathrm{crit}}$. Once the stability of the planetary and exomoon orbits are validated, the TTV (rms) amplitude is calculated (Kipping 2009a, see their Equation A27).

Overall, we perform 10,000 draws from the PDFs representing the binary system parameters and 50 draws of the exomoon parameters (mass and separation) assuming an Earth-analog orbiting the primary star with a prograde-orbiting satellite. The corresponding results are shown in the top left panel of Figure 11. The process is repeated under the same condition, but with a retrograde-orbiting satellite and prograde/retrograde satellites for an Earth-analog orbiting the secondary star to generate the results presented in the other three panels. Figure 11 illustrates the estimated number of Earth-mass planets that could host satellites for a given TTV assuming a Solar-mass primary star. The highest number of systems occurs for widely separated binary systems (a_p/a_bin ∼ 0.001) and low-TTV amplitudes (∼0.3 minutes). Higher TTV amplitudes are possible, but a large survey would likely be necessary to ensure a detection. A retrograde-orbiting exomoon around the primary star could produce a larger TTV amplitude (∼30 minutes) due to the a larger outer stability limit, while the maximum TTV amplitude in the other cases is ≲10 minutes. Counterintuitively, planets orbiting the secondary star (B) have a slightly better chance of hosting exomoons because a distance of 1 Earth Flux is closer to the host star and less likely to be significantly perturbed by the more massive star A. These results scale with the primary star mass, where more massive primaries (1.2 M_⊙) allow for larger TTV amplitudes and less massive primaries (0.8 M_⊙) restrict the parameter space due to more significant tidal interactions between the exomoon and host star. The primary star of α Cen AB would be a good candidate for searching for TTV inducing exomoons if transiting Earth analogs were present. However, surveys of α Cen AB for planets are difficult because of pixel saturation in photometric observations (Demory et al. 2015) and astrophysical noise in radial velocity observations (Zhao et al. 2018).

Zoom In Zoom Out Reset image size