Orbital Stability and Secular Dynamics of the Proxima Centauri Planetary System

The two innermost planets of the Proxima Centauri system are separated by just 0.02 au, inducing strong gravitational interactions between them. We assess this interaction by leveraging fast orbital stability indicators and find that orbital stability is very likely if the initial eccentricities of planets b and d are less than ∼0.2, but cannot confirm stability at larger values. We find that stability is not strongly affected by the true masses of the planets or by the distant planet c. However, mutual inclinations between 95° and 142° often result in unstable motion. We further explore the long-term evolution of the orbits in these stable regions of parameter space and find that circularization can take over 5 Gyr. This tidal evolution could support surface energy fluxes in excess of 1 W m−2 for over 1 Gyr, possibly affecting planet b’s habitability.


INTRODUCTION
Proxima Centauri is the nearest stellar neighbor to the Sun and is known to host at least 3 planets, including one ∼Earth-mass planet in the habitable zone (Anglada-Escudé et al. 2016;Damasso et al. 2020;Suárez Mascareño et al. 2020;Faria et al. 2022).Given the system's proximity and the possible habitability of the middle planet, b, this system is uniquely valuable in the study of exoplanet formation, evolution, and habitability.Herein we describe a series of experiments to test the orbital stability of the Proxima Centauri system over a range of plausible orbital architectures.
The orbital evolution of Proxima's planetary system could affect the habitability of Proxima b in many ways.Orbital oscillations could affect the seasons and trigger snowball states (Spiegel et al. 2009;Armstrong et al. 2014;Deitrick et al. 2018b;Wilhelm et al. 2022), while tidal damping can drive long term decay of orbital semimajor axis and eccentricity (Barnes et al. 2008;Barnes 2017).Additionally, tidal friction in Proxima b's interior could create large amounts of energy that affect the Corresponding author: Joseph R. Livesey jrlivesey@wisc.edugeneration of magnetic fields (Driscoll & Barnes 2015), vigorous volcanic activity (Jackson et al. 2008; Barnes et al. 2009), and possibly even drive a runaway greenhouse (Barnes et al. 2013).Thus, defining the physically plausible range of orbital architectures in this system could significantly improve our understanding of Proxima b's habitability.
Not surprisingly, then, this study is not the first to consider the stability of this system.In particular, Meng et al. (2019) estimated the parameter space that permits orbital stability, but they used an older architecture that included only planets b and c.In this study, we update their results by considering the currently observed planetary system.
In particular, the recently detected planet d (Faria et al. 2022) demands a reanalysis of this system's stability.The orbits of Proxima d and b are sufficiently near so that gravitational perturbations between the two could drive the two into instability, i.e., the ejection of a planet from the system.Such cases should therefore be deemed implausible starting conditions for the system, i.e., evolutionary models can safely avoid those unstable parameter combinations.Using semi-analytic and N -body simulations, we systematically searched for unstable configurations of this system.We also identify architectures that give rise to chaotic dynamics, which do not necessarily imply instability.Identifying these chaotic configurations is useful for investigating longterm behavior, as the dynamics of systems exhibiting orbital chaos are ill-described by a disturbing function containing only secular terms (e.g., that employed in DistOrb; see §2.5 and Barnes et al. 2020).
In §2, we outline the theoretical models we use to assess stability and the parameter space we will explore.In §3, we present stability maps of the system obtained through thousands of simulations varying the initial conditions.In §4, we discuss the major outcomes of this work and implications for future studies of this planetary system.Note that code for reproducing the results presented in this paper is available online. 1

METHODS
In this section we describe theories of orbital stability, the numerical methods we employed to identify unstable cases, and the range of initial conditions we explored.Throughout this study we consider only the three-body system consisting of the host star and planets d and b.We tested the role of planet c, with an orbital period of ∼ 1,900 days, and found that it has negligible impact on the stability of the system, so we will not present detailed results of its role.
Proxima Centauri is a member of the α Centauri triple star system, orbiting the binary α Centauri AB at a distance of about 8700 au (Kervella et al. 2017).Asteroseismological analysis has revealed that the binary is 5.3 ± 0.3 Gyr old (Joyce & Chaboyer 2018).We assume that Proxima Centauri formed in the same cluster as α Centauri AB and therefore has the same age.While Proxima's planetary system could be modified by the binary (Barnes et al. 2016), we will ignore that possibility in our studies.

Stability Criteria
In this study we evaluate the Lagrange stability of the Proxima system.In a Lagrange stable planetary system, all planets remain gravitationally bound to the star and the ordering of the planets is preserved.There is no analytic formula to determine whether a system is Lagrange stable, and so we resort to N -body simulations.
We first consider, however, the possibility that the three-body system is Hill stable (Szebehely & McKenzie 1977;Marchal & Bozis 1982;Gladman 1993;Barnes & Greenberg 2006), which is determined analytically, but can only approximate orbital stability.The definition of Hill stability is that the ordering of the bodies 1 https://github.com/jrlivesey/ProximaStabilitycannot change.In other words, if the outer planet escapes to infinity, the system would still be Hill stable.Although this definition is very limited, it does reduce to a single equation that can be quickly evaluated.In the present study, we assess the agreement between the "Hill stability boundary" of the system in orbital element space and the stability boundary obtained through simulations.Such an analysis has been performed previously for systems of two ∼Jupiter-mass planets (Barnes & Greenberg 2006, 2007), but not, to our knowledge, for systems of terrestrial planets.
A system is guaranteed to be Hill stable if it satisfies the inequality β > β crit for and (2) where m ⋆ , m 1 , and m 2 are the masses of the star, inner planet, and outer planet, respectively; M = i m i is the total mass of the system; M = (1/2) i̸ =j m i m j ; C is the total orbital energy of the system; and L is its total angular momentum (Marchal & Bozis 1982).A three-body system for which β/β crit < 1 is likely but not certain to be Hill unstable.
To evaluate Hill stability, we generated 10 6 possible orbital configurations and calculated β/β crit .We randomly sampled the longitudes of ascending node, arguments of periastron, and mean anomalies of both planets from U[0, 360) • .The inclinations of the planets were sampled from U[0, 15] • .The inclination of the planets relative to the sky plane, used to calculate their true masses, was sampled from U[10, 90] • .For both of the orbital eccentricities, we used e ∼ U[0, 0.9].We assumed the orbits were coplanar for these calculations.

MEGNO
A hybrid approach to estimating orbital stability is to compute a quantity called the MEGNO (Mean Exponential Growth of Nearby Orbits; Cincotta & Simó 2000), which measures the divergence of two trajectories in phase space.If δ(t) denotes the separation between two initially infinitesimally close particles in phase space, then the MEGNO Y is defined as and its time-averaged value as For periodic (stable) evolution, ⟨Y ⟩ → 2 as t → ∞.For a chaotic system, ⟨Y ⟩ → (γ/2)t, where γ is the maximum Lyapunov characteristic number (Cincotta & Simó 2000;Goździewski et al. 2001;Morbidelli 2002;Cincotta et al. 2003;Hinse et al. 2010).The MEGNO is defined similarly to the Lyapunov characteristic number, but is weighted with time, amplifying stochastic behavior in the evolution and allowing for an earlier detection of chaos (Cincotta & Simó 2000;Goździewski et al. 2001;Morbidelli 2002).We follow standard practice, and label the motion as chaotic if |⟨Y ⟩ − 2| > 10 −2 , indicating the system may be unstable on Gyr timescales.

General Relativistic Effects
Relativistic (GR) effects can modify the precessional frequency of the orbit by an amount where n is the mean motion of the planet in question (Sterne 1939).If we apply this relation to the most extreme case studied here, i.e., that with e = 0.9 for both planets, we obtain GR precession periods of about 22,000 years and 82,000 years for planets d planet b, respectively.Since previous work has shown that N -body integrations up to 10 6 orbits of the outermost body are sufficient to identify most unstable configurations of a planetary system (Barnes & Quinn 2004).We therefore run our simulations for about 10 6 orbits of planet b, which is ∼ 31,000 yr.Therefore, GR will only induce one extra circulation over the length of our simulation and should therefore not be a significant effect.To decrease computational cost, we therefore do not include a post-Newtonian correction to our model.By the same token, we may safely assume that the rate of precession due to the Lense-Thirring effect is also negligible (Will & Nordtvedt 1972;Iorio 2008).

N-Body Simulations
Table 1 lists the current best fits and uncertainties for the physical and orbital elements of the system.We assume the host star has a mass of 0.12M ⊙ (Kervella et al. 2017).However, we are also interested in the range of plausible architectures at the end of planet formation as tidal damping can cause long-term changes to the planetary orbits (e.g., Barnes et al. 2008;Barnes 2017;Meng et al. 2019).We will thus consider a wider parameter space than allowed by observational uncertainties, focusing on larger eccentricities and semi-major axes.
To gauge the dependence of the system's stability on the system architecture, we performed three sets of 10,000 N -body simulations with the publicly available REBOUND code 2 (Rein & Liu 2012;Rein & Tamayo 2015;Rein & Spiegel 2015;Rein & Tamayo 2016).In Set I we varied the initial eccentricities of both planets; in Set II we varied the planets' masses; and in Set III we varied the initial inclinations with respect to the fundamental plane.
In addition to evaluating the MEGNO, in our Set I simulations we identify the parameter space in which the orbits cross, and in which we therefore anticipate chaos.As there are two planets included in our Nbody simulations, planet d's apocenter coincides with planet b's pericenter, and thus the orbits cross where We also incorporate the two-planet chaos criterion given by Hadden & Lithwick (2018) in terms of the relative eccentricity, which in our case simplifies somewhat to In our Set II simulations, we test the stability of the system with combinations of masses up to 10M ⊕ for both planets.Previous studies have specified a critical orbital separation in a two-planet system with low e and i, at which the system becomes unstable.The critical separation is ∆a = 2 √ 3R H , where is the Hill radius (Hasegawa & Nakazawa 1990;Gladman 1993;Chambers et al. 1996).Primed quantities refer to the outer body.In the most extreme case considered here, both masses are 10M ⊕ , with which ∆a = 36.3RH .We therefore expect every experiment in Set II to result in stable motion; we perform them and include the results here for completeness.
The mutual inclination (or relative inclination) Ψ between the two planets' orbits is given by cos Note that because there are only two planets included in our N -body simulations, cos(Ω − Ω ′ ) ≃ −1 for our purposes.The resulting distribution of simulated cases is weighted toward higher values of Ψ, since nearly coplanar, prograde orbits are already expected to be stable.For each set of simulations, we identify the Hill stability boundary in the relevant parameter space as an additional analytic stability criterion.We compute this boundary for each set because the Hill stability parameter depends upon the eccentricities, masses, and inclinations of both planets according to Eq. ( 1).
The simulations are integrated for ∼ 10 6 orbits of Proxima b.This timescale is longer than the conventional integration times used for identifying chaotic motion with the MEGNO, which are typically ∼ 103−4 times the longest dynamical timescale in the system (Cincotta & Simó 2000;Hinse et al. 2010;Meng et al. 2019).We perform our integrations over this longer timescale because we find some cases that are surprisingly stable (see Section 3.1).
The initial conditions for all of our REBOUND simulations are provided in Table 2.The semi-major axes used are slightly greater than the measured values to account for tidal effects over time.The approximate mass values used in these simulations for both planets were obtained by setting i obs = 133 • , the mean observational inclination of Proxima c obtained astrometrically by Benedict & McArthur (2020).
We performed two varieties of N -body simulations with REBOUND.We first used WHFAST, a symplectic Wisdom-Holman integrator (see Rein & Tamayo 2015;Wisdom & Holman 1991) to integrate the simulations.WHFAST is fast and efficient, but does not fare well when the planets' orbits are strongly perturbed (i.e., when they become significantly non-Keplerian).In cases where the relative energy error exceeds 10 −6 , we re-run the simulations using IAS15, a 15th order Gauss-Radau integrator (Rein & Spiegel 2015).IAS15 employs adaptive timestepping to resolve close encounters.At the end of each simulation, we find ⟨Y ⟩ using built-in functions in REBOUND that employ variational equations to calculate the trajectories of nearby orbits (Rein & Tamayo 2016).

Secular Orbital Evolution
The N -body experiments discussed in the previous section yield a set of points in initial parameter space that give rise to quasi-periodic motion, where ⟨Y ⟩ → 2, and are thus suitable starting points for long-term simulations of the secular evolution.
We performed six such simulations using the modular planetary evolution code VPLanet 3 (Barnes et al. 2020).
In particular, we used the EqTide and DistOrb modules to simultaneously simulate evolution of the orbits due to equilibrium stellar tides (the "constant phase lag" model of Goldreich & Soter 1966) and due to the disturbing function parameterizing gravitational interactions between the planets (Murray & Dermott 1999).For detailed descriptions and validation of these models, the reader is directed to Barnes et al. (2020) and Deitrick et al. (2018a).
Each body in the system is given a value of k 2 and Q. k 2 is the Love number of degree 2 and Q is the "tidal quality factor" that parameterizes the rate of tidal dissipation in a body's interior.Earth has Q = 12 (Williams et al. 1978), but a terrestrial planet without large oceans may have a quality factor an order of magnitude larger.In each of our simulations, k 2 = 0.5 for the star, k 2 = 0.3 for the planets.The star's Q is set to 10 6 , while the planet's Q is sampled in the range U[10, 10 3 ].Planetary radii are necessary for calculating the tidal torques, and are determined according to the power law R ∝ m 0.274 (Sotin et al. 2007).We can only estimate these radii, as no planet in this system is transiting (Jenkins et al. 2019).The values we use fall within the range established by the examination of Proxima b's internal structure given in Brugger et al. (2016).In addition to evaluating the time dependence of the orbital elements, we calculate the tidal heating of planet b throughout each evolution according to Ėtide = −( Ėorb + Ėrot ), where Initial conditions were sampled randomly from the permitted region of parameter space using VSPACE, a tool for initializing VPLanet's parameter sweeps.They are presented in Table 4.In these simulations, we include planet c and assume for the sake of simplicity that Proxima c's orbit begins with zero inclination and does not evolve tidally.Therefore, values of R, Ω, and Q are not provided for this planet.

Orbital Stability
We find that the orbits of the two inner planets are Hill stable in just 7.2% of the parameter space we explored.The Hill stable configurations are those with the lowest eccentricities for both planets and the boundary is shown with the black curve in Figure 1.
The results of Set I are also shown, with colors denoting values of ⟨Y ⟩.Blue points indicate periodic (stable) configurations (|⟨Y ⟩ − 2| ≤ 10 −2 ), white points are chaotic (|⟨Y ⟩ − 2| > 10 −2 , but both planets remained bound), and red points are cases in which one of the planets escaped to infinity.A clear trend is present in which larger eccentricities lead to decreasing stability.There are many cases in which ⟨Y ⟩ < 2 -or even becomes negative -at the end of a simulation.Such trajectories lie in the neighborhood of stable periodic orbits   (e.g., Breiter et al. 2005).There are cases in which the variational equations fail to converge at the end of a simulation, in which case REBOUND returns no value for the MEGNO.This outcome means that the separation in phase space has grown very large, and thus we label these cases as chaotic without having an exact MEGNO value.
We find that most periodic configurations lie within the Hill stability boundary, with e d ≲ 0.25 and e b ≲ 0.2.Most configurations with e d ≲ 0.55 and e b ≲ 0.45 give rise to motion that is chaotic, but does not become unstable within the length of these simulations.Beyond these boundaries, the system is very likely to be unstable, with one of the planets escaping to infinity within the length of these simulations.
We also find "islands" of configurations at higher initial eccentricities that are mostly classified as chaotic, but with some periodic results, e.g., near (0.55, 0.4) and (0.85, 0.55).The observational uncertainties suggest the system cannot currently exist in these Hill unstable islands (see Table 1).We return to this possibility in §4.
The outlying stable regions shown in Figure 1 could correspond to mean motion resonances, however, we find no evidence that MMRs are responsible for stabilizing these configurations.For a selection of simulated high-e stable cases, we calculated the time evolution of resonant arguments, which are for the j 1 : j 2 MMR, where the coefficients satisfy the d'Alembert relation i j i = 0. We found no libration in ϕ up to j 1 = 30.We also re-ran five cases for 1 Myr, and found that all cases initially designated as periodic were reclassified as chaotic (i.e., ⟨Y ⟩ ̸ → 2).We therefore conclude that these outlying islands are in fact unstable and that the MEGNO method was not able to identify instabilities over the first 10 6 orbits.This conjecture is supported by the fact that these anomalous cases lie in the crossing orbits region of parameter space, and thus chaotic evolution is expected.Note that simulating the system for its age would require at least 5,000 CPU hours per simulation with WHFAST on a modern workstation, and so is currently intractable.
Figure 2 shows the eccentricity evolution for three example cases: a periodic case, a chaotic case, and an unstable case.The locations of these cases in Figure 1 are marked by stars.Recall that Set I assumes coplanar orbits, so there is no evolution with respect to the inclinations.The stable case shows classic periodic motion in which both eccentricities oscillate with just one frequency.The chaotic case appears regular, but contains at least two frequencies, probably resulting from the larger eccentricities.The unstable case ejects planet d in less than 250 years.
Next we turn to Set II, which focused on the role of planetary masses on stability.In each of these simulations, the initial orbits are circular and coplanar.We found that every configuration we tested resulted in periodic motion.We conclude that stability is not strongly influenced by the planetary masses, which is consistent with the m 1/3 dependence in the Hill stability criterion.
Finally, we present results for Set III, which focused on the inclinations.As shown in Figure 3, we find that the motion is periodic for Ψ < 50 • and Ψ > 150 • .With few exceptions, systems in between evolve chaotically.80% of systems with 95 • < Ψ < 142 • resulted in ejections.This finding is consistent with previous results that show stability increases as the orbits approach coplanarity (e.g., Meng et al. 2019).Our obtained distribution of chaotic configurations in mutual inclination space, with a wide chaotic region in the middle skewed toward greater Ψ, agrees with the stability analyses of two-planet systems performed in Veras & Ford (2010) in terms of the maximum eccentricity variation.Our Set III data are presented again in two dimensions -i d and i b -in Figure 4. Figure 5 shows the eccentricity evolution for three example cases of various initial mutual inclination.

Long-term Orbital Evolution
Using initial orbital configurations that give rise to regular motion, we simulated the long-term orbital evolution of the system.One result is shown in Figure 6.Notably, in this scenario tidal evolution persists longer than the 7 Gyr for which the simulation was run, and much longer than the expected value for the age of the system.
The evolution of Proxima d and b in terms of e and a in each of the VPLanet simulations is displayed in Figure 7. Five of the six secular evolutions result in tidal heating that remains above the value for modern Earth for up to 7 Gyr.In the remaining simulation, Proxima b has Q ≃ 16.This result suggests that if Proxima b has Q ∼ 10 (on the order of Earth), then the surface flux of tidal heating is not significant after the first few Gyr because the orbit has circularized.If, however, Q ∼ 10 2 , then tidal interactions render a first-order effect on the heat delivered to the world's surface over its entire history.The time-averaged MEGNO for this case diverges significantly from 2, signifying that the motion is chaotic.At least 3 distinct frequencies are seen to be present in both planets' eccentricities, which is greater than the number predicted by secular theory.An unstable case originating from e d = 0.45 and e b = 0.3 is shown in panel (c).The ejection of planet d from the system is evident at t = 222 yr, where e d > 1.
We found cases in which the system is not Hill stable, and yet exhibits quasi-periodic motion.However, most quasi-periodic configurations of the system reside within the Hill stability boundary, as we should expect.Our analysis also identified regions where our simulations did not produce ejections, despite having crossing orbits initially.Although these configurations appear stable, a subset of simulations re-run for a longer integration time results in instabilities, suggesting that most of the region where β < β crit is in fact unstable.
The separated islands of stable motion in Figure 1, however, are very likely inaccessible.The large region of instability separating them from the low eccentricity "mainland" is probably too large for the system to successfully cross.As shown in Figure 2, instabilities can arise within a few centuries, which is short compared to the tidal damping timescale.We thus conclude that Proxima b and d could only have formed in the low eccentricity region of stability.
Our simulations of secular orbital evolution in the system reveal that, beginning from the "stability mainland" of parameter space determined with our N -body integrations, the tidal orbital evolution of Proxima b can easily continue beyond 5 Gyr.In each case simulated with VPLanet, we see the eccentricity behaves as a damped harmonic oscillator, with secular amplitude decreasing rapidly with time.As such, we find no new constraint on the present-day eccentricity of this world, which has a value at +1σ of 0.06 according to Faria et al. (2022).This result suggests that tidal heating in Proxima b is likely ongoing, providing an additional source of energy to drive surface and atmospheric processes, perhaps including active biology.On the other hand, the increased volcanism could outgas water that is ul- timately destroyed by high energy photons, resulting in a desiccated, and therefore uninhabitable world (Luger & Barnes 2015).Regardless, our research shows that tidal heating should not be ignored for planet b unless the upper limit of its eccentricity can be significantly reduced.
Models of the evolution and habitability of Proxima b can use our results to target their studies in the physically allowed regions we have identified in Figures 1 and  4. Climate studies (e.g.Turbet et al. 2016;Del Genio et al. 2017;Colose et al. 2021) could consider eccentricities for Proxima b up to 0.45, but should probably focus on the blue triangle in the lower left corner of Figure 1.Studies of the planetary interior could also consider the role of tidal heating over the same range (see e.g., Henning et al. 2009;Barnes et al. 2013;Driscoll & Barnes 2015;Barnes et al. 2016).These results may also prove useful in dynamical modeling of the Proxima system from radial velocity data, providing limits on the maximum achievable orbital eccentricities and inclinations of planets d and b.     (Moore et al. 2007;Munk & Wunsch 1998).
dissipated from the orbit and rotation of the planet, respectively. 4Here, n and ω are the orbital and rotational frequencies, ζ = sgn(2ω − 2n), and

4
Our Equations 9 and 10 are simplified from the CPL tidal heating equations inBarnes et al. (2020, see their Appendix E).We set ω/n = 1 at the start of each simulation, and assume that ω/n ≃ 1 throughout.Thus, we take sgn(2ω − 3n) = −1, sgn(2ω − n) = +1 and sgn(n) = +1.Proxima b has zero obliquity in each of our simulations, and so we discard terms related to obliquity tides.

Figure 1 .
Figure 1.Panel (a) shows a stability map of Proxima d and b, with both eccentricities varying between 0 and 0.9.Blue cells indicate periodic evolution according to the MEGNO criterion, red cells indicate one planet was ejected, and white cells indicate chaotic (but not necessarily unstable) evolution.The solid black curve marks the Hill stability boundary in this plane, with Hill stable cases lying below the curve.The dashed-dotted black line indicates the limit of crossing orbits, i.e., above this line the apocenter distance of planet d is initially greater than the pericenter distance of planet b.The dotted line indicates the stability limit of Hadden & Lithwick (2018).Yellow stars indicate the simulations shown in detail in Figure 2. Green dots indicate the locations of the secular dynamics simulations detailed in Section 2.5.Grey triangles indicate the locations of the longer (1 Myr) integrations described in Section 3.1.Panel (b) shows a MEGNO map of the same parameter space.Most points in this map are similarly colored, as these simulations concluded with ⟨Y ⟩ ≥ 10 or with an ejection, in which case an arbitrarily high MEGNO value was assigned.

Figure 2 .
Figure 2. Example evolutions of Proxima d and b.Periodic evolution originating from the initial conditions e d = e b = 0.05 is shown in panel (a).An evolution originating from e d = 0.3 and e b = 0.2 is shown in panel (b).The time-averaged MEGNO for this case diverges significantly from 2, signifying that the motion is chaotic.At least 3 distinct frequencies are seen to be present in both planets' eccentricities, which is greater than the number predicted by secular theory.An unstable case originating from e d = 0.45 and e b = 0.3 is shown in panel (c).The ejection of planet d from the system is evident at t = 222 yr, where e d > 1.

Figure 3 .
Figure 3. MEGNO values of the system over all possible initial mutual inclinations.The dashed black line indicates ⟨Y ⟩ = 2.Blue dots are cases with periodic or quasi-periodic motion, grey dots are cases with chaotic motion, and the red ticks indicate values of Ψ for which the system is unstable (and therefore any MEGNO value is meaningless).All cases but one between 50 and 140 • resulted in chaotic motion.165 simulations resulted in ejections, all with mutual inclinations between 95 and 142 • .

Figure 4 .
Figure 4. Similar to Figure 1, but with the initial inclinations varied instead of eccentricities.The black lines in panel (a) indicate the Hill stability boundary.Yellow stars indicate the simulations presented in detail in Figure 5, and the green dots point out two combinations of inclinations used in the secular dynamics simulations.

Figure 7 .
Figure 7. Secular orbital evolution of Proxima d in each simulated case is displayed in panels (a) and (b), and Proxima b in panels (c) and (d), with the tidal heat flux through the planet b's surface shown in panel (e).Three simulations end with a tidal heat flux of ∼ 0.1 W m −2 , two end with a near-Earth value of ∼ 0.01 W m −2 , and one with a far lower value.The present tidal heat flux at Io's and Earth's surfaces are indicated for reference(Moore et al. 2007;Munk & Wunsch 1998).

Table 1 .
Best-fit parameters for the Proxima Centauri system.

Table 2 .
Initial conditions for the REBOUND simulations.

Table 3 .
Distributions from which the initial conditions of our VPLanet simulations are sampled.