1 Introduction↩︎

Despite the rapid progress in discovering extrasolar planets, the exoplanet population of the solar neighborhood remains mostly unexplored: for example, around the approximately 1,500 stars known within 15 pc [1], which may host several thousand planets, only 159 are known. Given that planets are intrinsically faint, those that are the closest to us (i.e., orbit the closest stars) are set to play key roles in the future of exoplanet studies and, especially, in searches for atmospheric biosignatures [2], [3]. Therefore, it is strategically important to explore nearby planetary systems and to identify those that are the most likely to harbor potentially habitable worlds. While the information on specific planetary systems almost always remains incomplete and often uncertain, the emerging body of information on exoplanet demographics and on planet formation and evolution can be combined in the future to improve the understanding of individual planetary systems [4].

The goal of this paper is to study the planetary system of (also known as 82 Eridani, , and ; not to be confused with \(\epsilon\) Eridani), the fifth-closest Sun-like star, via a novel integrative analysis. Through this analysis, we predict the orbits and properties of yet-undetected planets. Such predictions can help the interpretation of planet candidates identified at low significance levels and can also guide follow-up observations of the planetary system.

At a distance of \(6.04\) pc, is one of the nearest and brightest stars in the sky that is known to host super-Earth planets [5]. The first evidence for planets in the system was found in radial velocity data by [6], who confirmed two planets and suggested the presence of another. Follow-up observations by [7] provided additional evidence for several other planets, with up to six potential planets and candidates known in the system.

The planetary system around not only has an astronomical significance, but it also has a special place in science fiction and visions of humanity’s future. It is featured in various novels and short stories as a system with habitable worlds which foster intelligent life (c.f., Anderson’s Orbit Unlimited and Pournelle’s short story He Fell into a Dark Hole, or Kube-McDowell’s story that envisioned it as a home to a primitive human colony). These stories predated the discovery of planets in the system, but modern methods finally enable exploring the planetary architecture of the system and the properties of its worlds.

In this paper, we integrate specific information on the planetary system (constraints on planets and their orbits) with population-level information on exoplanets (exoplanet demographics) to provide a robust assessment of the planetary system. Specifically, we use the DYNAMITE¹ [8] integrated analysis software package, which was recently successfully utilized in multiple studies (see below). DYNAMITE asks the question: If an additional planet exists in a planetary system, what are the probability distributions of its orbital period, eccentricity, inclination, and planet size? To answer this question, DYNAMITE uses robust trends identified in the Kepler exoplanet demographics data (orbital period distribution, planet-size distribution, etc., based on the \(\sim\)​2,400 exoplanets that form the Kepler population) with specific data for a given single exoplanet system (detected planets and constraints on their orbits and sizes). Based on this information, DYNAMITE uses a Monte Carlo approach to map the likelihood of different planetary architectures, also considering orbital dynamical stability and allowing for freedom of statistical model choice.

By now, DYNAMITE has been used successfully in a multitude of studies: [8] first demonstrated the capabilities of DYNAMITE on 45 TESS-discovered transiting planetary systems to predict yet undiscovered planets. In that study, when a known transiting planet – Kepler-154 f – was deliberately hidden from the algorithm, DYNAMITE not only retrieved the planet but also predicted its size and inclination accurately. When the known non-transiting planet Kepler-20 g was removed from the Kepler-20 system, DYNAMITE again recovered the planet. When removing multiple planets from the TOI-174 system, DYNAMITE again indicated the likelihood of an additional planet at the location of the purposely removed planets. Even when the periods of other planets are altered within \(3\sigma\), the predictions are mostly unaffected by this, and only slightly shifted if all the planetary periods were changed in the same direction. Lastly, when analyzing the HD 219134 system, another planet with a period \(\sim\) 12 days is predicted between current planet c and candidate f, along with an additional planet external to planet candidate g (see also Dietrich et al. 2021, submitted).

When DYNAMITE was applied to the \(\tau\) Ceti system [9], it provided contextual and statistical evidence for unconfirmed planets b, c, and d, which were found in one earlier study but have marginal evidence over a decade of observations [10], [11]. In addition, DYNAMITE also predicted another yet-unobserved planet in the habitable zone. [12] used DYNAMITE on the K2-138 planetary system and accurately predicted planets c and e while analyzing different scenarios in their study. In tests focusing on the inner Solar System, DYNAMITE also successfully predicted the positions and sizes of Venus and Earth when one of them was removed from the inner Solar System [13]. While DYNAMITE provides a robust way of predicting the likely positions where additional planets could exist in a system, we emphasize that the analysis is based on the available exoplanet demographics data. Therefore, the validity of extrapolating from these results to planet mass \(< 1 M_\oplus\) or \(> 4 M_\oplus\) and to non-planetary objects (like planetesimal belts) could not be ascertained. Consequently, our analysis and predictions from DYNAMITE are likely to be correct for typical planetary systems but not expected to hold for rare planets or outliers among planetary architectures.

In addition to applying DYNAMITE to , in this study we also present an additional module that statistically determines the possible equilibrium and surface temperatures of the planets studied. We also compare the temperature results for the system and the inner Solar System. Though the model we use for this specific task is simple, it still provides some useful insights into the possible temperatures of the planets in the system.

This paper is organized as follows. Section 2 delineates the current knowledge of the system. We describe the DYNAMITE package and our implementation of it, the predictive temperature analysis, and the eccentricity analysis in Section 3. Section 4 summarizes the results of our analysis of the planetary architecture, the equilibrium and surface temperatures of the planets, and the results of our eccentricity analysis. We discuss the results of the architectural analysis (including the likely number of planets and their eccentricities) in Section 5. In Section 6, we motivate our prior on the planetary Bond albedo distribution and then discuss the results of the predicted equilibrium and surface temperature likelihood distributions, as well as the potential for habitable worlds in the system. Finally, in Section 7, we discuss the potential for habitable worlds in the system.

2 The e Eridani planetary System↩︎

2.1 Stellar Parameters↩︎

is a G8V, broadly Sun-like star located at a distance of \(6.0414 \pm 0.0029\) pc [14], [15] from Earth, with a luminosity of \(0.656 \pm 0.003 L_\odot\) [16]. Its key stellar parameters are summarized in Table 1.

Notes: (a) [6], (b) [17], (c) [18], (d) [16], (e) [15], (f) [19]

[17] calculated the mass of to be \(0.813_{-0.012}^{+0.018} M_\odot\) based on stellar evolution models. The stellar radius of \(0.92 \pm 0.02 R_\odot\) was calculated by [18] through a dense grid of evolutionary tracks. The calculated stellar radius value lies within the range of stellar radii calculated by [20] for this star based on Gaia parallaxes and the Stefan–Boltzmann Law. We adopted the derived luminosity, metallicity, and effective temperature values from [16], who carried out a detailed spectroscopic analysis of the HARPS Guaranteed Time Observations (GTO) “high precision” data. We note that a \(4.68\%\) higher luminosity than we adopted from [16] was calculated by [15], but this slightly higher value would not appreciably impact the boundaries of the habitable zone or any of our results. Though ’s stellar spectral classification is not yet fully settled, we note that – based on the values of mass, radius, luminosity, and temperature – is more likely to be a late G-type star than an early K-type star. However, our analysis relies on the derived luminosity and effective temperature values and not on the spectral type. The age of is likely old, but its exact value remains uncertain: \(5.76\) Gyr is reported by , an age \(> 12.08\) Gyr is found by [18], and \(13.5\) Gyr is derived by [19]. The latter two ages have been determined by isochrone fitting, while used the activity-rotation-age calibration method to calculate the age, following the method given by [21]. While such a large disparity in the star’s age could affect the potential habitability of the system, the age itself is not utilized in our statistical analysis.

2.2 Planetary system↩︎

In 2011, announced the discovery of two planets orbiting with orbital periods of 18 days (planet b) and 90 days (planet d). These planets were found using radial velocity (RV) measurements and were identified as likely super-Earths, given their measured \(m \sin i\) values. The analysis also found weak evidence for an additional planet, planet candidate c, at a 40-day period. In a follow-up study by , six different signals were found at periods of 11.9 d (planet g), 18.33 d (confirming planet b), 43.17 d (planet c), 88.9 d (confirming planet d), 147.02 d (planet e), and 331.41 d (planet f). The three additional signals found in the [7] study resulted from an extended set of observations, more sophisticated data analysis, and superior noise model. Of the six radial velocity peaks reported, three (belonging to planets b, d, and e) were strong while the other three needed further analysis to confirm them. This newer analysis found marginal evidence for the original planet candidate c from at a slightly longer 43-day period.

The planetary signal at an orbital period of 12 days corresponding to planet g is relatively weak, but it consistently appears in the majority of the datasets. Lastly, the planetary signal at 331.41-day period, planet f, is estimated to lie in the habitable zone of the system, as calculated by [22], [23] for a 5 \(M_ \oplus\) planet and 1 \(M_ \oplus\) planet around . state that the signal for planet f needs more observation to be confirmed, as the signal is close to the annual sampling period from Earth. The proximity of this orbit to a 2:1 resonance with planets b and d, while not actually residing in the resonant orbit, provides some support for its planetary nature, as [24] found that the most common period ratio between neighboring pairs of planets in the Kepler population was \(1.9\). Similarly, [25] found that 2:1 resonances between exoplanets are not intrinsically rare.

The eccentricities of the planets in this system remain poorly constrained. found that the uncertainties for the eccentricities were of similar order to the values and therefore fixed all eccentricities to 0 (circular orbits). allowed for both circular and Keplerian (non-zero eccentricity) orbits and found that moderate eccentricities (\(\sim\)​0.1-0.3) fit their data best. However, they also state this is likely due to instrumental noise and fitting bias, and so the actual eccentricities are likely lower: Indeed, the orbits of planets d and e would cross and very likely be unstable given the currently reported values. Later in this paper we explore the possible range of eccentricities that lead to stable systems.

We summarize the relevant parameters of the six planets and planet candidates in Table 2.

Note: These planetary parameters have been collected from [7]

.

Although not used in our study, for completeness we note that the system was also found to host a debris disk, revealed by infrared excess detected with the Spitzer Space Telescope [26]. Follow-up observations with the Herschel Space Observatory marginally resolved the disk emission and the disk’s inner and outer radii were estimated to be around 20 and 30 au, with the disk relatively depleted beyond 30 au [27]. Modeling of the disk spectral energy distribution argued for a cold disk (T\(_{eq}=80^{+70}_{-30}~K\)) with blackbody disk radius about 10 au, based on disk temperature [27]. The debris disk provides evidence for a planetesimal belt exterior to the known planets, similarly to the Solar System’s Kuiper Belt. In the future, studies of this planetesimal belt may provide constraints on the outer planetary system in the future.

3 Analysis↩︎

3.1 The DYNAMITE Package↩︎

Here we briefly explain the functionality of the DYNAMITE software package and the mode that we use in our study. As stated earlier, DYNAMITE amalgamates various statistical distributions identified in the Kepler exoplanet population data and combines that general information with the specific data of a given planetary system. It then integrates this with a dynamical stability criterion and tests the probability density functions (via Monte Carlo injections) for orbital period, planetary radius, and planetary inclination. Finally, DYNAMITE also considers the orbital stability of the different planetary architectures. This combined information allows the algorithm to estimate the likelihoods for the parameters of one additional unknown planet in the system. The basic assumptions we used to constrain the probability distribution in DYNAMITE are:

1. The period, planet radius, and orbital inclination are independent of each other.

2. The orbital periods considered range from 0.5 days to 2 years, corresponding to the approximate range for which strong constraints exists on planetary architectures and planet properties from the Kepler mission-discovered exoplanet population.

3. The planet radius range considered extends from 0.5 to 5 \(R_\oplus\). Here the lower limit comes from the Kepler population and the upper limit empirically contains a large majority of the TESS multi-planet systems. [8]. Planets larger than 5 \(R_\oplus\) tend to be rare and may have distinct orbital period, size, and mass distributions [28].

4. The range of allowed planetary inclinations is 0\(^{\circ}\)– 180\(^{\circ}\).

5. The dynamical stability is determined either via N-body simulations and the spectral fraction analysis from [29] or via the mutual Hill parameter for neighboring pairs of planets \(\Delta_c \geq 8\) [30], depending on user input.

For a detailed discussion of these assumptions, please see [8]. In our analysis, we used the period ratio model from [24], which assumes a broken power law for location of the first planet in period space (with the break at 12 days) and equal period ratios between pairs of planets. Therefore, this model inserts planets roughly symmetrically in gaps as a consequence of this chosen distribution, with the highest probability occurring when the period ratios are closest to the Kepler population period ratio mean of 1.9.

The inclination distribution, in the case of no known planetary inclinations, is isotropic but limited to non-transiting values for each planet (as the planets are known not to transit). The planetary radius distribution is a Lognormal distribution around a central cluster, which is fit to the radii of the known planets, or if the planet radii are unknown, from the assumed radii derived from the known mass or m sin i values and the [31] M-R relationship). After injecting a planet, we calculate the dynamical stability parameter and compare it to the threshold value to see if the system would stably exist with that planet addition.

Through 10\(^6\) Monte Carlo iterations, we build up the likelihoods for the additional planet parameters. While the results are reliant on the specific data for the given system being accurate, they are not sensitive to the exact values. We utilize the default parameter set for each planet from the NASA Exoplanet Archive², and the uncertainties on most of these parameters are low enough to enable a robust analysis. We explore the impact of uncertainties on the DYNAMITE output predictions in Section 5.1 and Appendix 10.

3.2 Planetary Configuration Hypotheses↩︎

Given the uncertain nature of three planet candidates in the system, we proceeded in our analysis by constructing hypotheses that explore all possible combinations of the planet candidates being genuine or not. For each hypothesis, we provide DYNAMITE predictions. As explained later, we found that the predicted planetary orbits that results from these different hypotheses are very similar. In the main body of the current paper we present four hypotheses, while in the Appendix 10 we show that results of the additional hypotheses. The hypotheses are constructed by starting from the three confirmed planets, and by adding planet candidates one by one, in the order of decreasing evidence for their existence. Thus, our first hypothesis (H[a]) includes only the three confirmed planets (b, d, e), and our second hypothesis (H[b]) adds to this set the planet candidate for which evidence is the most convincing (planet c, detected in two papers). Then, in our third hypothesis (H[c]), we included planet candidate f as a genuine planet (stronger candidate signal detected by ). Finally, we assume that all planet candidates are genuine planets (Hypothesis H[d]).

We note that, in the above sequence of hypotheses, predictions from each of the hypotheses motivates the next one, as described (i.e., H[a] predicts a planetary architecture consistent with H[b], etc.). As discussed in §5, our predictions for the planetary architecture are robust and do not display strong dependence on the presence or absence of any single planet. Our analysis of the planetary architecture assesses different scenarios to reflect the yet unsettled nature of the candidate planets. We organize these possibilities into four hypotheses, corresponding to four different potential states of the planetary system given current constraints:

• Hypothesis H[a]: Only the three planets b, d, and e are taken as known planets. None of the three planet candidates (g, c, f) are genuine detections.

• Hypothesis H[b]: In addition to the three planets (b, d, and e), planet candidate c is assumed to be a genuine planet. Planet candidates f and g are not considered to be real detections.

• Hypothesis H[c]: Planets b, c, d, e, and f are all genuine detections, but the peak corresponding to planet g is not a detection.

• Hypothesis H[d]: All six planet detections (g, b, c, d, e and f) proposed by are real planets.

We do not fix the inclinations of the planets and planet candidates, as those values are unknown due to the inability to measure them with RV observations alone, although we can constrain the inclination upper bound via the determination that these planets do not transit. We give the entire range of inclination possibilities to DYNAMITE, thus lowering the chances of biasing the model towards any specific planetary architecture type. The orbital eccentricity of the planets was also input based on multiplicity-dependent statistical distributions from the Kepler population [32]. Non-zero eccentricities and the statistical distribution of orbital eccentricities are explored further in [13].

We also predict the planet radius and mass of additional planets for each hypothesis. We use a clustered planet radius model from [28] that fits the currently known planet radii to a Lognormal distribution, from which the values are drawn for the Monte Carlo iterations of the additional planets. These planet radii are also converted into mass values via the mass-radius relationship from [31], specifically using the “volatile-rich" power-law for planets in the mass range \(5-25 M_\oplus\), with planets below \(5 M_\oplus\) assumed to be rocky and planets above \(25 M_\oplus\) assumed to contain a large gaseous envelope. We expect the planets in this system of sufficient mass orbit far enough away from that they would accrete and retain a significant gaseous envelope, instead of it being stripped away via photoevaporation [33].

3.3 Surface Temperature Assessment↩︎

With 3-6 planets possible in the planetary system of , the potential equilibrium and surface temperatures of these worlds are of great interest. As available data are scarce, we will follow a statistical approach to explore the potential range of temperatures on these worlds. We do so by expanding our DYNAMITE modeling framework with statistical predictions of the equilibrium and approximate surface temperatures. Our analysis should be taken only as an initial exploratory assessment, until more constraining data on the planets are available; nevertheless, as shown later in the manuscript, our analysis provides interesting insights into the planets’ possible natures. In the following text we will describe the methodology used to provide the initial, exploratory assessment of the potential range of planetary temperatures for the predicted planets.

The calculation of \(T_{eq}\) is based on two assumptions: first, the planet (with its atmosphere) is a blackbody, and second, the complete system is in radiative equilibrium with its surroundings. This is a very good approximation of the average temperature (i.e., disregarding day/night side variations) for an airless body (e.g., the Moon, Mars), and is accurate within 10% for planets without an atmospheric greenhouse effect. For , we adopt the \({a_{planet}}\) (the semi-major axis) and \(P_{planet}\) (orbital period) values for the planets and planet candidates from the recent RV study by .

For every orbital period from the Monte Carlo iterations from DYNAMITE, we draw an albedo from a normal distribution with mean \(\mu = 0.3\) and standard deviation \(\sigma = 0.1\) and calculate the radiative equilibrium temperature of the planets via Eq. 1 (\(A_{B}\) is the Bond albedo, \(T_{star}\) and \(R_{Star}\) are the surface temperature and radius of the star respectively, and \(a_{planet}\) is planet’s semi-major axis). \[\label{eq:t} T_{eq} = T_{star}(1-A_{B})^\frac{1}{4}\sqrt{\frac{R_{star}}{2a_{planet}}},\tag{1}\] A detailed discussion of our choice of a Gaussian distribution for the Bond albedo with the above mentioned parameters is given in Appendix 9.

We note here that the available data on exoplanet albedos is very sparse, and there remain uncertainties on the mass-radius conversion. These may lead to systematic offsets in our albedo model. Nevertheless, our results are not very sensitive to the exact albedo values, in part because our assumed albedo distribution is relatively wide. For example, the assumption that the albedo of planet b in the system is 0.6, instead of 0.3, results in a temperature decrease of \(\sim13\%\) from the equilibrium temperature calculated with albedo of 0.3. Moreover, as the albedo distribution we adopted is relatively broad and is derived using the bodies in the Solar system, few objects are expected to lie outside of this range.

Perhaps the most important limitation of our simple model is the lack of an atmosphere and the associated greenhouse effect. To address this, we calculated the surface equilibrium temperatures using a leaky greenhouse atmosphere model [34], wherein we assume a single-layered atmosphere (which approximates to a blackbody) that absorbs a fraction \(\alpha\) of the radiation from the planet. Our single-layer leaky greenhouse model is a reasonable first approximation for planets without heavy atmospheres and very strong greenhouse warming. Additionally, more comprehensive models could provide more realistic predictions of the surface temperatures, but would – in turn – also require assumptions on the atmospheric compositions and cloud properties which are unconstrained for most exoplanets. Therefore, a simple model with just one free parameter is a better choice. We use a normal distribution for \(\alpha\) (\(\mu = 0.5\) and \(\sigma = 0.1\)) so that nearly \(100\%\) of the values lie between 0 and 1. Due to the isotropic nature of the thermal emission, the atmosphere then radiates equal amounts of energy towards the planet as well as away from it. The energy balance of the stellar irradiation, surface emission, and single-layer greenhouse absorber can then be described as follows:

\[\label{eq:322} T_{surface} = \left(\frac{2}{2 - \alpha}\right)^\frac{1}{4}T_{eq},\tag{2}\]

where \(T_{\text{surface}}\) is the surface temperature.

To report the most likely equilibrium and surface temperatures for planet candidates, we select the median along with the 16th and 84th percentiles of orbital period distribution from DYNAMITE and convert those into the mean temperatures and their standard deviations. Thus, our temperature measurements mostly depend on the uncertainties in the predictions of the location of planet candidates by DYNAMITE; the uncertainties of the orbital period are (relatively) much smaller than the uncertainties typical to the other orbital parameters or planetary masses. The results are plotted as a 2-D histogram in Figure 4 and Figure 5.

For statistically predicting the equilibrium and possible surface temperatures of planets (b, d, and e) in the system, we plot another histogram using a linear orbital period distribution centered at the planet’s proposed orbital period and extends as per the uncertainties proposed by . With the orbital period prior, we use the same approach of calculating relative likelihood of surface temperatures with normally drawn albedos between 0 and 1 with \(\mu = 0.3\) and \(\sigma = 0.1\).

3.4 The Eccentricity Assessment↩︎

With multiple planets on close orbits, the analysis of the long term stability of the system can provide important insights. acknowledged that the planetary eccentricities in their best-fit multi-planet radial velocity models were, in fact, moderate-to-high, which may cause instability in the system. Our goal here is to evaluate the stability of the system with the eccentricities provided by and provide updated bounds on the eccentricities of the planets in the system. In order to analyse the stability of the system, we use REBOUND [35], [36] software package to run N-body simulations with two different eccentricity distributions: (1) a Normal distribution for each planet’s eccentricity with parameters provided by and (2) a Lognormal distribution for each planet’s eccentricity with parameters provided by [32].

The dynamical stability criterion that we use for each randomly drawn eccentricity checks the system for four different conditions:

1. The orbits of neighboring planets crossing each other before the N-body integration begins.

2. Close encounters and ejections during N-body integration.

3. The orbits of neighboring planets crossing each other during N-body integration.

4. Spectral Fraction [29], as explained below.

The spectral fraction analysis executes short N-body integrations (\(\sim 5\times 10^{6}\) orbits of the innermost planet) to predict the long-term (\(\sim 5\times 10^{9}\) orbits of the innermost planet) stability of the system. This is done by determining the “spectral fraction”, or the fraction of frequencies having power \(\geq5\%\) of the peak frequency in the power spectrum of the angular momentum deficit (AMD, the difference between the total angular momentum and planar circular angular momentum). If the spectral fraction is \(\leq5\%\), then the system is assumed to be likely stable long-term. Else, the system is assumed to be unstable long-term.

After assessing the likely planetary architectures for using DYNAMITE, we identify two limiting architectures for the system – Hypotheses H[a] and H[d] – which drive the loosest and the tightest constraints on the planetary eccentricities. Therefore, we explore these two specific architectures – the 3-planet system and the 6-planet system – to help us analyze the potential planetary eccentricities for this system. For each simulation, we iterate over 1000 randomly drawn eccentricities from the chosen distributions mentioned above and report the number of stable and unstable iterations. For every stable and unstable iteration, we calculate the combined average eccentricity of all the planets involved in the run. We then report the combined mean eccentricity for the planets in the system and associated standard deviation of all the cases that resulted in a stable and an unstable configuration in the simulation.

4 Results↩︎

4.1 System Architecture Analysis↩︎

4.1.0.1 Hypothesis a: planets b, d, and e

In our first hypothesis, Hypothesis H[a], we only placed the planets b, d, and e in the system. We found two roughly equally significant peaks in the probability distribution, H[a]–PxP–1 (i.e., Predicted eXoPlanet 1 in Hypothesis H[a]) at 10.9 [7.7, 12] (\(16^{th}\) and \(84^{th}\) percentiles) days (interior to planet b) and H[a]–PxP–2 at 40.7 [34.2, 52.4] days between planets b and d. A third, slightly larger probability distribution peak is identified as H[a]–PxP–3 at 271 [269, 520] days, exterior to planet e.

4.1.0.2 Hypothesis b: Accepting planet candidate c

We considered planet candidate c as a real planet in the system. In our analysis of this hypothesis, we recovered peaks H[a]–PxP–1 and H[a]–PxP–3 as H[b]–PxP–1 and H[b]–PxP–2, respectively. In addition to these two peaks, we also found two small regions of non-zero probability in the period space, between planets b and c, and planets c and d, but the injection probability of any additional planet in those gaps is relatively low. Therefore, this hypothesis predicts two planets, H[b]–PxP–1 and H[b]–PxP–2 at 10.9 and 271 days, respectively.

4.1.0.3 Hypothesis c: Accepting planet candidates c and f

In Hypothesis H[c], we considered planet candidate f and c to be genuine planets along with planets (b, d and e), due to the support given by H[b]–PxP–2. In Hypothesis H[b], the distribution of potential orbital periods of the outermost planet is wider as compared to other predictions because DYNAMITE only has one neighbouring planet to use in constraining the orbit. However, as this predicted orbit also matches the period of planet candidate f closely, this prediction provides support for existence of planet candidate f under the assumption of lower planetary eccentricities. In this case, we found a new peak for H[c]–PxP–2, the most exterior likelihood peak, at 611 [549, 733] days, and we also recovered H[b]–PxP–1 as H[c]–PxP–1. The newly found peak, H[c]–PxP–2, hints at the presence of an additional planet in the system in the habitable zone of this system, near the outer edge. However, as this peak is again an extrapolation, lying beyond the longest-period confirmed planet, the predictions are here less constraining. Furthermore, we also found 3 regions between planets b and c, c and d, and e and f where the probability for an exoplanet to exist is comparatively low but non-zero.

4.1.0.4 Hypothesis d: Additional planets

Finally in Hypothesis H[d], we assumed (given the support from the previous hypotheses) that all the planet candidates (g, c and f) are real and added them in the system along with the planets (b, d and e). In Hypothesis H[c], the inner peak is better constrained due to two reasons: (1) the presence of an inner known planet and (2) the broken power law on the period followed by the first planet from the orbital period model. As the peak of the broken power law occurs at \(\sim\)​12 days, very near the orbital period of planet candidate g, this provides support for the genuine nature of planet candidate g. In our analysis of this planetary architecture, we retrieve the H[c]–PxP–2 peak at 611 days as H[d]–PxP–1, as well as four small regions in the period space with low injection probability interior to planet g and between planets b and c, c and d, and e and f. While there is no direct evidence for the existence of the predicted exoplanet H[d]–PxP–1 yet, our prediction is based on robust statistics from thousands of Kepler multi-planet systems. By construction, DYNAMITE predicts the most likely position in orbital period space for an additional planet to exist. When the system has low probability of finding an additional planet interior to the outermost system (i.e., when the system is dynamically packed like the system with all candidates treated as genuine planets), then the most likely orbit for an additional planet will be beyond the orbit of the outermost known planet. We also note that this prediction is only valid under the assumption that the eccentricities of the planets in the system are lower (see section 4.3) than previously reported.

For such an additional planet in the system, we find the most probable planet radius to be \(1.4^{+0.525}_{-0.365}R\oplus\), from fitting a Lognormal distribution to the most likely planet radii given the known planet minimum masses and the “volatile-rich" mass-radius relationship from [31]. From this, we find the most likely planet mass to be \(2.91^{+5.79}_{-1.88}M\oplus\). These values are estimated by DYNAMITE, given the minimum masses of the planets, the above mass-radius relationship, dynamical stability considerations with the six-planet system architecture, a uniform prior on the inclination assuming the planet is likely not transiting, and the distribution in planet radii assuming size clustering in a system [28].

4.2 Potential Equilibrium and Surface Temperatures↩︎

We estimated the equilibrium temperature of planet b to be \(641 \pm 23\) K, planet d to be \(379 \pm 14\) K, and planet e to be \(320 \pm 12\) K. The uncertainty in the equilibrium temperature of planets is the standard deviation of the temperature distribution calculated via a linear prior on the orbital period, with limits on the orbital period uncertainty provided by . With the help of the orbital period probability distribution from DYNAMITE, we estimated the equilibrium temperature of planet candidate g to be \(791 \pm 42\) K, planet candidate c to be \(486 \pm 26\) K, and planet candidate f to be \(241 \pm 17\) K. We calculated these temperatures assuming that only the planets b, d, and e are genuine, but even if a peak in the distribution is not occupied by a genuine planet, the other peaks remain unaffected.

Radiative equilibrium temperatures, however, may underestimate surface temperatures of bodies with significant atmospheres due to additional heat trapping by greenhouse gases. To account for this, we expanded our findings of equilibrium temperatures and calculated the possible surface temperatures of the planets and planet candidates in the system (see §3.3). These surface temperatures are calculated assuming thin atmospheres and low greenhouse warming, in this case a simple one-layer “leaky greenhouse” model. Based on this model, we estimate the possible surface equilibrium temperature of planet b to be \(689 \pm 27\) K, planet d to be \(407 \pm 16\) K, and planet e to be \(344 \pm 14\) K. These results are shown in Figure 6. For the planet candidates, we calculate the surface temperature of planet candidate g to be approximately \(849 \pm 52\) K, planet candidate c to be \(519 \pm 29\) K, and planet candidate f to be \(259 \pm 19\) K. These results can be found in Figure 4 and Figure 5. The uncertainties in the surface temperature are calculated in the same way as calculated for equilibrium temperature.

In Hypothesis H[d], where we assumed that all six planetary signals found by are genuine planets, we found a significant peak we call H[d]–PxP–1 at an orbital period of 611 days. When performing our temperature analysis on this planet candidate, we found that the equilibrium temperature of this planet would likely be \(196 \pm 9\) K. We also explored the potential greenhouse warming on this planet and estimate its surface temperature to be \(210 \pm 10\) K.

4.3 Eccentricity Analysis↩︎

In the following we report the fraction of stable and unstable orbital configuration iterations for four different hypotheses – system architecture hypotheses H[a] and H[d], each with both Normal and Lognormal eccentricity distributions. In essence, the normal distribution represents the solution reported by [7], while the Lognormal eccentricity distribution represents an agnostic approach based on typical planetary eccentricities. For each of our four cases, we report the combined mean eccentricity and its standard deviation. The results are also shown in Figure 13 in Appendix 12.

4.3.0.1 Hypothesis H[a] - Normal Distributions

Assuming that the system has only 3 planets – b, d, and e – and a Normal eccentricity distribution centered on the values reported by resulted in only \(1.2\%\) of iterations being stable configurations via the criteria established in §3.4. Planets in the stable configurations have a mean eccentricity of \(0.16 \pm 0.02\), whereas the planets in the unstable configurations have a mean eccentricity of \(0.24 \pm 0.04\).

4.3.0.2 Hypothesis H[a] - Lognormal Distributions

Under the above planetary architecture assumption but with a Lognormal eccentricity distribution from SysSim [32] instead, we observe that \(98.3\%\) of our N-body integrations yield stable configurations. Using this distribution, planets in the stable configurations have a mean eccentricity of \(0.05 \pm 0.02\) while planets in the unstable configurations have a mean eccentricity of \(0.14 \pm 0.03\).

4.3.0.3 Hypothesis H[d] - Normal Distributions

When we assumed that the system consists of 6 planets – g, b, c, d, e, and f – and that these planets have a Normal eccentricity distribution, we find that only one out of thousand iterations resulted in stable configurations. The eccentricities for each planet in the stable iteration were 0.069, 0.127, 0.161, 0.137, 0.15, and 0.064 for planets g, b, c, d, e, and f respectively. All the other 999 iterations resulted in unstable configurations.

4.3.0.4 Hypothesis H[d] - Lognormal Distributions

Lastly, for a 6 planet system with the planets having eccentricities as per Lognormal distribution resulted in \(98.9 \%\) of iterations being stable configurations. The mean eccentricity for all the planets in stable configurations was \(0.026 \pm 0.008\), while it was \(0.065 \pm 0.018\) for all the planets in unstable configurations.

5 Discussion - System Architecture↩︎

5.1 Sensitivity to Radial Velocity Solutions↩︎

In terms of planetary system-specific information, our analysis uses the results of planet detections from [6] and [7]. These detections are results of multi-parameter fits to the observed radial velocity modulations. To some extent, the results of such orbital fits depend on the underlying assumptions. A particular challenge in fits to radial velocity curves in multi-planet systems is the inter-dependence of the planets’ orbital parameters (for example, fitting for three planets may yield planets with different eccentricities than a fit assuming four planets). While these effects have been considered carefully in [7] – by increasing step-by-step the complexity of fits and relying on Bayesian Information Criteria to identify the preferred solutions – it is not trivial how such correlated uncertainties propagate into our analysis. Therefore, we will now discuss to what extent uncertainties in these studies and choices made by their authors may impact our results. The specific questions we will explore are as follows: How does the number of planets detected impact our results? and How do the radial velocity models (and their setup) impact our predictions?

Ultimately, our analysis is based on the veracity of the three confirmed planets (b, d, and e, see Table 2). In essence, these confirmed planets anchor the orbital architecture of the system. These three planets have been identified independently in two studies (b, and d) or detected in one state-of-the-art study at high significance level (e). Therefore, given the current knowledge on the system, the existence of these planets is safe to assume. As uncertainties remain on the orbital parameters of these planets, it is possible that their orbital periods, minimum mass, or eccentricities is somewhat different than currently inferred from the data. However, our results show that the planetary system architectures derived are insensitive to small- to moderate changes in the orbital parameters of these planets.

We note, however, that should one of these three planets turn out to be highly eccentric (not consistent with current orbital fits, [7]) or a spurious detection, then our analysis must be repeated and results may change. While such a development is unlikely in the face of the available data, it is prudent to discuss how our results may be affected by the existence and specific orbital parameters of the yet unconfirmed planet candidates (g, c, and f). Here, our analysis can provide a clear demonstration that the parameters or existence of these planet candidates does not have an effect on our prediction. This is the case because through the comprehensive set of hypotheses we tested (see §3.2 and Appendix 10), every combination of the genuine or spurious nature of these planet candidates has been explored. These repeated, independent assessments all converged to the same general planetary architecture, regardless of the existence or properties of the three candidate planets. Therefore, we conclude that our analysis is insensitive to the presence of the planet candidates, but would be affected if one of the three confirmed planets would, in fact, turn out to be a spurious detection.

5.2 Eccentricity Bounds↩︎

Our eccentricity analysis relies on the Normal and Lognormal eccentricity distributions parameters estimated by and [32], respectively. We tested the planetary system for the 3-planet and 6-planet architectures for both the above mentioned eccentricity distributions. In our analysis, we first integrated every planet’s orbit for \(\sim 5\times10^{6}\) orbits of the innermost planet and then predicted its long-term stability (\(\sim 5\times10^{9}\) innermost planet’s orbits) using spectral analysis. For the specific case of the system, our long term analysis thus corresponds to roughly 1.2 Gyr in the future. Any iteration that resulted in a collision, ejection, crossing of orbits during integration, or a high spectral fraction was considered as an unstable configuration. Stable pairs of planets with crossing orbits can exist, such as Neptune–Pluto (due to their 3:2 resonance) or potential planet pairs with relatively high mutual inclination and equal mass ratios [37], but in this system we are very unlikely to encounter such a situation (as these planets are likely non-resonant).

Our analysis suggests that even if the system has only 3 planets (b, d, and e), with the eccentricity parameters calculated by it is still highly unlikely that the system will remain dynamically stable over 1.2 billion years. If the system instead has 6 planets, then the likelihood of the system being dynamically stable over 1.2 billion years reduces even further to nearly zero. For this analysis, we used the minimum masses and the semi-major axes of the planets and planet candidates as calculated by . When we used the Lognormal parameters from [32] for our stability analysis, we found that more than 98% of iterations were stable for both system architectures tested. Thus, these results convey that the chances of the system being dynamically stable for at least a billion years with any planetary architecture are highly likely.

acknowledged that their eccentricity solutions were high and that factors like instrumental noise could be one of the reasons for such high eccentricities. Given that the system is thought to be older than 5.76 billion years, arguably the system is now – and has been – in a dynamically stable configuration. Our stability analysis, however, shows that stability is only possible if the planets indeed have lower eccentricities than those reported by . We estimate that if the system has three planets – b, d, and e – and the planet period and minimum mass parameters are close to those calculated by then each planet’s eccentricity would be of order \(0.05 \pm 0.02\). Similarly, if is a 6-planet system – with planets g, b, c, d, e, and f – then the most likely eccentricity for each planet would be of order \(0.026 \pm 0.008\).

We note that if the eccentricity vales are updated and better constrained in the future, then our analyses would have to be repeated. We also want to point out that our results from DYNAMITE analysis are valid under the assumption that the planets, and the candidates if they exits, will have lower eccentricity values than those estimated by in order for the planetary system to be dynamically stable.

6 Exploration of Planetary Temperatures↩︎

Here we explore the potential temperatures for the predicted exoplanets, considering the ranges of their predicted semi-major axes as well as a simple statistical prior on planetary Bond albedos. As little is known yet about the planets and predicted planets discussed here, the following discussion aims to provide only a simple, initial exploration of the potential temperatures of the planets in order to further our understanding of the potential natures of these worlds. The results presented in §4.2 are an extension of the Stefan-Boltzmann Law. We motivate and discuss our choice of Bond albedo prior in Appendix 9. Later, we also discuss the possibility of habitable planets in the system.

6.1 Equilibrium and Surface temperatures↩︎

For planets b and d, we estimated the equilibrium temperatures to be in the ranges of \(618-664\) K and \(365-393\) K, respectively. These ranges are \(\pm 1\sigma\) from the median values. Our results are fairly consistent with those of , stated above. For planet e, we estimate its equilibrium temperature to be in the range \(308-332\) K. We estimate the equilibrium temperature of planet candidate g to be in the range \(735-819\) K, planet candidate c to be \(455-509\) K, and planet candidate f to be \(224-258\) K. Our results for equilibrium temperature of planet candidate c is consistent with that of .

The most likely surface temperature of planet e lies in the range of \(330-358\) K, making it possible for planet e to harbor liquid water on its surface. The surface temperature of planet candidate f, estimated to be in the range of \(240-278\) K, is very similar to that of the Earth and makes it a candidate for harboring liquid water as well. If planet f has a mass greater than \(10 M_\oplus\), then based on the mass-radius relationship by [31], the planet will likely be gaseous rather than rocky, making it uninhabitable. We estimate that planet candidate g’s surface temperature will be in the range of \(788-882\) K, making its surface nearly as hot as that of Venus. Planet candidate H[d]–PxP–1 would have a similar equilibrium and surface temperature to Mars, so if it is rocky and massive enough to maintain a significant atmosphere to increase the greenhouse effect, it could also be potentially habitable. The surface temperatures of planets and candidates are shown in Figure 4 and Figure 5. In addition to applying our model to the planetary system – as a simple test of the new temperature prediction module – we also used it to estimate the surface temperature of Mars. Specifically, when Mars was hidden from the algorithm, based on the other inner Solar System planets, DYNAMITE not only "predicted" Mars’ orbit but also predicted a surface temperature to lie between \(184-220\) K. While running the DYNAMITE analysis, we excluded Jupiter from the inner Solar system and ran the simulation with Mercury, Venus, and Earth. In comparison, the measured temperature range is \(165 - 235\) K [38]. Thus, at least for planets without thick atmospheres and strong greenhouse effect, our approach provides a good initial estimate. A more detailed description of what we did is given in Appendix 11 and the results are shown in Figure 12.

Though our process for temperature estimates is simplistic, it is a reasonable approach given how little is known about the planets – detected or predicted – in the system. As planet properties are refined in the future, the single-layer leaky greenhouse model used in our study should be replaced with more comprehensive models.

6.2 Potential for Habitable Worlds↩︎

A key component of planetary surface habitability is the planet’s ability to host stable liquid water on its surface. Consequently, a simple constraint on the surface temperature of the planet might seem a valid first step towards characterizing a planet’s habitability. Surface temperatures of planet e and planet candidate f lie between 273 to 373 K, which makes them potentially suitable for hosting liquid water on their surfaces. However, even after this constraint, the planet might be uninhabitable due to several factors like insufficient surface gravity for atmospheric retention, or accretion of massive \(H_{2}-He\) envelope.

The minimum mass of planet candidate f, as determined by , is \(10.26^{+1.89}_{-1.47} M_{\oplus}\), which likely makes it a gaseous planet. Thus, planet candidate f, even after lying close to the inner edge of the habitable zone of the system, might not be habitable. [39] demonstrated that a necessary, but insufficient on its own, condition for a planet to be potentially habitable is that the equilibrium temperature of the planet must be below \(270\) K. They further show that in such a case, if the surface temperature of the planet remains below the critical temperature of water (\(T_{c} = 647\) K), the thermal emission of the planet will not exceed the greenhouse runaway limit of 300 \(Wm^{-2}\). In case the equilibrium temperature exceeds 270 K, either the planet will have insignificant water or its surface temperature will exceed 1400 K. In both cases, the planet will become uninhabitable, thus planet e (with an estimated equilibrium temperature of \(320 \pm 12\)) might not be potentially habitable.

As outlined by [40], the equilibrium temperature of a planet near the outer edge of the habitable zone lies in the range 175 – 200 K based on selective cloud coverage. Thus, we believe that the newly predicted planet H[d]–PxP–1 at 611 days would be a potentially habitable planet candidate if it has the right mass, has an equilibrium temperature between \(175 - 200\) K and can even lie in the habitable zone as determined by [22]. However, for this newly predicted planet to have the possibility of being habitable, it likely must have a significant amount of water, such that it can host liquid water for any temperature between 273 – 647 K, and it likely must be geologically active so that \(CO_{2}\) accumulates in the atmosphere as soon as the surface temperature drops below 273 K [39].

7 Biosignature Studies with Future Telescopes and Missions↩︎

The combination of three factors make it very likely that the system will play pivotal roles in the future of exoplanet exploration. First – as our study highlights – it hosts a rich planetary system with (possibly) six small planets. Second, is among the closest Sun-like stars and the closest system with such a large number of discovered/predicted exoplanets. Third, one – perhaps two – of its small planets would orbit in the habitable zone.

The very proximity of makes it an ideal target for near-future direct imaging studies: its habitable zone extends from 0.15 to 0.23 in projection (see §2.2). As a southern star (Declination \(-\)​43\(^\circ\)), it is an ideal target for both the European Extremely Large Telescope’s (E-ELT) PCS high-contrast, extreme adaptive optics imaging/spectrograph [41], as well as for the Giant Magellan Telescope’s high-contrast imaging systems (e.g., GMagAO-X, [42]). With its habitable zone planets seen at about 160 mas separation from , these instruments will likely be able to detect it in visible/near-infrared reflected/scattered light.

The detectability of the planets within the habitable zone will depend, fundamentally, on whether the given telescope’s inner working angle (IWA) will be smaller than the projected star–planet separation and, if so, whether the star–planet contrast will be detectable considering instrumental and astrophysical background signal. To explore which planets may be viable candidates for future instruments, in Table 3 we compare their maximum angular separation with the inner working angle typical to those instruments. For example, The LUVOIR mission concept’s ECLIPS coronagraph would provide continuous coverage at wavelengths ranging from 200 nm to 2.0  at a spectral resolving power of R=140 (in the visible) and R=70–200 (in the near-infrared). The inner working angle (IWA) of ECLIPS would be 4 \(\lambda/D\) (in the ultraviolet) and 3.5  \(\lambda\)/D (in the visible and near-infrared). Coronagraphic imagers on HabEx [3] and on large ground-based telescopes (e.g., [42] and [41]) will also allow very similar contrast-limited performance, but coupled to different diameter mirrors. Assuming a 4\(\lambda/D\) limit, for each planet we calculate the required minimum mirror diameter at 0.4  and at 1.5 \(\micron{}\). These minimum diameters are shown in Table 3. We determine the detectability (last three columns) by comparing the inner working angle for the given instrument to the maximum projected separation of the planet. Our assessment suggests that planets close to or within the habitable zone (f and H[d]-PxP-1) are excellent candidates for both space- and ground-based high-contrast imaging.

Furthermore, future thermal infrared imagers at the VLT – perhaps building on the NEAR concept [43], [44] – may also be able to detect some of the planets close to the habitable zone of . In fact, has also been identified among the highest-priority targets for the high-resolution thermal infrared imaging with the E-ELT’s METIS instrument [45].

Note: The required minimum mirror diameter assumes the planet is observed at maximum elongation and assumes an inner working angle of 4.0 \(\lambda\)

/D.

Finally, future interferometric missions focusing on habitable planets – such as the LIFE mission concept [46] – will likely focus on due to its relatively easy-to-observe yet rich planetary system and likely presence of potentially habitable planets.

8 Summary↩︎

We analyzed system for its multiplicity, the planetary architecture and orbital parameters of the planets, the equilibrium and surface temperature of the planets and potential candidates, and the dynamical stability of the system based on varying planetary eccentricities.

Starting from the assumption that the system consisted of only three planets – b, d, and e – the DYNAMITE analysis predicted three more potential candidates – g, c, and f – which had also been reported in earlier studies but could not be confirmed. No information about the potential candidates was used as an input while running the DYNAMITE analysis. The planetary parameters for these planet candidates were in good agreement with those reported in previous works ( and ). Subsequently, when all the 6 planets were used as an input, DYNAMITE predicted one additional planet in the habitable zone of the system (P = 611 d) with a likely radius and mass of \(1.25-5.16 R_\oplus\) and \(1.03-8.7 M_\oplus\), respectively.

Our analysis of the equilibrium and surface temperatures of the planets was based on an albedo prior derived from the celestial bodies in the Solar system and the assumption of thin atmosphere. We find that though planet candidate f lies in habitable zone, its likely gaseous nature will make it uninhabitable. We find that the newly predicted planet (H[d]-PxP-1) might be habitable.

From our eccentricity analysis, we find that if is a three-planet system with planets b, d, and e, then the combined mean eccentricity for the system to be stable is \(\sim0.05\). If the system is a six-planet system instead, then the combined mean eccentricity for the system to be stable is of order \(\sim 0.026\). In either case, we find that the eccentricity of each planet should be significantly lower than the value fitted to the RV data, as also proposed by [7].

As the planetary system’s stability necessitates a lower-than-reported eccentricity for the planets, our analysis is based on this assumption. If better constraints on the eccentricities become available in future, then our analysis could be repeated again with the updated values.

Our assessment of the system provides support for a seven-planet architecture with one planet candidate in the habitable zone of the system. Given its proximity and – in projection – large habitable zone as well as its rich inner planetary system, the nearby is likely to further ascend on the lists of planetary systems that are the most promising for exploration via direct imaging and for the search for life.

The results reported herein benefited from collaborations and/or information exchange within the program “Alien Earths” (supported by the National Aeronautics and Space Administration under Agreement No. 80NSSC21K0593) for NASA’s Nexus for Exoplanet System Science (NExSS) research coordination network sponsored by NASA’s Science Mission Directorate. We acknowledge use of the software packages NumPy [47], SciPy [48], and Matplotlib [49]. This paper includes data collected by the Kepler mission. Funding for the Kepler mission is provided by the NASA Science Mission Directorate. We obtain the data set from the NASA Exoplanet Archive [[50]]³ This dataset or service is made available by the NASA Exoplanet Science Institute at IPAC, which is operated by the California Institute of Technology under contract with the National Aeronautics and Space Administration. The citations in this paper have made use of NASA’s Astrophysics Data System Bibliographic Services. This research uses data from the Solar System Dynamics Small-Body Database, which is hosted by the Jet Propulsion Laboratory and sponsored by NASA under Contract NAS7-030010. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

9 Selection of Bond Albedos for planets↩︎

Our estimate of the equilibrium temperatures is informed and influenced by a prior Bond albedo distribution. Though the albedo distribution of exoplanets is not yet determined, we can provide a reasonable prior on the Bond albedos considering current observational constraints on exoplanets and data on Solar System bodies. [51] devised a statistical method to estimate the average albedo of an ensemble of exoplanets. They stacked multiple light curves of sub-Saturn planet candidates from the Kepler Object of Interest table, between \(2R_\oplus - 6R_\oplus\), and measured the average depth of the light curve, from which they estimated the average albedo of the ensemble. They considered two different cases: including the planet Kepler-10b and excluding it, such that Kepler-10b does not dominate the results due to its unusually deep secondary eclipse.

In a follow-up study, [52] expanded the sample size of the ensemble and divided the ensemble into three sub categories: \(1-2R_\oplus\), \(2- 4R_\oplus\), and \(4-6R_\oplus\). They then subdivided each class further into two different lists based on if the false-positive probability was less than \(1\%\). In addition, the analysis provided two cases, one case assuming that the planet completely redistributed the heat and the second assuming instantaneous re-radiation of heat by the planet. [53] used a different statistical approach for estimating the average albedo of Terran and Neptunian planet ensembles. Instead of using light curves and planet candidates, multiple phase curves of confirmed planets were stacked, in order to estimate the albedos of the ensemble from the resulting average phase curve. The resulting albedos from all these works can be found in Table 4.

Solar System bodies offer greater amount of information and with higher precision than current exoplanet studies, albeit for objects that may differ in nature from the exoplanets such as those found in the system. Nevertheless, we consider the approximate range and distribution of the Bond albedos for Solar System bodies as a general context for our prior. Our discussion starts with asteroids and then continues with moons, rocky planets, and giant planets.

[54] calculated phase integral values of Solar System asteroids as a function of their geometric albedos and classified all the asteroids in three main classes: low albedo, medium albedo, and high albedo asteroids. Bond albedo is defined as \(A_{B} = A_g q\), where \(A_g\) is the geometric albedo of the body and \(q\) is the phase integral. Therefore, the Bond albedos can be calculated for numerous asteroids listed in the dataset provided by the Jet Propulsion Laboratory through its Solar System Dynamics Small-Body Database⁴. On using the average phase integral values of the three asteroid classes, we find two discontinuities in the Bond Albedo distribution of asteroids. To remove these discontinuities in the Bond Albedo distribution, instead of using the average phase integral for different classes, we use the linear fit equation, q = 0.359 + 0.47p [54] to calculate the phase integral, from which we calculate the Bond albedo of each asteroid. Besides asteroids, we also collected Bond albedos of the eight planets and many satellites in our Solar System, which can be found in Table 5.

We compared the Bond albedos of various bodies in our Solar System with the constraints placed on Bond albedos by various studies (as presented in Table 4). Venus might be an outlier as it is not consistent with the results from the studies by [51], [53], and [52]. Earth’s albedo is consistent with [53], but it is only consistent with the results of [52] under the assumptions of complete heat redistribution and include Kepler-10b, along with those of [51] which exclude Kepler-10b. Mars’ albedo is almost consistent with [53]. Lastly, the albedos of Uranus and Neptune are consistent with both [53] and [52], but are only consistent with those results of [51] which exclude Kepler-10b. When using DYNAMITE to predict the size of the planets around , we find that the radii of planets lie roughly in \(1-2R_\oplus\) range, and likely to be Terran or Neptunian in nature. As a result, we chose the upper limits placed on the Bond albedo by [53], as it is a broad range which meaningfully captures maximum data points and matches the likely natures of the planets. With these upper limits, nearly 69% and 99.8% of albedos drawn from our adopted Gaussian distribution fall below \(A_{B} = 0.35\) and \(A_{B} = 0.63\). Thus, we conclude that the Gaussian distribution of albedos is a reasonable prior for our initial, exploratory estimate of the equilibrium temperatures of the predicted planets in the system.

From these albedos, we calculated the equilibrium temperature of various bodies using two simple assumptions: the celestial bodies radiate as a blackbody, and they are in radiative equilibrium with their surroundings. We also used the Leaky Greenhouse Model (see §4.2) to calculate equilibrium temperature of those bodies which are larger than Mercury in size. We show the Bond albedos and Equilibrium Temperatures of all these bodies in Figure 3, as well as a histogram showing the Bond albedo distribution of various bodies. For bodies of the size of Mercury or larger, the median Bond albedo is 0.3 and the standard deviation is 0.178. On the right hand side of this figure, we show our adopted Gaussian Distribution of Bond Albedos centered at \(\mu = 0.3\) with \(\sigma = 0.1\).

Note

: Albedos of Venus, Earth, Mars, Saturn, Moon, Io, Callisto, Ganymede, Eris, Borelly and 11 Parthenope did not have any uncertainity attached.

10 Additional Hypotheses↩︎

The four hypotheses presented in §3 do not capture all the possible configurations of confirmed planets and unconfirmed planet candidates in the system. Therefore, we present here all the remaining possible combinations of the three confirmed planets and three planet candidates in the system, which is summarized in Table 6. As before, the starting point of our DYNAMITE analysis for all scenarios are the three confirmed planets. When we use DYNAMITE with planets b, d, and e, there are three strong likelihood predictions near the locations of the unconfirmed planet candidates g, c, and f. We then progress through each of the possible planet candidates, assuming one additional candidate is a genuine planet. Finally, we show what the DYNAMITE predictions give for the full 6-planet system.

11 DYNAMITE analysis on Inner Solar System and Temperature Estimation for Mars↩︎

In addition to analysing the planetary system, we also ran DYNAMITE on the Inner Solar System while deliberately excluding Mars from the system. We found three peaks in the orbital period space distribution plot out of which the inner two peaks are significantly smaller than the outer peak. The first peak (at 43.5 days) lies interior to Mercury, the second peak (at 141 days) lies between Mercury and Earth, while the last peak (at 673 days) lies exterior to Earth. The outermost peak is less constrained as DYNAMITE tries to match the Kepler period ratios by adding planets exterior to Earth. The outermost peak is the only significant peak and its location matches to that of Mars’ orbital period (\(\sim 687\) days – an excellent 2% match). This gives another line of evidence that the model accurately retrieves the intentionally removed planet.

Using the methods explained earlier in §3.3, we calculated the likelihood surface temperature of Mars to be \(202 \pm 18\) K (\(\pm\)​2\(\sigma\) range approx. 166 K – 238 K). Keeping in mind that Mars has a fairly complex climate and shows weak Greenhouse effect, our calculated surface temperature for Mars is an excellent match with the measured temperature range of the Martian surface \(165-235\) K [38]. Therefore, the predicted surface temperature is in excellent agreement with the actual surface temperature range. Given the typical limitations on the atmospheric characterization of exoplanets, our model can provide useful initial exploratory estimates of the surface temperature ranges that may be typical to them.

12 Dynamical Stability analysis of the e Eridani System↩︎

Here we present the results of our dynamical stability analysis of planetary system as a function planetary eccentricities. We analyzed the three-planet and six-planet system with: (1) Normal eccentricity distribution following the RV solutions from [7] and (2) Lognormal eccentricity distributions from [32] (see Section 3.4 for more details). We calculated the combined average eccentricity for each of the four cases. The results for each run are shown in Figure 13.

References↩︎

---

1. https://github.com/JeremyDietrich/dynamite↩︎

2. https://exoplanetarchive.ipac.caltech.edu/↩︎

3. Accessed on 2022-01-23 at 09:00, returning 1 row.)↩︎

4. https://ssd.jpl.nasa.gov/sbdb_query.cgi↩︎