GJ 357 d: Potentially Habitable World or Agent of Chaos?

As described in Section 1, the GJ 357 system consists of three known planets orbiting a low-mass star. For the analysis in this work, we adopt the stellar and planetary parameters of Luque et al. (2019). The host star has a spectral classification of M2.5V, with a mass of M_⋆ = 0.342 M_⊙, a radius of R_⋆ = 0.337 R_⊙, an effective temperature of T_eff = 3505 K, and a luminosity of L_⋆ = 0.01591 L_⊙. These properties allow us to calculate the HZ of the system, including the conservative HZ (CHZ) and the optimistic extension to the HZ (OHZ) based upon the assumption that Venus and Mars had surface liquid water in their past, described in detail by Kane et al. (2016). We calculate distance ranges of 0.131–0.254 au and 0.103–0.268 au for the CHZ and OHZ, respectively. The extent of the HZ and the orbits of the known planets are shown in Figure 1, where the CHZ is shown in light green and the OHZ is shown in dark green, and the semimajor axes of the planetary orbits are 0.35 au, 0.061 au, and 0.204 au for the b, c, and d planets, respectively. It is worth noting that Jenkins et al. (2019) refer to the innermost planets of the system as "c" and "b" in order of increasing semimajor axis, whereas Luque et al. (2019) refer to those same planets as "b" and "c." Here, we adopt the naming convention of Luque et al. (2019), and thus the three planets are referred to as "b," "c," and "d" in order of increasing semimajor axis.

Zoom In Zoom Out Reset image size

There are various components of the GJ 357 system that remain unconstrained. The planetary orbits shown in Figure 1 are assumed to be circular, and indeed statistical studies have found that smaller planets tend to have relatively low eccentricities (Kane et al. 2012; Van Eylen & Albrecht 2015). Luque et al. (2019) considered eccentric orbits in the preliminary analysis of the RV data, but assumed circular orbits for their final model due to comparable fits between eccentric and noneccentric cases, and the computational burden of including eccentricity as a free parameter. Jenkins et al. (2019) also explored nonzero eccentricities, which included a 1σ eccentricity for planet d consistent with ∼0.1, but fixed circular orbits in the final analysis. Thus, an eccentric orbit for planet d is allowable with the present data for the system.

The true architecture of the system is only known to the extent that the sensitivity of the observational data allows, including whether there may be further planets within the system. For example, an additional planet of Earth-mass that lies within the HZ would be challenging to detect with the present data set. We calculate that an Earth-mass planet at the inner and outer edges of the OHZ would have RV semiamplitudes of 0.47 m s⁻¹ and 0.29 m s⁻¹, respectively, which fall below the 2 m s⁻¹ rms precision of the utilized spectrographs. We therefore find that an additional Earth-mass planet in the HZ would likely remain undetectable with the current RV data.

According to Luque et al. (2019), planets c and d have minimum masses of 3.4 and 6.1 Earth masses, respectively. Since planets c and d are not presently known to transit, the true mass of the planets may be significantly higher than the minimum masses. A dynamical analysis of the system performed by Luque et al. (2019) did not place significant constraints on the inclination of planets c and d. However, the dynamical analysis performed by Jenkins et al. (2019), whose RV analysis derived slightly higher planetary masses than those found by Luque et al. (2019), determined that the mass of planet d lies in the range of 7.2–11.2 Earth masses. The resolution on whether or not planet d transits the host star is an important component of understanding the true nature of this HZ planet.

The analysis performed by Jenkins et al. (2019) and Luque et al. (2019) used Sector 8 of TESS photometry, which occurred during the TESS Prime Mission. Since then, GJ 357 has been observed again during sectors 35 and 62. The time-series photometry observed by TESS was obtained through the Mikulski Archive for Space Telescopes (MAST Team 2021). We utilize the 2 minutes presearch data conditioning simple aperture photometry (PDCSAP), which was processed by the Science Processing Operations Center pipeline (Jenkins et al. 2016). Figure 2 shows the light curve of GJ 357 for all three TESS sectors, where the transits of GJ 357 b can be plainly seen in the light curve every ∼4 days. The transit depth is 1095 ppm, which translates to a radius for planet b of R_p = 1.217 R_⊕. Given that the average rms scatter for the shown three TESS sectors is 620 ppm, and that planets c and d are more massive than planet b, transits of the outer two planets should be visible in the data if they occur. Indeed, planets of size 1.5 R_⊕ and 2.0 R_⊕ would produce transit depths of 1665 ppm and 2960 ppm, respectively. Assuming the lowest of these radii values, transits of planets c and d are ruled out at 2.7σ per TESS measurement, which is equivalent to that expected for a grazing transit. Based on the system properties described in Section 2.1, we estimate central transit durations of ∼2.0 and ∼3.5 hr for planets c and d, respectively. Thus, central transits for planets c and d are ruled out at a significance of ∼21σ and ∼28σ, respectively.

Zoom In Zoom Out Reset image size

Given the expected transit depths and photometric precision, the three TESS Sectors shown in Figure 2 are sufficient to confirm that planet c (whose orbital period is 9.12 days) does not transit. Planet d has an orbital period of 55.661 days, and so its alignment with a particular TESS sector is less obvious than for planet c. However, the RV orbit shows that an inferior conjunction passage missed Sector 8 entirely, as noted by both Jenkins et al. (2019) and Luque et al. (2019). Using the planet d orbital ephemeris of Luque et al. (2019), we calculate BJD times of inferior conjunction of 2459272.34 and 2459995.93, which occur during the sector 35 and 62 observing windows, respectively. These inferior conjunction times are indicated by the vertical dashed lines shown in Figure 2. The inferior conjunction for sector 35 occurs directly after the telemetry data gap, and reveals no evidence for a transit. The inferior conjunction for sector 62 also falls within the TESS photometry, and in fact falls at the same time as a planet b transit, which would have resulted in a syzygy transit event had planet d transited (Luger et al. 2017; Veras & Breedt 2017). The lack of observable signature at either location verifies that planet d does not transit the host star.

As noted in Section 2.1, the confirmation that planet d does not transit means that the planetary mass may be significantly higher than the minimum mass of 6.1 M_⊕. There have been numerous derivations of mass–radius relationships for exoplanets that estimate the upper limits for a terrestrial body (Dressing et al. 2015; Rogers 2015; Chen & Kipping 2017). The empirical relationship between planet mass and radius derived by Chen & Kipping (2017) detected a transition from terrestrial into "Neptunian" planets with a greater volatile inventory at a boundary of ∼2 M_⊕. Without a radius measurement, there is a great deal of degeneracy regarding the bulk properties of a planet that possibly lies within the terrestrial regime, even if the true mass of the planet is known (Valencia et al. 2007; Dorn et al. 2015; Zeng et al. 2016). The subsolar metallicity of GJ 357 (Luque et al. 2019), combined with the relatively high mass of planet d, suggests the planet lies outside the nominal rocky planet zone (Unterborn et al. 2023). Furthermore, Kopparapu et al. (2014), who assume that planets with masses larger than 5 M_⊕ are not rocky, found an increasing rate for the outgoing longwave radiation with planet mass due to the smaller atmospheric column depth, decreasing greenhouse warming and increasing the width of the HZ at the inner edge. With all of these considerations in mind, it is difficult to state with any certainty what the true nature of planet d is, but the minimum planetary mass allows for a large parameter space where a habitable scenario is increasingly unlikely.

We adopt the methodology described by Kane (2019), Kane et al. (2021c), in which the Mercury Integrator Package (Chambers 1999) was applied with a hybrid symplectic/Bulirsch-Stoer integrator with a Jacobi coordinate system (Wisdom & Holman 1991; Wisdom 2006). Each simulation was integrated for 10⁶ yr, equivalent to ∼6.5 × 10⁶ orbits of planet d, and with a time step of 0.1 day to ensure adequate resolution of planet–planet encounters that involve the innermost planet. We used the stellar and planetary properties provided by Luque et al. (2019) and the HZ boundaries calculated in Section 2.1. We tested the orbital stability within the HZ by injecting an Earth-mass planet in a circular orbit that is coplanar with planets b and c, the latter of which is assumed to have a near edge on orbit that allows the minimum planet mass to be adopted as an approximation of the true mass. The new planet was injected at semimajor axes within the range 0.1–0.27 au and in steps of 0.001, encompassing the full OHZ range of 0.103–0.268 au. Additionally, each semimajor axis step incorporated initial evenly spaced mean anomalies of 60°, 180°, and 300° for the injected planet.

Because the orbit of planet d is poorly constrained, our simulations were conducted for three specific orbital configurations of planet d. First, we considered the case of planet d being near coplanar with the other planets such that, like planet c, the minimum planet mass (6.1 M_⊕) is a reasonable approximation of the true mass. Second, we considered the case of the planet d orbit being significantly inclined with respect to the orbital plane of the other planets, producing a planet mass of 10.0 M_⊕, which is within the dynamical limits found by Jenkins et al. (2019). Third, we considered the case where the orbital inclination and planet mass of planet d is the same as the second case, but now with a slight eccentricity of 0.1. The combination of these three cases resulted in ∼1500 simulations in total.

The outcome for each of the simulations described in Section 3.1 was assessed based on the survival of the injected planet, where nonsurvival can mean the planet was either ejected from the system or lost to the gravitational well of the host star. Note that the Mercury Integrator Package assumes point masses, and so planet–planet collisions, though possible, are not considered in the simulations. The results of the simulations are summarized in Figure 3, which contains three panels that represent the results for each of the three orbital configuration cases for planet d mass and eccentricity, shown near the top of each panel. Each panel shows as a solid line the survival time of the injected planet (as a percentage of the full 10⁶ yr integration) as a function of the semimajor axis of the injected planet. The HZ is depicted as for Figure 1, with the CHZ shown in light green and the OHZ shown in dark green. The semimajor axis of planet d is indicated by the vertical dashed line. We detected no cases in which any of the three known planets did not survive the simulations, largely due to their substantial mass compared with the injected planet.

Zoom In Zoom Out Reset image size

For all three cases, there is an island of stability around the semimajor axis of planet d, allowing for the possibility of Trojan planets in similar orbits (Páez & Efthymiopoulos 2015). Remarkably, such Trojan planetary orbits can maintain long-term stability (Cresswell & Nelson 2009; Schwarz et al. 2009), although eccentricity of the primary planet reduces this stable region (Dvorak et al. 2004), as seen in the bottom panel of Figure 3. For the first case (M_p,d = 6.1 M_⊕ and e_d = 0.0; top panel), the presence of planet d clears a substantial region around the orbit, and 18% of the HZ is rendered unstable. For the second case (M_p,d = 10.0 M_⊕ and e_d = 0.0; middle panel), the instability regions surrounding the orbit of planet d slightly expand to occupy 21% of the HZ. For the third case (M_p,d = 10.0 M_⊕ and e_d = 0.1; bottom panel), the introduction of a relatively small eccentricity to the planet d orbit greatly increases the instability with the HZ, resulting in 60% of the HZ being unstable for the injected planet. Instability of the injected planet primarily arises though gravitational perturbations from the known planets, particularly at mea motion resonance (MMR) locations, that increase the eccentricity of the injected planet. Such eccentricity increases are often lead to more frequent perturbations that culminate in the ejection of the planet from the system. These simulation results show that, for the planet mass range explored, increasing the eccentricity of planet d has a larger effect on the instability within the HZ than increasing the planet mass. However, regions of the HZ where the injected planet remains in the system does not guarantee that the planet's orbit is conducive toward potential habitability.

If the injected terrestrial planet in our simulations survived the 10⁶ yr integration time, there may remain orbital consequences from interacting with the other planets in the system. In general, compact planetary architectures benefit from stability enabled by the planets' relatively small Hill radii, since that scales linearly with semimajor axis. The results shown in Section 3.2 demonstrate that increasing planet mass, and therefore Hill radius, gradually increases the region of instability surrounding the planet. However, increasing eccentricity has a far greater effect on system dynamics through conservation of angular momentum, and even stable configurations can inherit significant eccentricity evolution cycles.

Another feature shown in Figure 3 is the maximum eccentricity achieved by the injected planet, indicated as black dots in each of the panels. These are shown for the cases where the planet survives the full integration time. The maximum eccentricity values shown for the first two architecture cases (top and middle panels) show that the injected planet eccentricities remain low when the initial conditions of planet d assume a circular orbit. Exceptions to this include slight eccentricity increases at locations of MMR, such as 0.129 au and 0.155 au, corresponding to 2:1 and 3:2 MMR with planet d, respectively. More significant exceptions are those close to the instability regions surrounding planet d, where planetary orbits lie at the edge of chaotic instability. An example of this is shown in Figure 4, which provides the eccentricity evolution for all four planets in the system in the case where M_p,d = 6.1 M_⊕, e_d = 0.0, and the injected planet has a semimajor axis of 0.225 au. The data shown in Figure 4 are the first 10⁶ yr of an extended 10⁷ yr integration conducted to explore the longer-term stability for this particular architecture. The interaction with planet d quickly raises the eccentricity of the injected planet into a quasichaotic state where it remains up until ∼0.5 × 10⁶ yr, even influencing the eccentricity evolution of planets b and c. Beyond ∼0.5 × 10⁶ yr, the injected planet enters into a stable periodic exchange of angular momentum with planet d. Though the injected planet survives the 10⁷ yr integration, there is no guarantee that the system will retain stability beyond the simulated period.

Zoom In Zoom Out Reset image size

The third architecture case (bottom panel of Figure 3) shows a much larger excitation of the injected planet eccentricities, where the maximum eccentricity increases with increasing semimajor axis and remains at ∼0.2 for the majority of 100% survival simulations. We conducted 10⁷ yr simulations for several initial semimajor axis values of the injected planet and found that the planet either did not survive the full simulation or showed increasing signs of chaotic behavior until the end of the simulation. Thus, many of the injected planets in the third architecture case are unlikely to maintain long-term stability beyond the time frame of the simulations shown in Figure 3, particularly at the locations of MMR, due to the chaotic nature of the induced eccentric orbits.