Predicting multiple planet stability and habitable zone companions in the TESS era

TL;DR

This paper introduces a rapid, predictive method to assess the dynamical stability of multi-planet systems and their habitable zones, aiding the search for stable Earth-like companions in the TESS era.

Contribution

The authors develop a novel numerical tool based on test particle simulations to efficiently evaluate stability regions and habitable zone companions in exoplanet systems.

Findings

Test particles survive in stable, unperturbed or resonant orbits.

Stable regions produce a unique (a,e) signature.

The method can assess stability without exhaustive simulations.

Abstract

We present an approach that is able to both rapidly assess the dynamical stability of multiple planet systems, and determine whether an exoplanet system would be capable of hosting a dynamically stable Earth-mass companion in its habitable zone. We conduct a suite of numerical simulations using a swarm of massless test particles in the vicinity of the orbit of a massive planet, in order to develop a predictive tool which can be used to achieve these desired outcomes. In this work, we outline both the numerical methods we used to develop the tool, and demonstrate its use. We find that the test particles survive in systems either because they are unperturbed due to being so far removed from the massive planet, or due to being trapped in stable mean motion resonant orbits with the massive planet. The resulting unexcited test particle swarm produces a unique signature in (a,e) space that represents the stable regions within the system. We are able to scale and translate this stability signature, and combine several together in order to conservatively assess the dynamical stability of newly discovered multiple planet systems. We also assess the stability of a system's habitable zone and determine whether an Earth-mass companion could remain on a stable orbit, without the need for exhaustive numerical simulations.

Tables

Table 1. Table 1 : The orbital parameters of the massive planet used in each set of simulations.

Parameter	Set 1	Set 2	Set 3
$M$ ($M_{\oplus }$)	$1,2,4,\dots ,1024$	32.0	32.0
$a$ (au)	1.0	1.0	1.0
$e$	0.0	$0,0.05,0.1,\dots ,0.3$	0.0
$i$ ($°$)	0.0	0.0	$0,2.5,5,\dots ,10$
$\Omega$ ($°$)	0.0	0.0	0.0
$\omega$ ($°$)	0.0	0.0	0.0
$M$ ($°$)	0.0	0.0	0.0

Table 2. Table 2 : Multiple exoplanet systems stability signature predictions compared with the stability of the system determined from detailed numerical simulations of Agnew et al. ( 2018b ) . This table is intended to demonstrate the agreement between our approach and numerical simulations, but should not be seen as a definitive assessment of those system’s stability.

System	Numerical	Signature	Agreement?
	Result	Prediction
24 Sex	Unstable	Unstable	✓
47 UMa	Stable	Potentially Stable	-
BD-08 2823	Stable	Stable	✓
HD 10180	Stable	Potentially Stable	-
HD 108874	Stable	Potentially Stable	-
HD 113538	Stable	Potentially Stable	-
HD 134987	Stable	Potentially Stable	-
HD 141399	Stable	Potentially Stable	-
HD 142	Stable	Potentially Stable	-
HD 159868	Stable	Potentially Stable	-
HD 1605	Stable	Potentially Stable	-
HD 160691	Unstable	Potentially Stable	-
HD 187123	Stable	Potentially Stable	-
HD 200964	Unstable	Unstable	✓
HD 219134	Stable	Stable	✓
HD 33844	Unstable	Unstable	✓
HD 47186	Stable	Stable	✓
HD 4732	Stable	Potentially Stable	-
HD 5319	Unstable	Unstable	✓
HD 60532	Stable	Potentially Stable	-
HD 9446	Stable	Potentially Stable	-
HIP 65407	Stable	Stable	✓
HIP 67851	Stable	Potentially Stable	-
TYC 1422-614-1	Stable	Potentially Stable	-
XO-2 S	Stable	Potentially Stable	-

Table 3. Table 3 : Summary of stellar mass, planetary mass and planetary semi-major axes of the planet/s in the HD 5319 and Kepler 16 systems as they are currently known.

	HD 5319	b	c
Mass	$1.51M_{\odot }$	$1.76M_{\mathrm{Jup}}$	$1.15M_{\mathrm{Jup}}$
Normalised Mass	$1M_{\odot }$	$1.166M_{\mathrm{Jup}}$	$0.762M_{\mathrm{Jup}}$
Semi-major axis [au]	-	$1.6697$	$2.071$
	Kepler-16	b
Mass	$0.69M_{\odot }$	$0.333M_{\mathrm{Jup}}$	-
Normalised Mass	$1M_{\odot }$	$0.483M_{\mathrm{Jup}}$	-
Semi-major axis [au]	-	$0.7048$	-

Equations

$$a_{\text{pl}}=\sqrt[3]{\frac{G{M}_{\star }{T}_{\text{pl}}^2}{4{\pi }^2}},$$

$$a\propto T^{2/3}.$$

$$a\propto M_{\star }^{1/3}.$$

$$\{(a,f_{\mathrm{norm}}(a)):a\in n^{1/3}A\},$$

$$f_{\mathrm{norm}}(n^{1/3}a)=f(a).$$

$$k$$

$$=a_{\mathrm{pl}}/n^{1/3}.$$

$$\{(a,f_{\mathrm{planet}}(a))$$

$$f_{\mathrm{planet}}(a_{\mathrm{pl}}a)=f(a).$$

$$\frac{1}{1.51}(1.76+1.15){\mathrm{M}}_{\mathrm{Jup}}=1.927{\mathrm{M}}_{\mathrm{Jup}}=612.8{\mathrm{M}}_{\oplus }$$

$$a\in [0.215,4.642]$$

$$a\in a_{\mathrm{pl}}[0.215,4.642]$$

$$a\in 1.6697\cdot [0.215,4.642]$$

$$⟹a\in [0.359,7.751]$$

$$a\in 2.071\cdot [0.215,4.642]$$

$$⟹a\in [0.445,9.614]$$

$$\frac{1}{0.69}0.333{\mathrm{M}}_{\mathrm{Jup}}=0.483{\mathrm{M}}_{\mathrm{Jup}}=153.4{\mathrm{M}}_{\oplus }$$

$$a\in [0.215,4.642]$$

$$a\in a_{\mathrm{pl}}[0.215,4.642]$$

$$a\in 0.7048\cdot [0.215,4.642]$$

$$⟹a\in [0.152,3.272]$$

Full text

††pubyear: 2017††pagerange: Predicting multiple planet stability and habitable zone companions in the TESS era–C

1 Introduction

The search for potentially habitable worlds is an area of immense interest to the exoplanetary science community. Since the first discoveries of exoplanets orbiting main sequence stars (Campbell et al., 1988; Latham et al., 1989; Mayor & Queloz, 1995), planet search surveys have endeavoured to discover the degree to which the Solar system is unique, and to understand how common (or rare) are planets like the Earth (e.g. Howard et al., 2010a; Wittenmyer et al., 2011). The launch of the Kepler spacecraft in 2009 led to a great explosion in the number of known exoplanets (e.g. Borucki et al., 2010; Sullivan et al., 2015; Morton et al., 2016; Dressing et al., 2017). Kepler carried out the first census of the Exoplanet era, discovering more than two thousand planets111As of 13 September 2018 (Kepler and K2 mission site, https://www.nasa.gov/mission_pages/kepler/main/index.html).. Kepler’s results offer the first insight into the true ubiquity of planets, and the frequency with which Earth-size planets within the Habitable Zone (HZ) can be found orbiting Sun-like stars (Catanzarite & Shao, 2011; Petigura et al., 2013; Foreman-Mackey et al., 2014). The Transiting Exoplanet Survey Satellite (TESS) seeks to continue this trend in exoplanet discoveries (Ricker et al., 2014; Sullivan et al., 2015; Barclay et al., 2018), concentrating on planet detection around bright stars that are more amenable to follow-up spectroscopy (Kempton et al., 2018).

In addition to the rapid increase in the known exoplanet population, the generational improvement of instruments being used for detection, confirmation, and observational follow-up, has recently allowed for planets to be detected with masses comparable to that of the Earth, albeit on orbits that place them far closer to their host stars than the distance between the Earth and the Sun (e.g. Vogt et al., 2015; Wright et al., 2016; Anglada-Escudé et al., 2016; Gillon et al., 2017). Whilst such Earth-mass planets generate larger, more easily detectable radial velocity signals, current and near-future instruments (e.g. the ESPRESSO and CODEX spectrographs) seek to detect planets inducing Doppler wobbles as low as $0.1~{}\textrm{ m s}^{-1}$ and $~{}0.01~{}\textrm{m s}^{-1}$ respectively (Pasquini et al., 2010; Pepe et al., 2014; González Hernández et al., 2017). At such small detection limits, those spectrographs may well be able to detect Earth-mass planets that orbit Sun-like stars at a distance that would place them within the star’s HZ (Agnew et al., 2017; Agnew et al., 2018a, b). With the advent of these new high-precision spectrographs, as well as the launch of the next generation of space telescopes (such as the James Webb Space Telescope, JWST), it will be possible to perform observational follow-up in far greater detail. As a result, it is timely to consider methods by which we might prioritise the observations needed in order to detect Earth-size planets (Horner & Jones, 2010).

Whilst the detection of a Solar system analogue would be considered the holy grail in the search for an exoplanet that truly mimics Earth, finding something that resembles the Solar system itself beyond a handful of similarities has so far proven to be particularly challenging (Boisse et al., 2012; Wittenmyer et al., 2014b, 2016; Kipping et al., 2016; Rowan et al., 2016; Agnew et al., 2018a). By relaxing the criteria of our search to instead look for multiple planet systems where at least one planet is comparable in mass to Earth and resides in its host star’s HZ (Kasting et al., 1993; Kopparapu et al., 2013, 2014; Kane et al., 2016), we will still yield exoplanetary systems that share several similarities with our own system, and are still candidates for further study from the perspective of planetary habitability. Additionally, we consider the notion that the observational bias inherent to several detection methods (Wittenmyer et al., 2011; Dumusque et al., 2012) can be interpreted to suggest that systems with massive, giant planets may also coexist with smaller, rocky exoplanets that so far have been undetectable due to detection limits (Agnew et al., 2017). Indeed, given the challenges involved in finding habitable exo-Earths, it might be the case that such planets exist within known exoplanetary systems, lurking below our current threshold for detectability. Future studies of those systems might allow such planets to be discovered, if they exist – which serves as an additional motivation for the development of a method by which we can prioritise systems as targets for the search for Earth-like planets: systems that host dynamically stable HZs.

The standard method used to determine whether a system could host unseen planetary companions is by assessing the system’s overall dynamical stability. This can be done analytically (e.g. Giuppone et al., 2013; Laskar & Petit, 2017), or numerically (e.g. Raymond & Barnes, 2005; Jones et al., 2005; Rivera & Haghighipour, 2007; Jones & Sleep, 2010; Wittenmyer et al., 2013; Agnew et al., 2017). Such studies have resulted in the development of several methods that allow exoplanetary systems to be assessed by planetary architectures as they are observed today (Giuppone et al., 2013; Carrera et al., 2016; Matsumura et al., 2016; Agnew et al., 2018a). It has been frequently demonstrated that a variety of resonant mechanisms are often integral in determining the dynamical stability of a system (e.g. Wittenmyer et al., 2012; Gallardo, 2014; Kane, 2015; Gallardo et al., 2016; Mills et al., 2016; Luger et al., 2017; Delisle, 2017; Agnew et al., 2017; Agnew et al., 2018b). Typically, such studies examine a single exoplanetary system, and study it in some depth to determine whether it is truly dynamically stable, and whether other planets could lurk undetected within it. Such efforts are extremely computationally intensive, which means that few systems can readily be studied in such detail.

Here, we consider whether it is instead possible to use detailed $n$-body simulations to build a more general predictive tool that could allow researchers to quickly assess the potential stability (or instability) of any given system without the need for their own suite of numerical simulations. This would identify systems for which further observational investigation might be needed. In doing so, we develop a tool by which one can: 1) rapidly and conservatively assess the dynamical stability of a newly discovered multiple planet system to identify those where further observation is required to better constrain orbital parameters, and 2) identify the regions in a given exoplanetary system that could host as-yet undiscovered planets, based on their dynamical interaction with the known planets in the system. We present a suite of simulations that have identified the stable regions around generic planetary systems that can be used to examine the stability of specific systems. Ultimately, we wish to demonstrate the use cases of our predictive tool to rapidly assess new systems found with TESS.

In section 2 we introduce the use cases that we have in mind when developing our predictive tool. We outline the fundamental approach of our method in section 3, and demonstrate how the general stability signatures we compute can be normalised and translated to fit any single planetary system discovered. In section 4 we show how the mass, eccentricity and inclination of a massive body influences its stability signature. We follow this with various examples of how our method can be used to infer system stability or to determine where stable resonant HZ companions may exist in section 5, and summarise our findings in section 6.

2 Dynamical Predictions

The predictive tool we present was developed to be applied to exoplanetary systems discovered by TESS with two goals in mind: 1) to determine the dynamical feasibility of newly discovered multiple planet systems, and 2) to predict the regions of stability in TESS systems where another unseen planet may exist.

The first use case is to provide a tool that may be incorporated into the Exoplanet Follow-up Observing Program for TESS (ExoFOP-TESS)222https://heasarc.gsfc.nasa.gov/docs/tess/followup.html. When a multiple planet system is discovered, we can utilise our tool to dynamically assess the stability of the system given the inferred orbital parameters. While our approach will miss more complex destabilising behaviour (such as the influence of secular resonant interactions as demonstrated by Agnew et al., 2018b), it can conservatively assess the dynamical stability of a system “on the fly”. Demonstrating instability using a conservative approach would suggest further observations are required to better constrain the orbital parameters of the planets, or to reassess the number of planets in the system in the cases of potential eccentricity harmonics and aliasing (Anglada-Escudé et al., 2010; Anglada-Escudé & Dawson, 2010; Wittenmyer et al., 2013).

The second use case is to assist in the search for potentially habitable Earth-size planets. Using $n$-body simulations, we can make dynamical predictions of regions in the ($a$,$e$) parameter space where additional planets will definitely be unstable. This formed the core of our previous work in this series (Agnew et al., 2017; Agnew et al., 2018a, b), where we were able to identify systems that could potentially host dynamically stable HZ planets. Here, we extend that work to develop a stability mapping process that will enable future studies to quickly identify the regions in newly discovered exoplanetary systems where planets can definitely be ruled out, on dynamical grounds, as well as those regions where additional planets could only exist under very specific conditions (such as when trapped in a mutual mean-motion resonance with another planet, e.g. Howard et al., 2010b; Robertson et al., 2012; Wittenmyer et al., 2014a).

To develop such a tool, we performed a large number of detailed $n$-body simulations to determine how the unstable regions centred on a given planet’s orbit are influenced by its mass, orbital eccentricity and inclination in order to produce a scale free template that can be applied to any system.

3 Method

We have developed our predictive tool by conducting thorough, high resolution test particle (TP) simulations. In each simulation, a swarm of TPs are distributed randomly throughout a region in ($a$,$e$) space around the orbit of a massive planet. At the end of the simulation, the distribution of surviving TPs reveals evacuated regions (which correspond to areas of instability), and regions where TPs remain unperturbed, and so remain on similar orbits to those held at the beginning of the integrations. Since our simulations are primarily focused on facilitating the search for habitable worlds, we chose to enforce a maximum initial eccentricity on the orbits of TPs of 0.3, based on studies that suggest that the habitability of a given Earth-like world would be significantly reduced for higher eccentricities (Williams & Pollard, 2002; Jones et al., 2005).

The goal of these simulations is to determine the stable, unperturbed regions around a massive body, which we refer to as that body’s “stability signature”. More specifically, we put forward two definitions: the optimistic stability signature being the curve that bounds the maximum eccentricity of the unexcited TPs, and the conservative stability signature that bounds the minimum eccentricity of the excited TPs. The key reasons for considering the signature as a curve rather than a 2-dimensional area in ($a$,$e$) space are: 1) in the simplest case (the massive body moving on a circular orbit) only those TPs with sufficiently eccentric orbits have apsides that enter within the region of instability nearby to the massive planet (generally some multiple of its Hill radius), and so the regions below the curve (i.e. the TPs with less eccentric orbits) will be stable, and 2) multiple curves can be combined and presented in a “look-up map” within which we can interpolate to find the curves for planetary masses we do not simulate explicitly. This is explained in greater detail in section 4.

By determining the stability signatures for a range of planet masses, the signatures can be used as a scale free template that can be applied to any system in order to assess its dynamical stability without the need to run numerical simulations. Our proposed exploration of the planet mass parameter space and the resulting stability signatures will also enable researchers to predict the regions of a given exoplanetary system where additional planets are dynamically feasible, and those where no planets are likely to be found. We anticipate that this method will provide a useful filtering tool by which systems can be examined to quickly predict whether they are dynamicaly feasible, and which might be able to host a potentially habitable exoplanet. We present here the general simulations we carried out from which we determined the stability signatures and how the signatures from these simulations can be used to create the stability template for any specific system.

In this work, we consider the stability of additional bodies moving on co-planar orbits in circular, single planet systems. While it is highly unlikely that a planet will have a perfectly circular orbit, our intent is for this to be a conservative assessment of the stability of the system. Our approach can be expanded and refined to also consider eccentric and inclined orbits, which we will pursue in future work.

3.1 Numerical approach

Two observables when detecting an exoplanet are stellar mass, $M_{\star}$ (obtained indirectly from the luminosity of the star, $L_{\star}$), and the orbital period of a planet, $T_{\textrm{pl}}$ (obtained directly). We can use these to calculate the semi-major axis of the planet by

[TABLE]

where $G$ is the universal Gravitational constant, and hence we can infer the distance of a planet from the observed orbital period. For any fixed $M_{\star}$ value, the semi-major axis will scale with orbital period according to the power law

[TABLE]

Kopparapu et al. (2014) put forward a method to calculate the HZ of a system using $M_{\star}$ and the planetary mass, $M_{\textrm{pl}}$. In general, a minimum mass for $M_{\textrm{pl}}$ is determined via the radial velocity method. Taken in concert with the above, this means that from the observed parameters $L_{\star}$ and $T_{\textrm{pl}}$, we can infer $M_{\star}$ and $a_{\textrm{pl}}$, measure $M_{\textrm{pl}}$, and hence compute the HZ around the star.

By numerically simulating a range of mass ratios ($\mu=M_{\textrm{pl}}/M_{\star}$) while keeping constant the mass of the star ($M_{\star}=\mathrm{M}_{\odot}$) and semi-major axis of the planet ($a_{\mathrm{pl}}=1$ au), we produce the stability signature for each of the simulated masses and can consolidate them to make a “look-up map”. To apply our general stability signatures to a newly discovered system, one would first determine the mass ratio between the planet and the star in that system, and use that mass ratio to interpolate between our simulated stability signatures in our look-up map. This yields the signature around a planet of any mass $M_{\textrm{pl}}$. Once this step is complete, one then translates the signature obtained from the nominal $a_{\mathrm{pl}}=1$ au to a location closer to, or farther from, the host star to match the discovered planet. We can then compare the stability signatures of all planets in multiple planet systems with one another to conservatively assess dynamical stability, as well as determine which systems can host hypothetical exo-Earths within their HZ, without the need to run numerical simulations.

This approach limits the stability constraint to be multiples of the orbital period rather than number of years. Our numerical simulations are for a planet orbiting at $a_{\textrm{pl}}=1$ au in a $M_{\star}=1~{}\mathrm{M}_{\odot}$ system (i.e. $T_{\textrm{pl}}=1$ year). As we simulate our systems for $10^{7}$ years (i.e. $10^{7}$ orbits), if we were to use our stability signature for a planet that orbits at a semi-major axis with an orbital period of $0.1$ years (i.e. $1/10^{\textrm{th}}$ of the simulated orbital period) we can only refer to the planet as being stable for $10^{6}$ years (i.e. $10^{7}$ orbits).

3.2 Normalisation of system

As our general simulations were performed for a range of mass ratios, but for a fixed stellar mass, normalisation is important so that we are still able to determine the stability signature around a newly discovered system regardless of its stellar mass.

As per our definition of the stability signature being the curve that bounds the maximum eccentricity of the unexcited TPs, let us consider a hypothetical function that maps semi-major axis, $a$, to this curve, $\{(a,f(a)):a\in A\}$, where $A$ is the domain over which we simulated. Our simulations all use a stellar mass of $M_{\star}=1~{}\mathrm{M}_{\odot}$, and a planet orbiting at $a_{\mathrm{pl}}=1$ au (i.e. $T_{\mathrm{pl}}=1$ yr). To normalise a system, we first determine what the semi-major axis is that matches an orbital period of one year for a system with a different stellar mass. From equation 1, the semi-major axis will scale with stellar mass according to the power law

[TABLE]

For a given system discovered with stellar mass $M_{\star}=n\mathrm{M}_{\odot}$, we can scale the masses $M_{\star}$ and $M_{\mathrm{pl}}$ by $1/n$, and compute that the one year orbital period will occur at $a_{\mathrm{1year}}=n^{1/3}$ au. We demonstrate in Figure 1 how the stable signatures are shown to be identical when the mass ratio, $\mu$, and orbital period, $T_{\textrm{pl}}$, remain constant. Thus, when an exoplanet of mass $M_{\mathrm{pl}}$ is discovered orbiting a star of mass $M_{\star}$, it is possible to first normalise the system by scaling the planetary and stellar masses while retaining a constant mass ratio.

Once normalised, we can then interpolate between the planetary masses we simulated in order to obtain the stability signature for any planetary mass (within the upper and lower $\mu$ we simulate). In terms of function notation, a new function for the normalised stability signature $f_{\mathrm{norm}}$ will have the scaled domain of $f$, while still mapping to the unscaled stability signature values $f(a)$, that is,

[TABLE]

where $f_{\mathrm{norm}}(a)$ is the function mapping the scaled domain to the original signature curve, given by

[TABLE]

We must then determine how the domain will vary for a function that is translated from the semi-major axis where an orbital period of one year occurs, $a_{\mathrm{1year}}$, to the semi-major axis of the discovered exoplanet, $a_{\mathrm{pl}}$.

3.3 Translatability of signature

As our general simulations are run for a planet at a fixed semi-major axis of $a_{\mathrm{pl}}=1~{}\mathrm{au}$, translatability is important so that we are still able to determine the stability signature around a newly discovered system regardless of the semi-major axis of the planet. We translate the normalised signature by simply taking the ratio between the semi-major axis of the planet ($a_{\mathrm{pl}}$) and the one year semi-major axis ($a_{\mathrm{1year}}$). We express this ratio as

[TABLE]

The discovered exoplanet’s stability signature at its semi-major axis is then obtained by multiplying the normalised signature by $k$. In terms of function notation, a new function for the planet’s stability signature $f_{\mathrm{planet}}$ will have the translated domain of $f_{\mathrm{norm}}$, while still mapping to the unscaled stability signature values $f(a)$, that is,

[TABLE]

where $f_{\mathrm{planet}}(a)$ is the function mapping the scaled and translated domain to the original signature curve, that is,

[TABLE]

Thus, from the observed parameters $L_{\star}$ and $T_{\textrm{pl}}$, we can infer $M_{\star}$ and $a_{\textrm{pl}}$, and from radial velocity measurements we can determine $M_{\textrm{pl}}$. With these parameters we are able to: 1) normalise the system by scaling $M_{\star}$ and $M_{\textrm{pl}}$, 2) interpolate between the planetary masses in our look-up map in order to obtain the normalised signature, and 3) translate the normalised signature to the semi-major axis of the detected planet to find its stability signature without the need for additional numerical simulations. Being able to obtain the stability signature of any discovered planetary system allows us to: 1) combine the stability signatures of all planets within a multiple planet system in order to rapidly and conservatively assess its dynamical stability, and 2) determine the stable, unperturbed regions around an exoplanet in order to constrain where a habitable terrestrial planet could exist. These applications will be discussed in depth in section 5.

3.4 Simulations

All our simulations model a two-body, star–planet system within which we randomly scatter $10^{5}$ massless TPs in order to test the stability of a hypothetical third body. The motion of two massive bodies and a third massless body constitutes the Restricted Three-Body Problem. The $10^{5}$ TPs represent $10^{5}$ possible orbital parameter configurations for this third body. These TPs are scattered in order to find the stable regions within the ($a$,$e$) parameter space within the vicinity of the orbit of a massive body. The central star for our simulations has mass $M_{\star}=1~{}\mathrm{M}_{\odot}$, and the semi-major axis of the planet is $a_{\mathrm{pl}}=1~{}\mathrm{au}$. We then carry out three distinct suites of simulations, in which we vary the planetary mass ($m_{\mathrm{pl}}$), eccentricity ($e_{\mathrm{pl}}$) or inclination ($i_{\mathrm{pl}}$). The orbital parameters of the planet for each suite of simulations is summarised in Table 1. The orbital parameters for the swarm of TPs are randomly generated within a fixed range for all three sets of simulations, and will be discussed in detail after first outlining how the planetary parameters are varied in each simulation.

In our first set of simulations, we explore the effect of planetary mass, $m_{\mathrm{pl}}$. Since the mass of the host star is kept constant, this allows us to examine the effect of the mass ratio $\mu=M_{\mathrm{pl}}/M_{\star}$ on the resulting stability signatures. In these simulations, we set the orbital parameters of the planet to those values shown for Set 1 in Table 1. The mass of the planet is then varied sampling the range $M_{\mathrm{pl}}=1$ to $1024\ \textrm{M}_{\oplus}$, with the mass increased incrementally by factors of 2 (i.e. $1,2,4,8,...$). We span these masses so as to appropriately cover the expected mass distribution of planets found with TESS (Sullivan et al., 2015).

In our second set of simulations we explore the effect of eccentricity. In these simulations we fix the planetary mass to $M_{\mathrm{pl}}=32\ \textrm{M}_{\oplus}$ as this falls within the dominant mass range of expected TESS findings (Sullivan et al., 2015). In these simulations, we set the orbital parameters of the planet to those values shown for Set 2 in Table 1. The eccentricity of the planet is then varied from $e_{\mathrm{pl}}=0.0$ to $0.3$ in steps of $0.05$ for each simulation. In our final set of simulations, we explore the effect of inclination. In these simulations we again fix the planetary mass to $M_{\mathrm{pl}}=32\ \textrm{M}_{\oplus}$. and set the orbital parameters of the planet to those values shown for Set 3 in Table 1. The inclination of the planet is then varied from $i_{\mathrm{pl}}=0.0\degree$ to $10.0\degree$ in steps of $2.5\degree$ for each simulation.

For all of our simulations we randomly scatter $10^{5}$ massless TPs throughout the system. As we are interested in obtaining the stability signatures in the vicinity of the orbit of a massive body, the TPs are distributed between orbits with periods in $10:1$ and $1:10$ commensurability with the planet (which for $a_{\mathrm{pl}}=1~{}\mathrm{au}$ is given by $0.215~{}\mathrm{au}\lesssim a\lesssim 4.642~{}\mathrm{au}$), and with eccentricities $0.0\leq e\leq 0.3$. It should also be noted that, by varying the values of $\omega_{\mathrm{tp}}$ and $M_{\mathrm{tp}}$ across the full [math] to $360\degree$ range, this allows for a fixed value of $\omega_{\mathrm{pl}}=0.0\degree$ and $\Omega_{\mathrm{pl}}=0.0\degree$ for sets 2 and 3 respectively to still cover the relevant parameter space.

3.5 Analysis

One of the key goals of this work was the development of a tool that can be used to assess the dynamical stability of a system without the need to run numerical simulations. We achieve this by running n-Body simulations with a swarm of massless TPs to produce a set of scalable templates, which reveals the unperturbed regions around a planet, which we refer to as that planet’s “stability signature”. More specifically, we determine two stability signatures: one that optimistically provides an upper bound above which TPs are shown to be excited, and a second that more conservatively provides a lower bound below which no TPs were shown to be excited.

For each simulation, we determine the system’s stability signatures based on the excitation of the TPs perturbed by the planet over a period of $10^{7}$ years. We consider a TP to be excited by the planet if its semi-major axis changes relative to its initial position by $\Delta a_{\mathrm{excited}}=(a_{\textrm{final}}-a_{\textrm{initial}})/a_{\textrm{initial}}>0.01$ over the course of the entire simulation. This is because ultimately we wish to use the signatures to predict the locations most likely to host additional planets in the system – particularly Earth-mass bodies. Since our simulations use massless TPs, they do not take into account mutual gravitational interactions between the planet and any companion – represented by the TPs – orbit. In order to use TPs to predict the stability of massive bodies with any confidence (such as demonstrated by Agnew et al., 2017; Agnew et al., 2018b), we only need to consider those TPs that are not excited (in other words, the gravitational influence of the known exoplanet is not destabilising it). The gravitational influence of the less massive Earth-mass planet is also unlikely to perturb the known exoplanet in this scenario. In the cases of our first and third sets of simulations, where the massive planet moves on a circular orbit, the Jacobi integral is constant as the system is a circular restricted three-body problem. The dynamics in this case differs from the restricted three-body problem that we see with eccentric orbits, and so the results of Set 2 will also differ significantly due to the difference in the dynamical evolution resulting from the additional eccentricity of the planet itself.

Thus, by considering the relative change in semi-major axis of TPs, and ignoring those that are excited by the gravitational influence of the planet (i.e. $\Delta a_{\mathrm{excited}}>0.01$), we are left with the most stable, unexcited bodies. From this, we can extract the optimistic stability signature by placing the unexcited TPs into bins in semi-major axis, and taking the maximum eccentricty of the TPs in each bin. In contrast, we obtain the conservative stability signature, below which we see no TP excitation, by placing the excited TPs into bins in semi-major axis, and taking the minimum eccentricity of the TPs in each bin. In both scenarios, we disregard outliers. We do this for two reasons. Firstly, tiny regions of stability seem to be unlikely places for planets to form and survive in all but the youngest stellar systems, and so such regions typically offer little prospect of predicting potential stable planet locations. In addition, we are interested in regions of stability that are represented by stable swarms of TPs, rather than individual outliers, as stable resonances are not infinitely thin, but have a measurable width in semi-major axis space. In this work, we ignore surviving TPs that are alone in a given bin. For bins that contain multiple surviving TPs, we consider the cumulative distribution function (CDF) of the eccentricities of the TPs. We ignore any TPs that have a final eccentricity in the top $2.5\%$ of the CDF within each bin. The converse process is true for ignoring outliers when determining the conservative stability signature, considering the excited TPs instead of the surviving TPs.

Plotting a curve connecting the eccentricity of the most eccentric, unexcited TP in each semi-major axis bin – and for which outliers have been removed – yields a signature like that shown by the yellow curve in the upper plots of Figure 3. Similarly, we can plot the curve of the least eccentric, excited TPs as shown by the green curve (middle plots). Figure 3(a) shows that at low masses ($\mu={8\textrm{M}_{\oplus}}/\textrm{M}_{\odot}$), other than the very obvious region that is cleared in the planet’s immediate vicinity, the optimistic signature is very difficult to extract, and as a result, the accuracy with which this can be achieved is somewhat limited. However, the curve of the excited TPs is much more clearly defined, and so this can still be acquired easily. Figure 3(b) shows that for more massive planets ($\mu={64\textrm{M}_{\oplus}}/\textrm{M}_{\odot}$), the outlined method for obtaining the optimistic stability signatures is much cleaner as the signature is more distinct. The optimistic stability signatures (the yellow curves in Figure 3) are thus the curve in ($a$,$e$) space, above which are the perturbed regions near a given planet. The conservative stability signatures (the green curves) are the curve in ($a$,$e$) space below which are the stable, unperturbed regions near a given planet. These can be combined as shown in the tri-colour figures (bottom plots) to yield a classification scheme where planets can be deemed stable, unstable, or potentially stable in between these two extremes where we have evidence of TPs being both excited and unexcited.

Using the approach outlined above, the stability signatures were extracted for all of our simulations and compiled to produce the desired look-up maps, and to compare how different parameters might impact dynamical stability of the systems. The resulting look-up map forms the basis of the predictive tool we developed in assessing system stability without further numerical simulations.

4 Results and Discussion

4.1 Effect of mass

The influence of the mass of the known planet on the stability in a given system is the main focus of this work, and as such is covered in the most detail. We present the compiled stability signatures and look-up maps in Figure 5.

Figure 4(a) shows the stability signature for all the planetary masses we simulated in Set 1 in Table 1. The stability signatures were extracted as outlined above, and specifically the signatures for $\mu={8\textrm{M}_{\oplus}}/\textrm{M}_{\odot}$ and $\mu={64\textrm{M}_{\oplus}}/\textrm{M}_{\odot}$ can be seen in Figure 3. The $x$-axis shows the semi-major axes between orbits with periods in $10:1$ and $1:10$ commensurability with the planet (which for $a_{\mathrm{pl}}=1~{}\mathrm{au}$ is given by $0.215~{}\mathrm{au}\lesssim a\lesssim 4.642~{}\mathrm{au}$), while the $y$-axis shows the eccentricity values of the stability signatures in Figure 4(a), and indicates the mass ratio of the planet in each simulation in Figures 4(b) and 5(a). As the maximum eccentricity we tested with our simulations was $0.3$, this means that the maximum value in these signatures correspond with very stable, unexcited regions for a body up to an eccentricity of $0.3$. While the maximum values could mean bodies with eccentricities higher than $0.3$ may also be stable at these locations, as we did not explore higher than $0.3$ we can only predict up to this value. We also overlay various other boundaries for the onset of chaos. These include simpler approximations, such as the boundaries of $1$, $3$ and $5$ Hill radii from the massive planet, as well as more developed analytic definitions. There have been several studies into the onset of chaos (e.g. Mustill & Wyatt, 2012; Giuppone et al., 2013; Deck et al., 2013; Petit et al., 2017; Hadden & Lithwick, 2018), and here we plot those results presented by Mustill & Wyatt (2012), Giuppone et al. (2013) and Deck et al. (2013).

We emphasise that for low mass planets the optimistic stability signatures are difficult to extract accurately without making significant assumptions or changes to the extraction methodology. This was demonstrated for $\mu={8\textrm{M}_{\oplus}}/\textrm{M}_{\odot}$ in Figure 3(a), and can be seen for other masses $\mu\leq{8\textrm{M}_{\oplus}}/\textrm{M}_{\odot}$ in Figure 4(a). Despite this, the conservative signature for such low mass scenarios can still be seen to be demonstrably stable for all but the regions located nearest to the planet, but this should be kept in mind for any systems with $\mu\leq{8\textrm{M}_{\oplus}}/\textrm{M}_{\odot}$.

Figure 4(a) shows that for systems with $\mu>{8\textrm{M}_{\oplus}}/\textrm{M}_{\odot}$, the optimistic stability signatures are much less noisy and hence more distinctly defined. The semi-major axes that correspond with stabilising mean motion resonances are also clearly visible (shown by the dashed magenta lines), as is the manner in which the widths of those resonances can vary with planetary mass. In contrast, Figure 4(a) shows that the conservative signals are much more distinctly defined even for low mass ratios, and so are useful for identifying stability.

Figures 4(b) and 5(a) shows the look-up maps we developed by representing the eccentricity values of the stability signature for each $\mu$ as shaded areas. The dark regions represent strongly stable regions, specifically where the maximum eccentricity of an unexcited TP is $0.3$ (the maximum value we used in our simulations). Conversely, the lighter regions represent the unstable regions, where the maximum eccentricity of an unexcited TP is zero or near zero, meaning that region in ($a$,$e$) space has been completely cleared of TPs. The stabilising resonances are again evident in Figure 4(b) (the vertical dark streaks embedded in the light white cloud, highlighted by the dashed magenta lines), as is the strong influence that the mass ratio has on the reach of the unstable region. We see similar features in Figure 5(a), but in this case the mean motion resonances can be seen to have a destabilising effect on the TPs, specifically the 2:1 and 1:2 MMRs. In general, the various analytic definitions for the onset of chaos bounds or straddles the unstable regions we present to some extent in both the optimistic and conservative cases, suggesting there is use to both the optimistic and conservative definitions depending on the desired application. It also highlights that a more robust definition of excitation utilising the derivations in these works may yield better, more precisely defined signatures. Generally, these analytic criteria derive the boundary of the onset of chaos to be related to the relative difference in semi-major axis between the two bodies. As such, monitoring this throughout each simulation may yield a far more robust method by which to obtain the stability signatures. The results of each simulation from which the stability signature was extracted can be found in Figures 15.

In addition to the details that are evident in Figures 4(b) and 5(a), the look-up map can also be used to assess dynamical stability in a system, with specific emphasis on the conservative stability contour map. We can interpolate between the masses we simulated to obtain the signature for a planet of any mass (within the maximum and minimum mass ratios of $\mu=1-{1024\textrm{M}_{\oplus}}/\textrm{M}_{\odot}$). We look at how to do this to achieve our two goals – assessing dynamical stability in multiple planet systems, and how to identify systems where an exo-Earth may exist in the HZ of systems – in section 5.

4.2 Effect of eccentricity

Whilst we focus primarily on co-planar, circular orbits, our method can be further developed and refined to extend to eccentric orbits. In this work, the eccentricity look-up map we generate is limited to a single mass ratio ($\mu=32\textrm{M}_{\oplus}/\textrm{M}_{\odot}$) and is not used like the ones presented in Figure 5. However, it does allow us to explore the influence of planetary eccentricity on the stability of a system. Figure 6 shows the optimistic and conservative stability signatures, and the look-up maps for a mass ratio of $\mu=32\textrm{M}_{\oplus}/\textrm{M}_{\odot}$. The $x$-axis shows the semi-major axis range simulated to obtain the stability signatures, i.e. between orbits with periods in $10:1$ and $1:10$ commensurability with the planet (which for $a_{\mathrm{pl}}=1~{}\mathrm{au}$ is given by $0.215~{}\mathrm{au}\lesssim a\lesssim 4.642~{}\mathrm{au}$), while the $y$-axis shows the eccentricity values of the stability signatures in Figure 6(a), and indicates the eccentricity of the planet in each simulation in Figures 6(b) and 6(c). We only explore between $e=0$ and $e=0.3$, given that existing literature shows that multiple planet systems (especially where one is a potentially habitable terrestrial planet) are uncommon for high eccentricity orbits ($e>0.4$). This has been explained as being most likely attributable to planet-planet scattering during planetary formation and evolution (Matsumura et al., 2016; Carrera et al., 2016; Agnew et al., 2017; Zinzi & Turrini, 2017).

The stability signatures resulting from these simulations prove more challenging to extract since the TPs become far more excited. This can be seen by plotting the excitation of the TPs in ($a$,$e$) space in Figure 15 and Figure 16 for the circular and non-circular cases respectively. While we do see the theoretical work bounding the unstable, chaotic region in Figure 6(c), there is a noticeable deviation with the optimistic signature (as evident in Figure 6(b)). This highlights two important notions: 1) as mentioned earlier, the case for investigating the effect of eccentricity means the system is no longer a circular restricted three-body problem. As such, the Jacobi integral is not conserved and so we see significantly different dynamics than in the circular case when investigating the effect of mass, and 2) the method by which we determine excitation is not well suited when investigating eccentric orbits. Since we see qualitative agreement in Figures 4(b), 5(a) and 6(c) with different analytic derivations for the onset of chaos (e.g. Mustill & Wyatt, 2012; Deck et al., 2013; Giuppone et al., 2013), a modification of the excitation criteria that utilises some of the research presented in these works should be incorporated in future work that seeks to include the eccentricity parameter space. This may yield criteria that are more suitable for the eccentric cases.

Figure 6 demonstrates the destabilising influence planetary eccentricity has on the dynamical stability of nearby objects – a result that matches what has been found in previous studies of proposed multiple planet systems (Carrera et al., 2016; Agnew et al., 2017). Figure 6 shows how rapidly the regions nearby a planet are destabilised as the planet moves from a circular to an eccentric orbit. The semi-major axes that correspond to the locations of potential stabilising mean motion resonances are also far less pronounced, as seen by the less distinct dark vertical streaks in Figure 6(b). Even more impressive is the range of destabilisation in semi-major axis in Figure 6(c), reaching out in both directions to the locations of the $3:1$ inner resonance and the $1:4$ outer resonance with the planet’s orbit in the highest eccentricity case we examined ($e_{\mathrm{pl}}=0.3$). In contrast, the unstable regions in the circular case only reached to the $5:3$ inner, and just beyond the $1:2$ outer respectively, re-enforcing the significance of even moderate orbital eccentricities on multiple system stability (Zinzi & Turrini, 2017).

4.3 Effect of inclination

Similar to our investigation into the effect of eccentricity on the stability signatures of a planet, here we explore the influence of planetary inclination on the stability of a system. Figure 7 shows the various stability signatures and the look-up maps for a mass ratio of $\mu=32\textrm{M}_{\oplus}/\textrm{M}_{\odot}$. The $x$-axis shows the semi-major axes that we simulated to obtain the stability signatures, while the $y$-axis shows the eccentricity values of the stability signatures in Figure 7(a), and indicates the inclination of the planet in each simulation in Figures 7(b) and 7(c).

What is immediately apparent from both Figures 7(a) and 7(b) is that shallow mutual inclinations (i.e. $i\leq 10\degree$) have very little influence on the stability signatures of a system. Figure 7(a) shows that the stability signatures vary very little outside of the inherent stochastic variations expected using a randomly distributed swarm of TPs. Figure 7(b) shows that the stabilising influence of mean motion resonances remains more or less consistent across the range of inclinations explored, as is evident by the near uniform vertical dark streaks. Similarly, the destabilising effect of the 2:1 and 1:2 mean motion resonances are evident across all inclinations, with the conservative signature in Figure 7(c) showing depletion in regions that align with the locations of the 2:1 and 1:2 MMRs.

While the conditions here are those for the circular restricted three-body problem, it is important to note that the dynamics will change still further when the third body is treated as a massive, rather than massless object – in other words when we move from considering the restricted three-body problem to the unrestricted three-body problem. The mutual gravitational interactions that would exist between massive bodies, rather than a massive body and massless TPs, is expected to yield more complicated dynamical behaviour. In such a scenario, inclination would have a much more significant influence on dynamical stability. However, previous studies have shown that multiple planet systems are likely to exist on shallow, near co-planar orbits (Lissauer et al., 2011a, b; Fang & Margot, 2012; Figueira et al., 2012; Fabrycky et al., 2014).

5 Applications

The first batch of TESS science observations has already resulted in the detection of two planets (Huang et al., 2018; Vanderspek et al., 2018), and this is expected to grow to several thousand throughout the lifetime of the survey (Ricker et al., 2014; Sullivan et al., 2015; Barclay et al., 2018). Here, we demonstrate how our approach can be applied for our two proposed use cases for newly discovered TESS systems.

The first use case is to conservatively assess the dynamical stability of multiple planet systems. The approach we take is to assess the stability of each planet separately to provide insight into the multiple planet system as a whole. We detail this in section 5.1. The assessments presented here are conducted using the best-fit orbital parameters inferred from observations to test the robustness of our approach. We follow this with a demonstration of how to assess stability across the region covered by the uncertainties of the orbital parameters of planets, and so how our method can provide a means to assess the stability of a system with the currently inferred planetary parameters, and if more observations are needed to better constrain the true nature of the system (e.g. Anglada-Escudé et al., 2010; Anglada-Escudé & Dawson, 2010; Wittenmyer et al., 2013; Horner et al., 2014). As our look-up map in Figures 4(b) and 5(a) are only for circular orbits, and does not take into account secular interactions which can occur over longer timescales, this is a conservative assessment that can be used as a quick, first order check of a system, enabling a rapid assessment of systems that are likely to be dynamically unstable with their current nominal best-fit orbits.

The secondary use case allows for the rapid identification of single planet (and certain multiple planet) systems where additional planets can maintain stable orbits within the HZ (with potential Earth-mass planets being of particular interest here). In this way, the expected return of several thousand newly discovered systems by TESS can be quickly assessed to identify those which could contain as yet undetected exo-Earths. We detail how we achieve this in section 5.2.

5.1 Multiple Planet Stability

To assess the stability of multiple planet systems, the stability of each planet needs to be assessed separately. Thus, for any planet, $P_{1}$, we must investigate the gravitational influence that the other planets, $P_{2}$, $P_{3}$, $P_{4}$ …, $P_{n}$ have on $P_{1}$. The same assessment must be conducted for each planet in the system, investigating the gravitational influence that the other planets $P_{1}$, $P_{3}$, $P_{4}$, …, $P_{n}$ have on $P_{2}$; the gravitational influence that the other planets $P_{1}$, $P_{2}$, $P_{4}$, …, $P_{n}$ have on $P_{3}$; and so on.

Let us consider a three planet system. Starting with $P_{1}$, we first obtain the stability signature of the other planets by combining their masses with $P_{1}$ (i.e. with respective masses $M_{2}^{\prime}=M_{1}+M_{2}$ and $M_{3}^{\prime}=M_{1}+M_{3}$), interpolating within our look-up maps, and translating the signatures to each planet’s semi-major axis $a_{2}$ and $a_{3}$ respectively. We combine these signatures to get an optimistic and a conservative signature by taking the minimum value of the $P_{2}$ and $P_{3}$ stability signatures at each semi-major axis location. We then plot these combined stability signatures, as well as where $P_{1}$ is located, in ($a$, $e$) parameter space. We can then infer if $P_{1}$ would be stable, unstable, or potentially stable based on where it lies within the combined stability signature. The inference of stability is determined as outlined in section 3: below the conservative line corresponds with the unexcited, stable regions around a planet; above the optimistic line corresponds with the excited regions around a planet; and in between corresponding with potentially stable regions. This process is then repeated for $P_{2}$ and $P_{3}$. Once all planets in the multiple system have been assessed, if any planet is located above the optimistic stability signature, the system is considered potentially unstable. Conversely, if all of them fall below the conservative stability signature, the system can be considered stable. A worked example of this process for HD 5319 can be found in Appendix A.

We test the robustness our stability signature predictions with the dynamical stability of systems found through numerical simulations. Agnew et al. (2018b) conducted such a study, simulating all multiple planet systems with a gas giant within the then known exoplanet population for $10^{8}$ years in order to assess their dynamical stability. Figure 8 shows the HD 5319 system, which Agnew et al. (2018b) found to be dynamically unstable. We can see that the red planet is located above the optimistic stability signature of the blue planet, meaning it resides in the excited, unstable regions within the system, suggesting the system overall is unstable. In Figure 9 we show a similar assessment of the 47 UMa system which was found to be potentially dynamically stable by Agnew et al. (2018b). We follow the method outlined earlier for assessing the stability of a three body system by combining stability signatures, and we can see that all three planets are located within the potentially stable region.

We can compare our stability signature predictions with all the multiple planet systems that Agnew et al. (2018b) tested numerically. As they sought to identify systems that may host hidden exo-Earths in the HZ, they had additional criteria relating to the habitability of each system that ultimately yielded a selection of 51 multiple planet systems from the then known exoplanet population. As our stability signatures only cover mass ratios $\mu\leq 1024\textrm{M}_{\oplus}/\textrm{M}_{\odot}$, and $e\leq 0.3$, this places a limitation on which systems we can test our predictions against. Of the 51 systems Agnew et al. (2018b) tested numerically, 25 systems fall within the mass ratio and eccentricity ranges we used here, and so we we can directly compare our stability signature predictions with these 25 systems. Our comparison is summarised in Table 2.

We find that our stability signatures yield agreement on the dynamical stability of a system in $32\%$ (8/25) of the systems tested numerically by Agnew et al. (2018b), no strong disagreement in any of the systems tested, and only one case of disagreement with a system found to be unstable numerically (HD 160691). This is a promising result, as our predictions are using stability signatures for circular orbits, and as highlighted in section 4.2, higher eccentricity orbits will create less stable signatures. This means that the discrepancy in the numerical and signature predictions found for the HD 160691 (the only system Agnew et al. (2018b) found to be unstable that our method did not also predict to be unstable) may be reconciled, as the signature being used to predict stability would be more stable (due to it being circular) than what it should be in reality ($e>0$). It should be noted that our stability signatures do not take into account higher order secular interactions, and so stable predictions are inherently less conclusive because of these not being incorporated. There is also agreement in the potentially stable system predictions, although for multiple planet stability it is the unstable predictions that are more useful.

By demonstrating the robustness of our predictions, we can now present the use case for assessing the stability of a planetary system. We do this by carrying out Monte Carlo (MC) simulations for a system, randomly sampling within the accepted range of values. For each simulation, we can determine a stability metric by assigning a value of [math] for an unstable system, $1$ for a stable system, and linearly interpolating the metric between [math] and $1$ between the two signatures for the potentially stable systems (taking the minimum interpolated value of all planets in a given multiple system). Taking the mean of all the MC simulations for a system yields a measure of how stable the system is with the current planetary parameters. As an example, we consider HIP 65407, for which planet HIP 65407 b has $0.172\leq a_{\mathrm{b}}\leq 0.182$ au and $0.396\leq M_{\mathrm{b}}\leq 0.46$ $\mathrm{M}_{\mathrm{J}}$ and HIP 65407 c has $0.308\leq a_{\mathrm{c}}\leq 0.324$ au and $0.73\leq M_{\mathrm{c}}\leq 0.838$ $\mathrm{M}_{\mathrm{J}}$. We run 100,000 MC simulations and yield a stability metric of $92\%$, suggesting this system should not be prioritised when determining which systems require additional observations to better constrain the planetary parameters.

It must be re-iterated that this is a conservative prediction, and the implications of some dynamical interactions may be missed. Specifically, this assessment only compares stability between pairs of planets. For systems with more than two planets, this method will not include the effect of multiple-body interactions. Such interactions could potentially destabilise a planet-pair that would, on their own, be mutually stable. As a result, it seems likely that truly stable solutions would require planets to be somewhat more widely spaced in such systems than our two-planet stability maps might otherwise suggest (e.g. Pu & Wu, 2015, 2016; Tamayo et al., 2016). Additionally, previous studies have demonstrated that the ratio of two planet’s masses has very little influence on stability. Instead, it is the cumulative mass of the planets that impacts most upon their stability (Petit et al., 2017; Hadden & Lithwick, 2018). As such, our planetary predictions between large mass planets is optimistic as the application is more suited to identifying less massive, Earth-mass companions. The most appropriate use case for our approach in assessing planetary stability is in identifying unstable systems by assessing many permutations of planetary parameters across the constrained error bars, and using that assessment to inform observers as to whether to gain more data to better constrain them.

5.2 Predicting HZ companions

Generally, detailed numerical simulations are the means by which to identify the stable and unstable regions within a planetary system (Barnes & Raymond, 2004; Raymond & Barnes, 2005; Kane, 2015). Such simulations are computationally expensive and so other methods to more rapidly predict stability within a system would be particularly useful. Here, we demonstrate how we can utilise our approach to identify where additional planets can maintain stable orbits within the HZ – with a specific focus on Earth-mass planets – in lieu of computationally expensive simulations. It is essential to test the robustness of our approach, and we do so by comparing our predictions with the standard numerical approach. Agnew et al. (2017), Agnew et al. (2018a) and Agnew et al. (2018b) have performed such simulations to assess HZ stability (with both massless TPs and $1~{}\textrm{M}_{\oplus}$ bodies) for a large sample of single Jovian planet systems. Here we draw upon the results of their simulations to compare with our predictions.

5.2.1 TP companion in HZ

We first compare our predictions with systems that were simulated numerically with massless TPs. Agnew et al. (2018a) conducted a large suite of simulations for 182 single Jovian planet systems using $5000$ TPs randomly distributed throughout the HZ, and simulated each system for $10^{7}$ yrs. Depending on the orbital parameters of the Jovian planet, these simulations took days of computational time. In contrast, our predictions can be performed in seconds. We compare the TPs that were not removed from the system by ejection or collision at the end of the simulation with the predicted unperturbed, stable regions below the stability signature using our method to assess the robustness of our approach for massless bodies.

As our stability signatures have only been determined for mass ratios $\mu\leq 1024\textrm{M}_{\oplus}/\textrm{M}_{\odot}$ and $e\leq 0.3$, this places limitations on which systems we can compare with. Agnew et al. (2018a) present 10 near-circular systems (i.e. $e_{\mathrm{pl}}<0.05$) that have mass ratios which fall within this range, and for which they have explored the stable regions in the HZ using a swarm of massless TPs. These systems are ideal candidates to demonstrate how the stability signatures can predict regions of HZ stability. To do this, we: 1) normalise the system by scaling the mass of the star and the planet, 2) interpolate between the masses on the look-up map to obtain the optimistic stability signature of each planet, 3) translate the signature of each planet to its semi-major axis, and 4) overlay the stability signature over ($a$,$e$) alongside the TPs that survived to the end of the numerical simulations ran by Agnew et al. (2018a). We can then examine if the TPs align with the area beneath the stability signature which corresponds with the unexcited, stable regions in the system. A worked example of this process for Kepler-16 can be found in Appendix A. Figure 10 shows the stability signatures and surviving TPs as found by Agnew et al. (2018a) for Kepler-16 and kappa Cr B, two systems where the known planet interacts significantly with the HZ. All 10 of the near-circular systems can be seen in Appendix B.

In these figures, the coloured dots represent the final position and the relative change in semi-major axis of the TPs that remained at the end of the simulation, the green curve is the stability signature, the shaded green area is the HZ that exists beneath the signature, and the blue point is the massive planet. It can be seen that there is strong agreement between the stability signatures and the surviving TPs for each system as was determined numerically by Agnew et al. (2018a). Especially so in the case of Kepler-16 as it is the most circular of the two systems. kappa Cr B also shows strong agreement in that a large majority of the TPs fall below the stability signature, but as the orbit of the planet is slightly less circular (0.0069 in Kepler-16 vs 0.044 in kappa Cr B) we see some excitation of the TPs to eccentric orbits greater than $e>0.3$. Regardless, for near-circular systems we see that the stability signature is a strong predictor of stable regions, and that further development to incorporate eccentric orbits has the potential to yield even stronger predictions for less circular orbits.

5.2.2 Earth-mass companion in HZ

As one of the goals of our work is to identify where Earth-mass planets could maintain stable orbits, we need to demonstrate the stability predictions are not just valid for massless TPs, but also for massive bodies. We have the same mass ratio and eccentricity constraint as we did for the TP predictions, and so we similarly focus on systems with mass ratios $\mu\leq 1024\textrm{M}_{\oplus}/\textrm{M}_{\odot}$ where the planet moves on near-circular orbits. Agnew et al. (2017) conducted a suite of simulations to explore the stable regions in the HZ of 15 single Jovian planet systems by sweeping a $1~{}\textrm{M}_{\oplus}$ body over the ($a$,$e$) parameter space. For each system, 20,400 individual simulations were run, where the orbital parameters of the $1~{}\textrm{M}_{\oplus}$ body were gradually varied to cover the desired parameter space. Such thorough numerical simulations take days to weeks to complete, whereas our predictions can be performed in seconds. We compare the $1~{}\textrm{M}_{\oplus}$ bodies that were not removed from the system by ejection or collision at the end of the simulations with the predicted unperturbed, stable regions below the stability signature using our method to assess the robustness of our approach for Earth-mass bodies.

None of those systems studied by Agnew et al. (2017) would be considered near-circular, so we just look at those that have mass ratios that fall within the $\mu$ range and with an $e_{\mathrm{pl}}\leq 0.3$. Three systems fulfill these requirements, HD 132563 B, HD 147513 and HD 34445, and hence are candidates to demonstrate how the stability signatures can predict HZ stability. While these systems have relatively low eccentricities ($0.22$, $0.26$ and $0.27$ respectively), they are not what one would consider near-circular. As such we also include a multiple planet system where the planet in closest proximity to the HZ is much closer to circular. That system is 47 UMa, for which Agnew et al. (2018b) carried out a similar massive body parameter sweep. In this case the separations between planets in $47$ UMa are such that the predictions are still useful. To perform these HZ predictions we perform the same steps as outlined in section 5.2.1, but instead of plotting the TPs that survived to the end of the numerical simulations, we plot the $1~{}\textrm{M}_{\oplus}$ bodies that survived to the end of the numerical simulations performed by Agnew et al. (2017) and Agnew et al. (2018b). Figure 11 shows the stability signatures and surviving $1~{}\textrm{M}_{\oplus}$ bodies found by Agnew et al. (2017) and Agnew et al. (2018b) for HD 132563 B, HD 147513, HD 34445 and 47 UMa.

An additional step required for 47 UMa is that the stability signatures of all planets must be combined. At any semi-major axis location, the stability signature of each planet will have a corresponding eccentricity representing the maximum eccentricity of the unexcited region. As we assess whether one planet would excite any other, the smallest of these eccentricity values will ultimately bound the unexcited region at that semi-major axis location. As such, we can take the lowest value of all the stability signatures at each semi-major axis location to generate a combined stability signature. The combined signature for 47 UMa is that shown in the bottom panel of Figure 11.

In these figures, the coloured dots represent the final position and the relative change in semi-major axis of the $1~{}\textrm{M}_{\oplus}$ bodies that remained at the end of the simulation, the green curve is the stability signature, the shaded green area is the HZ that exists beneath the signature, and the blue, red and yellow points are the massive planets.

As the known massive planet in HD 132563 B, HD 147513 and HD 34445 is eccentric, the stable bodies do not fit the predictions as accurately as the near-circular cases. We can see that while the massive bodies tend to align along the semi-major axis, they are excited to higher eccentricity orbits. This matches what we see in the eccentric simulations as described in section 4.2 and with what can be seen in Appendix C. However, it does demonstrate that the predictive ability of the stability signatures is still valid for a $1~{}\textrm{M}_{\oplus}$ body. In the case of 47 UMa, as the planets are adequately separated the stability signature of the planet nearest to the HZ can be used to predict stability within it. That planet moves on a much more circular orbit, and so we can see that the stable $1~{}\textrm{M}_{\oplus}$ bodies align very closely with the area below the stability signature. This agreement of the stability signature with $1~{}\textrm{M}_{\oplus}$ bodies shows that our method can be used effectively to identify where stable, HZ Earth-mass planets could maintain stable orbits in known exoplanet systems without the need to run exhaustive numerical simulations.

It is important to note that Figures 10 and 11 demonstrate that several bodies do remain stable at the end of the numerical simulations carried out by Agnew et al. (2017), Agnew et al. (2018a) and Agnew et al. (2018b) that have eccentricities greater than $0.3$. The limit of $e\leq 0.3$ that we enforce here to find the stability signatures and develop our look-up map is the result of our focus on habitability (Williams & Pollard, 2002; Jones et al., 2005). As such, it should be kept in mind that our method only predicts the unperturbed regions in ($a$,$e$) space with $e\leq 0.3$, even though it is possible for orbits with higher eccentricities to be stable.

6 Conclusions

Numerical simulations are integral to assessing the detailed dynamics of planetary systems. However, alternative methods to classify system stability are beneficial in ensuring computational resources are efficiently allocated to assess only the most complicated systems. In this work, we have presented an alternative approach in assessing the stability of newly discovered exoplanet systems, such as those that will be found in coming years by TESS, and its associated follow-up programs. This includes the dynamical stability of a multiple planet system with the best-fit orbital parameters, the overall stability of the HZ of a system, and whether an Earth-size planet could maintain a stable orbit within the HZ. The key findings of our work are as follows:

• •

Mass ratio ($\mu=M_{\mathrm{pl}}/M_{\mathrm{\star}}$) and orbital eccentricity are very influential in determining a system’s overall dynamical stability. In particular, we find that even moderate orbital eccentricities can prove to be destabilising, a finding that re-enforces the results of several studies that highlight the inverse relationship between multiplicity and eccentricity (Carrera et al., 2016; Agnew et al., 2017; Zinzi & Turrini, 2017).

• •

“Stability signatures” can be obtained using our method for each planet in a multiple planet system, and these signatures can be used to assess the stability of each planet, and to determine the overall dynamical stability of a particular set of planetary parameters for a multiple planet system. Comparing our predictions with previously run numerical simulations using best-fit orbital parameters (Agnew et al., 2018b), we found our approach yields no strong disagreement in any of the systems assessed, and agreement in $32\%$ of cases.

• •

The stability of a planetary system can then be investigated by carrying out the stability signature assessment for many different permutations of planetary and orbital parameters. By considering the parameters over their respective error ranges, one can suggest whether more observations should be taken in order to better constrain the system to more stable orbital parameters.

• •

The stability signature interpolated from our results also proves effective at predicting the stability of massless TPs in near circular systems. Comparing the signature of several single planet systems with those TPs that were found to be stable with numerical simulations by Agnew et al. (2018a) shows strong agreement.

• •

The stability signature is also good at predicting the stability of a $1~{}\textrm{M}_{\oplus}$ planet in the HZ of low eccentricity systems. Comparing the stability signature of several single planet systems with stable $1~{}\textrm{M}_{\oplus}$ bodies found numerically by Agnew et al. (2017), we can see good agreement, with discrepancies being due to the stable $1~{}\textrm{M}_{\oplus}$ bodies being more excited to higher eccentricities as a result of the existing single planet not being on a near-circular orbit (an assumption of our model).

• •

The stability signature is also very good at predicting the stability of a $1~{}\textrm{M}_{\oplus}$ planet in multiple planet systems where the separation between known planets is such that the HZ only interacts with the nearest planet. $47$ UMa provides such a case where the planet nearest to the HZ is near circular, and here we again see strong agreement between the stability signature and the numerical simulations performed by Agnew et al. (2018b).

Our work has focused on the simplest, near circular, co-planar case, and so there is room to refine our approach, which we intend to do in the future. We have demonstrated that our method shows a high degree of success for low eccentricity systems, and for multiple planet systems with large orbital separations.

This approach will be particularly useful for systems discovered with TESS, allowing the system stability to be assessed for the best-fit orbital parameters and informing whether additional observations should be made to further constrain the orbits. It can also be used to rapidly predict which newly discovered systems may have dynamically stable Earth-size planets orbiting in their HZ. As our predictive tool does not require further numerical simulations, they can be incorporated into the Exoplanet Follow-up Observing Program for TESS (ExoFOP-TESS) to provide these insights as planets are discovered.

Acknowledgements

We wish to thank the anonymous referee for their thoughtful report, and helpful comments and recommendations that have improved the paper. We wish to thank Jessie Christiansen for helpful discussions regarding TESS and ExoFOP-TESS. MTA was supported by an Australian Postgraduate Award (APA). This work was performed on the gSTAR national facility at Swinburne University of Technology. gSTAR is funded by Swinburne and the Australian Government’s Education Investment Fund. This research has made use of the Exoplanet Orbit Database, the Exoplanet Data Explorer at exoplanets.org and the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program.

Appendix A Worked Examples

A.1 Multiple Planet Stability

Here we show a worked example for assessing the dynamical stability of HD 5319. The current known arrangement of HD 5319 is summarised in Table 3. We follow the steps as they are outlined in section 5.1.

We first normalise the system. For a star with mass $M_{\star}=n\mathrm{M}_{\odot}$, we scale the masses of both the star, $M_{\star}$, and the combined mass of the planet pair, $M_{\mathrm{pl,1}}^{\prime}=_{\mathrm{pl,2}}^{\prime}=M_{\mathrm{pl,1}}+M_{\mathrm{pl,2}}$, by $1/n$. This yields a normalised planet mass of

[TABLE]

With the normalised mass, we are now able to obtain the stability signature of the each planet (which is the same in the two planet case) from our look-up maps by interpolating between the masses we have simulated. Figure 12(a) shows where the combined normalised mass of the planets “slices” across our two maps. The shading of the look-up map corresponds with the maximum unexcited eccentricity for the optimistic stability signature, and the minimum excited eccentricity for the conservative stability signature, each slice of the contour-map will have a corresponding curve. These are shown in Figure 12(b).

Having obtained the stability signatures for the planets, we next need to translate the signatures to the semi-major axes of each planet. The domain over which our simulations were run is

[TABLE]

As per equation 3.3, the domain of our translated signature will be

[TABLE]

This yields a domain for planet b of

[TABLE]

[TABLE]

and a domain for planet c of

[TABLE]

[TABLE]

The translated signatures for each planet can be seen in Figure 12(c).

To assess stability, we identify where each planet is located relative to the stability signatures of the other. Figure 12(c) demonstrates that in the case of HD 5319, planet c falls outside the optimistic stability region of planet b, and so is excited and so it may be unstable.

A.2 Predicting HZ companions

Here we show a worked example of predicting HZ companions in Kepler-16. The current known arrangement of Kepler-16 is summarised in Table 3. We follow the steps as they are outlined in section 5.1.

We first normalise the system. For a star with mass $M_{\star}=n\mathrm{M}_{\odot}$, we scale the masses of both the star, $M_{\star}$, and planet, $M_{\mathrm{pl}}$, by $1/n$. This yields normalised mass for planet b of

[TABLE]

With the normalised mass, we are now able to obtain the stability signature for planet b from our look-up map by interpolating between the masses we have simulated. Figure 13(a) shows where the normalised mass of the planet “slices” across our map. For predicting HZ companions, we use the optimistic stability signature. The shading of the look-up map corresponds with the maximum unexcited eccentricity for the optimistic stability signature. This is shown in Figure 13(b).

Having obtained the stability signature of planet b, we next need to translate the signature to its semi-major axis. The domain over which our simulations were run is

[TABLE]

As per equation 3.3, the domain of our translated signature will be

[TABLE]

This yields a domain for planet b of

[TABLE]

[TABLE]

We then compute the HZ boundaries for the system. We use the conservative HZ definition for a $1~{}\mathrm{M}_{\oplus}$ as given in Kopparapu et al. (2014) for our HZ boundaries. Figure 13(c) shows the translated signature and the HZ that falls below it for Kepler-16.

Finally, using the results of the numerical simulations run by Agnew et al. (2018a), we can plot the final position of the TPs remaining at the end of their simulation in ($a$,$e$) space. Figure 13(d) demonstrates there is tight agreement between the stability signatures and the numerical simulations.

Appendix B HZ companion Figures

Appendix C Test particle Figures