Dynamics of a Probable Earth-mass Planet in the GJ 832 System

S. Satyal; J. Griffith; Z. E. Musielak

doi:10.3847/1538-4357/aa80e2

1. Introduction

The recent discovery of exoplanets has shown that the multi-planetary systems in compact orbits seem to be common in the Milky- Way Galaxy. For example, Kepler 186 (Quintana et al. 2014) is a five-planet system in which the farthest planet from the star, Kepler 186f, is located at 0.432(+0.171, −0.053) au (Torres et al. 2015) and is within the habitable zone of the host star. The nearest planet from the star, Kepler 186b, is at 0.0378 au and orbits the star every 3.88 days. Another such planetary system is Gliese 581 (Mayor et al. 2009), which is known to host three planets along with the two others that are not yet confirmed. The three known planets, Gliese 581b, c, and e, orbit the star within 0.07 au. Also, many binary star systems have been discovered that are known to host multiple circumprimary planets in S-type orbits, for example 55 Cancri b, c, d, e, and f (Fischer et al. 2008). With the addition of a third planet in the Kepler 47 system (J. A. Orosz et al. 2016, in preparation), it became the first multi-planet circumbinary system and opened a new chapter for us to understand the planet formation processes and the dynamical compactness of planetary orbits.

The GJ 832 planetary system (Wittenmyer et al. 2014) is another multi-planet system that is currently known to host two planets around an M-dwarf star and is located at a relatively close distance of 16.1 light years from Earth. GJ 832c (inner planet) orbits its star at a distance 0.163 ± 0.006 au away and is potentially a rocky planet with a mass 5.4 ± 0.95 M_⊕. This planet is located in the inner boundary of the habitable zone, but it is not expected to be habitable mostly due to its close proximity to the star and its possibility of having a dense atmosphere (Wittenmyer et al. 2014). Orbiting the same star, a distant planet GJ 832b (outer planet) was discovered by Bailey et al. (2009); which is a long-period, 3657 ± 104 day giant planet at 3.56 ± 0.28 au, with a mass (m sin i) of 0.64 ± 0.06 M_J. The planets orbit a main-sequence dwarf star, GJ 832, of spectral type M1.5V (Jenkins et al. 2006), with a mass of 0.45 ± 0.05 M_⊙ (Bonfils et al. 2013) and a temperature of 3472 K (Casagrande et al. 2008). This is a fairly old system. Guinan et al. (2016) estimated the age of the star to be 6 ± 1.5 Gyr by using their X-ray Activity-Age relationship. We have calculated the size of the stellar habitable zone by using the formula provided by Kopparapu (2013). Then, the orbital stability of additional planets with masses ≤15 M_⊕ is investigated within and outside the boundaries of this habitable zone.

The main goal of this article is to explore the gravitational effect of the outer planet on the orbital stability of the inner planet as well as on planets with masses ≤15 M_⊕ injected between the two known planets. In addition, the long-term stability and orbital configurations of the inner and the injected planets, with a concentration on the time evolution of their semimajor axis (a_pl), eccentricity (e_pl), and inclination (i_pl), are studied in the a_pl − e_pl and a_pl − i_pl phase-spaces, and time-series 2D plots.

We also used the integrated data from the time evolution of orbital parameters to generate synthetic radial velocity (RV) curves of the known planets as well as the added planets with masses varying from 1 M_⊕ to 15 M_⊕ in the system. Moreover, based on the maximum amplitude of the RV curve obtained from the observation of the inner planet, the approximate mass and distance from the star for the potential Earth-mass planet were computed using the RV signature of the Keplerian motion.

This paper is outlined as follows. In Section 2, we describe our numerical simulations, the results are presented and discussed in Section 3, and the paper is concluded with a brief summary of our main results in Section 4.

2. Numerical Simulations

The best-fit orbital parameters of the GJ 832 system as obtained from the original discovery papers (Bailey et al. 2009; Wittenmyer et al. 2014) are given in Table 1 and for the cases when the parameters are unknown, they are either set close to zero or set to a range with a fixed step size. Both of the known planets in the GJ 832 system were detected by the RV technique from which the orbital parameters were extracted by using the best-fit Least-Squares Keplerian Orbital Solutions. We used these parameters as the initial conditions for starting our numerical simulations. The best-fit orbital solutions include uncertainties in their respective parameters. For example, the semimajor axes of the inner and outer planets with uncertainty are 0.163 ± 0.006 au and 3.56 ± 0.28 au, and the masses (m sin i) are 5.4 ± 0.95 M_⊕ and 216 ± 28 M_⊕, respectively. In some of our phase-space analyses, we have considered the maximum uncertainty values, especially in the planets' eccentricities, semimajor axes, and masses. This has allowed us to see the gravitational effects on the long-term orbital stability of the inner planet as well as the injected planet for some extreme values of their orbital paramters.

Table 1. Best-fit Orbital Parameters of the GJ 832 System Obtained from Wittenmyer et al. (2014) and Bailey et al. (2009)

Parameters	GJ 832b	GJ 832c
m sin i (M_⊕)	216 [188, 245]	5.4 [4.45, 6.35]
Semimajor axis (a)	3.56 au [3.28, 3.84]	0.163 au [0.157, 0.169]
Eccentricity (e)	008 [0.02, 0.1]	018 [0.05, 0.31]
Inclination (i)	(0–90)°	(0–90)°
Longitude of the ascending node (Ω)	(10⁻⁵)°	(10⁻⁵)°
Argument of the periapsis (ω)	246° [224, 268]	100 [323, 57]
Mean anomaly (μ)	307° [285, 330]	165° [112, 218]
Period (P) (days)	3657 [3553, 3761]	35.68 [35.65, 35.71]

Note. Mass of the star = 0.45 ± 0.05 M_⊙. The initial values for inclination and longitude of the ascending node of the planets are set by us.

Download table as: ASCII Typeset image

We have considered the motion of the planets of masses, m_pl around the central star in the general elliptical as well as the circular cases. To calculate the planets' and the star's initial conditions (ICs) in terms of their position and velocity and start the integration processes, we used their best-fit orbital elements: semimajor axis (a), eccentricity (e), inclination (i), argument of periapsis (ω), ascending node (Ω), and mean anomaly (M), which were obtained from RV measurements (Wittenmyer et al. 2014). The initial orbital inclination of the inner and the injected Earth-mass planets are taken relative to the orbital plane of the star and the outer planet. Thus, any inclination we mention during our investigation is relative to the star–outer planet plane. To check the stability of the system for other limiting cases, we have used the upper limit of the best-fit orbital elements. Note that, in this paper, the stability of the system is defined in terms of its lifetime, which is based on the time a planet survives the total simulation time, collides with other bodies, or gets ejected from the system during the orbital integration period. The stability is also studied in term of the maximum eccentricity, e_max attained by the planetary orbits during their evolution processes. A system is considered stable when the integrated bodies survive the total simulation time and their e_max deviate least from the initial value; else, the system is considered unstable.

Using the orbital integration package MERCURY (Chambers & Migliorini 1997; Chambers 1999), the built-in Hybrid algorithm was used to integrate the orbits of the system in astro-centric coordinates. MERCURY was effective in monitoring the ejection or collision of the inner and the injected planets due to a close encounter with the star or the outer planet. While integrating the orbits, a time step of = 10⁻³ year/step was considered to obtain high precision data and minimize the error accumulation. The change in total energy and total angular momentum was calculated at each time step, which fell within the range of 10⁻¹⁶ to 10⁻¹³, respectively, during the total integration period of 10 Myr, and the range of 10⁻¹⁰ to 10⁻¹² during 1 Gyr. The data sampling (DSP) was done per day and per year for shorter integration periods and at every 100 kyr for a billion-year integration period. The lifetime maps and the maximum eccentricity (e_max) maps are generated for multiple (up to 14,400) initial conditions in a_pl, e_pl, and i_pl phase spaces, and they are simulated for 10 Myr. The billion-year simulations are performed for low-resolution phase-space maps and for some selected single initial conditions.

3. Results and Discussion

3.1. Dynamics of GJ 832c

The inner planet, GJ 832c, does not have well constrained orbital parameters (Table 1), including its orbital inclination (i_pl) and the longitude of ascending node (Ω_pl), which do not have the best-fit values. Therefore, to set up our initial conditions for the simulations, we have set Ω_pl to (10⁻⁵)^o and i_pl to a value between near co-planar (10⁻⁵)^o and 90^o. Other orbital parameters are set at their best-fit nominal values given in Table 1. Similarly, to set up the initial conditions for the outer planet, GJ 832b, its best-fit nominal values are considered when available; otherwise they are set close to zero. We have also performed simulations where the initial conditions are set to upper values of the uncertainty limit. The lower uncertainty values of the orbital parameters would have less of an effect on the stability of the phase-spaces between the known planets, which is our major region of interest for an additional planet; hence, no simulations are run for such cases. The ω_pl, Ω_pl, and M_pl are considered slow moving angles and have the least effect on the orbital stability. The influence of these angles is more significant in the resonant angle studies, which determines the libration or circulation of the phase angles (Murray & Dermott 1999). Also, see Satyal et al. (2014) for the analysis of these angles in the study of the chaotic dynamics of the planet in the HD 196885 binary system.

To investigate the orbital inclination of the inner planet, its orbits are integrated with 14,400 initial conditions (ICs) in varying i_pl − a_pl phase-space, and 8000 ICs in e_pl − a_pl phase-space. To perform our simulations, we only considered prograde orbits, where i_pl is sampled from (10⁻⁵)^o to 90^o with a step size of 0.5^o, a_pl is sampled from 0.1 to 4.0 au with a step size of 0.05 au, and finally e_pl is sampled from 10⁻⁵ to 1.0 with a step size of 0.05. Then, within a block of [(i_pl or e_pl), a_pl], each of the ICs mentioned above are set to evolve for up to 10 Myr and 1 Gyr. During this integration period, the close encounters, ejections, and collisions between the planets and the host star are allowed to occur, which marks the stoppage of the integration processes for those ICs. If the integrated orbit survives the total simulation time, then we consider it to be a stable orbit. However, in some cases, when the integrated bodies eject or collide during close encounters, displaying the instability of the system, we note the time of such events and use that time to create a global dynamical lifetimemap that displays dynamically stable or unstable regions.

3.1.1. Analysis of the lifetime map: GJ 832c

To explore the dynamics of the inner planet, its lifetime map (Figure 1) is created for multi ICs in i_pl and a_pl phase space. Each grid point of the map represents one IC. The color coding of each IC gives the survival time of the planet in the respective simulation; however, when a smaller number of data points is available, the interpolation method is used to construct data points in the neighborhood of the known points. The color codes are in the z-axis with the index given in the right-hand color bar. The dark blue color in the map indicates the survival of the planet for 10 Myr, which in this case corresponds to the total simulation time. Lighter colors (from blue to white in our color bar) represent unstable dynamical configurations, indicating that the planet was ejected from the system or collided with the star or the outer planet in less than 10 Myr. The vertical dashed lines in the figure, labeled GJ 832c and GJ 832b, represent the best-fit semimajor axes of the planets. At the best-fit location, GJ 832c remains in a stable orbit for the total simulation period with an orbital inclination as high as 70°. With longer simulation times (>10 Myr), this inclination regime may significantly be reduced, which we explore in the following sections. On the contrary, during the 10 Myr integration period the inner planet with an initial orbital inclination larger than 60°–70°, either collided with the star or was ejected from the system. This is because when the perturber (GJ832b) has an eccentric orbit, the inner planet's orbital inclination and eccentricity may reach extremely large values due to the Lidov–Kozai (Kozai 1962; Lidov 1962) effect (Ford et al. 2000; Naoz et al. 2013; Satyal et al. 2014). This effect, especially when the i_pl crosses a certain limit, eventually leads the system toward instability, which is observed in Figure 1.

3.1.2. Maximum Eccentricity in a_pl, e_pl, and i_pl Phase-Spaces: GJ 832c

The maximum eccentricity (e_max) map is shown in Figure 2 for the same phase space, i_pl − a_pl, as in Figure 1. The color coded z-axis in the map represents the maximum orbital eccentricity attained by the inner planet during the total integration period. This e_max map was obtained by integrating 14,400 ICs for 10 Myr and recording the maximum eccentricity for each state during their time evolution. Therefore, the e_max value for each of the states may or may not be equal to the final eccentricity attained at the end of the simulation. In addition, a major fraction of these maximum values, especially in the unstable regions and around its borders with the stable regions, is expected to rise with significantly increased simulation time (see Section 3.1.3 for 1 Gyr simulations).

The e_max map in Figure 2 shows the maximum eccentricity values that evolved from their nominal values during 10 Myr simulation, and for different ICs of the inner planet's orbital inclination which is plotted along the y-axis. The planet's e_max value increases with increasing i_pl along the best-fit a_pl value (0.163 au, white dashed line), observed in the figure in various color schemes (blue-green-orange-pink regions). The planet's e_max remains less than 0.5 in the regions where the i_pl is less than 45°. For higher i_pl values, e_pl quickly increases to 1. Thus, the map suggests that the likely i_pl value is less than 45°. Note that the blue stable region for i_pl > 80° and 1 < a_pl < 2 indicates e_max < 0.1. This is because the planet attained maximum eccentricity (∼1) in a time shorter than our data sampling period. Hence, the simulation algorithm recorded the initial eccentricity value instead. Furthermore, it is unlikely for the planet to have i_pl greater than 45° and still maintain a stable orbit because the observed trend in the map shows the e_max quickly increases to 0.8 and 0.9. Therefore, with increasing i_pl, the e_max value is only expected to rise for longer simulation time. Hence, the blue-green region is the most dynamically stable zone for the inner planet. The dynamics of the region between the two planets, especially GJ 832c, changes with an additional bodies injected between them, which we explore in Section 3.2.

The inner planet's e_max data shown in e_pl and a_pl phase-space (Figure 3) also reveals vertical structures corresponding to the location of mean-motion resonances with the outer planet and it complements, in terms of stability and instability regions, the e_max map in Figure 2 plotted for different phase-space. The y-axis displays different initial conditions for e_pl, while the color coded z-axis shows the e_max attained by the respective initial conditions. The e_max value demonstrates lesser variation from its nominal value, seen in the horizontal blue region along the red asterisk at 0.163 au in Figure 3. The best-fit e_pl value is along the best-fit a_pl values (vertical dashed lines) for both planets and are denoted by the red asterisks. When the e_pl is set at 0.31, the upper uncertainty value, it does not show any significant deviation during the full integration period, suggesting that the uncertainty in the e_pl cannot be constrained any further from this analysis.

3.1.3. Effects of Uncertainty Limits in the Stability of the System

To check how the uncertainty of the orbital elements in Table 1 affects our conclusions on the stability of the system, we integrated the orbits with the upper uncertainty values of the planets' masses, eccentricities, mean anomalies, and arguments of periapsis. The upper limit of the semimajor axis is chosen for the inner planet and the lower limit for the outer planet. This allows us to see how the stable phase space observed in the above sections changes when the outer giant planet, for its upper limit in mass, is placed closer to the inner planet.

The maximum eccentricity reached by the inner planet after 1 Gyr of orbital evolution is plotted in the e_max map shown in Figure 4. Compared to the previous maps, this map has lower resolution (∼1400 ICs in x and y axes). Nevertheless, the e_max for the semimajor axes between 0.169 au and 1.5 au shows no significant deviation from the assumed initial e_pl. However, beyond 1.5 au, the map shows a rise in the instability islands, where the e_max is found to reach values close to 1. For most of the cases when the inner planet's eccentricity reached more than 0.5, it was ejected from the system or collided with the outer planet.

Comparing the 10 Myr simulation map (Figure 2) with the 1 Gyr simulation map (Figure 4), we see that the planet's orbital inclination between 40° and 50° shows the eccentricity variation from 0.3 to 0.4 and 0.4 to 0.5, respectively. That is, a 0.1 eccentricity variation is observed in the same inclination regime for the increased simulation time. The unstable region, as expected, has extended further inward from the location of the outer planet, which is set at 3.28 au, the lower limit of its semimajor axis. Also, the observed resonance structures in Figure 2 have started to diffuse within the phase-space, which is primarily due to the lower resolution map and shifting of the high-mass outer planet closer toward the inner planet.

For the lower limit of the uncertainty values of the orbital elements and the planetary masses, we expect the stability region to remain unaltered, if not widened a little, when compared to the map shown in Figure 2, where the nominal values of the best-fit orbital parameters are used to integrate the system.

3.2. Dynamics of an Additional Planet

The known planetary configuration in this system shows a super Earth orbiting at a close proximity, 0.163 ± 0.006 au from the host star, while a gas giant orbits distantly at 3.56 ± 0.28 au (for reference, Mercury orbits the Sun at 0.39 au and Jupiter at 5.2 au). Therefore, the existence of other Earth-mass planet(s) (could be bigger or smaller than Earth) between the inner and the outer planets is a plausible scenario. As observed in Figures 1–4, a planet with mass equal to 5.4 M_⊕ and orbiting between 0.1 au and 2.0 au is dynamically stable for a wide range of initial orbital eccentricities and inclinations. However, it is also important to note that even inside 2.0 au unstable orbits also exist as, for example those near mean-motion resonances. To observe how the stability of the system changes with additional bodies, we injected a third planet (middle planet from herein), having 1 M_⊕ and studied its orbital dynamics in similar phase-spaces, which are discussed in Section 3.1.

Orbital parameters for the middle planet are chosen based on the stability zone observed in Figures 1–4. For each initial configuration, i_pl is varied from (10⁻⁵)°–90°, a_pl from 0.1 to 4 au, and e_pl, Ω_pl, and ω_pl are set close to zero, while M_pl is randomly chosen between (10⁻⁵)°–360°. The mass is first set at 1 M_⊕, and later raised up to 15 M_⊕ to observe the orbital variations within the stability zone. The nominal best-fit values of the orbital parameters are used for the inner and outer planets while integrating the middle planet. Later, in the second case, while studying the phase-space of the middle planet, the upper limit of the known planets' orbital parameters, including their masses, are considered and the system is integrated for up to 1 Gyr.

3.2.1. Lifetime and Maximum Eccentricity Maps: the Middle Planet

The lifetime map and the e_max map (Figures 5 and 6) of the injected middle planet with 1 M_⊕, are generated from its survival time in the orbits and the maximum eccentricity attained by the orbits during the total integrating period of 10 Myr, respectively. These maps indicate a wide stability region in the i_pl − a_pl phase-space. The blue region in the lifetime map, and the blue-green region in the e_max map, extends from 0.2 au to 2.2 au along the semimajor axis and 0° to 40° (on the average) along the inclination axis. The map shows sharp changes in the e_max values along the ∼40° i_pl value, seen along the horizontal strip of the cyan-blue marker. This is the Lidov–Kozai (Kozai 1962; Lidov 1962) resonance above which the planet is orbitally unstable. The e_max values in the region of phase-space, which is below 40° is less than 0.15 and all the orbits survive the total simulation time. Outside this stable zone, the middle planet's eccentricity was forced to 0.4 or higher in most of the cases, which caused it to either collide with the inner planet or get ejected from the system. In either case, the planet lost its orbital stability.

**Figure 5.** Global dynamical *lifetime* map of an Earth-mass middle planet in the GJ 832 system injected in the dynamically stable zone (observed in Figure 2) between the inner and the outer planet. The map represents the survival time of the middle planet for 14,400 ICs in i_pl − a_pl phase-space and simulated for 10 Myr. The index in the color bar indicates the survival time, where the lighter colors represent the instability (ejection or collision) and the dark blue color represents the stability (survival) up to the integration period. The vertical dashed lines are the locations of the best-fit semimajor axis of the two known planets. The dark blue region indicates the stable orbital configuration for the Earth-mass planet.
Download figure:
Standard image High-resolution image

**Figure 6.** Maximum eccentricity (e_max) map of the Earth-mass planet injected between the inner and outer planets, in the i_pl and a_pl phase-space and simulated for 10 Myr with 14,400 ICs. The index in the color bar indicates the e_max the orbit evolved into, starting from 10⁻⁵ during the total simulation time. The e_pl varies less from 0.2 to 2.2 au and below the 40° inclination regime. In any other cases, the planet suffered an ejection or collision with the inner or outer planets, especially when e_pl reaches a value greater than 0.5. The dark blue-green color represents the best-fit e_pl parameter from observation (0.18) and other light color represents the e_max value the planet attained for the respective choices of initial conditions in i_pl and a_pl. The vertical dashed lines are the best-fit semimajor axis of the two known planets, and the resonances due to the outer planet are similar to that observed in Figure 2.
Download figure:
Standard image High-resolution image

3.2.2. Long-term Orbital Stability: the Middle Planet

The lifetime and e_max maps discussed in Section 3.2.1 have shown a region where the 1 M_⊕ planet can maintain a stable orbit between the two known orbits, and for up to 10 Myr. For similar phase-space, we then integrated the system for 1 Gyr and raised the mass of the middle planet to 15 M_⊕. The primary reason we chose a 15 M_⊕ planet is because, if the system is stable for this mass, it is more likely to be stable for lower mass planets that could exist between the inner and outer planets. The other reason is discussed in Section 3.3 on how the maximum amplitude of the RV signal is produced by a 15 M_⊕ planet when the sin(i) is considered to be maximum (i ∼ 90°). In addition, we set the masses, eccentricity, and argument of periapsis of the inner and outer planet to their upper uncertainty limit as given in Table 1. Also, the orbit of the inner planet is moved outward to its maximum value (0.169 au) and the outer planet is moved inward to its minimum value (3.28 au) to minimize the space between them and see how that affects the dynamics of the middle planet.

Figure 7 shows the maximum eccentricity of the 15 M_⊕ middle planet in i_pl − a_pl phase-space after 1 Gyr orbital evolution with ∼1400 initial conditions in vertical and horizontal axes. The semimajor axes of the inner and outer planets are denoted by the black dashed lines. The color bar indicates the maximum eccentricity attained by the middle planet during the total integration period. The dark blue-green regions suggest a potentially stable zone, where the e_max is less than 0.2, while the other brighter colors represent potentially unstable zones where the e_max reaches close to 1. The eccentricity is found to deviate least from its initial value of 10⁻⁵ around the 1 au marker and remains less than 0.1 for i_pl up to ∼30°. The vertical structure observed close to 1.35 au corresponds to the location of the near 4:1 mean-motion resonance with the outer planet and in our case all selected orbits around this region are unstable ones.

The stability regions get narrower when compared to the phase-space map generated for the 1 M_⊕ planet observed in Figure 6. The inclination regime has reduced to ∼15° for the a_pl between 0.5 and 0.75 au. Also, the overall i_pl regime for the middle planet is less than ∼30° and extends roughly from 0.5 to 2.0 au. The map suggests that the best possible location for the middle planet would be in the proximity of the 1.0 au marker, where the eccentricity undergoes the least deviation from its initial value.

The phase-space around the GJ 832c at 0.169 au (black dashed line) has undergone the least eccentricity variation compared to the e_max map in Figure 6. This does not mean that the injected 15 M_⊕ planet is stable in that region. The maximum eccentricity has remained very close to ∼0.1 (represented by the blue color in the map) because the high-mass injected planet collides with the low-mass inner planet (5.4 M_⊕), changing the orbital configuration of the latter and making it dynamically unstable. During this collision, the maximum eccentricity attained by the 15 M_⊕ planet deviates least from its nominal value, and this is what we have plotted in the map. Also, in the higher inclination regimes, the inner planet displays more variation in its eccentricity compared to the middle planet, which we discuss in the next section.

3.2.3. Eccentricity and Inclination Time Series of the Planets

To confirm that the dark blue-green regions as seen in Figure 7, which we claim to be an orbitally stable zone for the middle planet, continues to remain in long-term stable orbits and to observe a time-series evolution of e_pl, we picked a few initial i_pl points along the 1 au mark (see Figure 7, red asterisks) and re-integrated for the selected ICs. The other orbital parameters that are used to set up the ICs are similar to the ones we discussed in Section 3.2.2. The masses of the known planets are set at their upper uncertainty limit of the best-fit values given in Table 1, and the injected middle planet's mass is set at 15 M_⊕. The time evolution of the middle planet's eccentricity for three different orbital orientations (0°, 30°, and 39°) are given in Figure 8. The initial e_pl and a_pl were set to 10⁻⁵ and 1 au, respectively, and the planet was allowed to evolve in the gravitational influence of the two known planets and the star.

Time-series e_pl evolution (Figure 8, top plots): For i_pl = 0° (Figure 8, top panel), the amplitude of the eccentricity oscillations remains near the initial values for the known planets and varies between 0 and 0.08 for the middle planet. The observed variation in the eccentricity time series is less significant until the initial inclination (i_o) is set above 30°. For i_o = 30°, the eccentricity time-series of the inner planet starts to display larger amplitude oscillations (0 to 0.1); however, no violent end is observed during the billion-year simulation time. Finally, for i_o = 39°, the inner and the middle planets' eccentricity time-series vary significantly in their amplitudes, displays chaotic evolution, and eventually evolve closer to 1. Then, the middle planet collides with the inner planet at around 1.081 Gyr leading to the system's instability.

Similarly, we looked at the i_pl evolution (Figure 8, bottom plots) of all three planets for three different cases when the middle planet's initial inclination (i_o) is set to 2°, 30°, and 39°, and eccentricity to 10⁻⁵. The amplitude of the i_pl oscillations start to rise significantly for higher i_o. For example, when i_o is set to 30° and 39°, the inner planet's inclination reaches 10° and 60°, respectively, from its near co-planar orbit and displays more chaotic orbits for the latter case. Recent studies (Barnes et al. 2016) of such similar chaotic orbits, with large variations in eccentricity and inclination, have shown stable orbits for up to 10 Gyr.

The Lidov–Kozai resonance occurs at 39 fdg 2 (see Kozai 1962; Lidov 1962; Innanen et al. 1997), beyond which the anti-correlation between the e_pl and i_pl excites the orbits into high eccentricities, significantly reducing the periastron distance and leading to a collisional path (for detailed Kozai resonance analysis, see Satyal et al. 2014). Also, a recent study of the planets in circumbinary orbits of α Centauri AB by Quarles & Lissauer (2016) shows how the Lidov–Kozai resonance beyond 40° limits the stability region. Therefore, we believe that the maximum inclination for the planets in the GJ 832 system, which are orbiting interior to the outer planet, is less than the critical angle of 39 fdg 2.

The above discussion is for the case when initial eccentricity (e_o) of the middle planet is set to near zero. Our separate tests show that the system looses its stability due to a collision between the inner and middle planets shortly after 700 Kyr of integration time when e_o is set above 0.4. For these sets of simulations, the initial orbital inclination was set close to zero in all cases, and the other orbital parameters remained the same as discussed above.

3.2.4. Orbital Resonances in the Presence of the Middle Planet

Figure 6 does not exhibit any orbital resonant structures corresponding to the expected locations of MMRs, especially in the outer regions of the inner planet's best-fit location at 0.163 au. We expect such structures to emerge with significantly longer simulation time and in a higher resolution map. In addition, the prominent vertical structures are observed beyond 1.0 au, which corresponds to the unstable regions near the locations of MMRs with the outer planet. These resonances at the 1.40 au, 1.70 au, 1.93 au, 2.00 au, 2.25 au, and 2.70 au are due to the 4:1, 3:1, 5:2, 7:3, 2:1, and 3:2 MMRs, respectively. Similar resonance structures are observed in the solar system's asteroid belt, where the Kirkwood gaps are in the 3:1, 5:2, 7:3, and 2:1 resonance with Jupiter (Moons & Morbidelli 1995). However, not all resonances are unstable. For example, Jupiter's moons (Io, Europa, and Ganymede) form a resonant system with 1:2:3 orbital resonance and maintain stable orbits.

The vertical structures between 1.5 au and 2.5 au corresponding to the location of mean-motion resonances fade away and shift inwards when the outer planet's location is moved to 3.28 au, as shown in Figure 7. This is primarily due to the change in the outer planet's orbital parameters and the reduced resolution of the phase-space map.

3.2.5. Orbital Stability of Test Particles beyond the Outer Planet

We performed a test particle simulation to check for a possible stable orbit of a planet in a region exterior to the outer planet. For this, we simulated 21,000 test particles (TPs) scattered between 0 and 8 au in the presence of the inner and outer planets. These planets are set at their best-fit semimajor axis location, but their masses are set to the upper uncertainty limit. The initial eccentricity and inclination of the TPs are set to 10⁻⁵ and (10⁻⁵)°–60°, respectively, and the other orbital parameters (ω, Ω, and M) are randomized between (10⁻⁵)°–360°. The TPs do not interact with each other, but they do interact with the known big bodies and evolve due to their gravitational influence.

In Figure 9, we have plotted the maximum eccentricity versus the semimajor axis of the TPs that survived the 10 Myr simulation time. At the beginning of the simulations, the TPs are assigned different initial inclinations (i_o), shown in the figure legend. For each i_o, the orbits of the TPs evolve in time and, based on their orbital configuration, they either collide with the planets/star, get ejected from the system due to close encounters, or continue to evolve in their orbits. The survival rate of the test particles is high between the known planets (denoted by two vertical dashed line) and exterior to the outer planet, in the region between 5 and 8 au. Since there are no known bodies beyond the outer planet to constrain the orbital configuration of the test particles, the outward stability region simply continues. Therefore, we chose the cut-off mark at 8 au.

Most of the surviving TPs between 0.163 and 3 au whose i_o is between 50° and 60° evolve into very high eccentric orbits even though their e_o was set close to zero. For higher starting i_o values, the gravitational effects from the existing planets and the star push the TPs into higher eccentric orbits. And, when the i_o is set at relatively less inclined orbits, between (10⁻⁵)°–40°, the majority of the test particles' maximum eccentricity remains below 0.2. This complements our analysis in Section 3.2.

The TPs beyond the outer planet maintain stable orbital configuration from 5 to 8 au for all i_o and their maximum eccentricity remains less than 0.3 for most of the cases. Few test particles are locked in mutual co-orbitals with the outer planet, seen along the dashed line at 3.56 au. The survival of the test particles in the phase-space between 5 and 8+ au suggests a potential dynamically stable region for additional bodies. However, any planet residing in this outer region may have minimal to no influence on the observed dynamics of the GJ 832b and the injected middle planet, unless that planet has sufficiently large mass. But, the radial velocity observations of GJ 832b excludes such a scenario. In this paper, we will not explore any further into the dynamical regions beyond the outer planet.

3.3. Analysis of the Synthetic RV Signal

The orbits of the inner and the middle planets were simulated for 1500 days, with a high data sampling period (1 per day). The middle planet was first set at 1 au with assigned masses of 1 M_⊕, 5 M_⊕, and 10 M_⊕ in near circular orbits and integrated separately for each mass. Then, using the integrated data, we generated a set of synthetic RV curves based on the RV signature of Keplerian motion given by Equation (1), adapted from (Seager 2010).

The radial velocity equation is given as

$\begin{eqnarray}&&{V}_{r}=\sqrt{G/({m}_{1}+{m}_{2})a(1-{e}^{2})}\cdot {m}_{2}\sin ({i}_{\mathrm{sky}})\\ &&\,\cdot \,[\cos (\omega +f)+e\cos \omega ],\end{eqnarray} \tag{ 1 }$

where G is the universal gravitational constant, m₁ is the stellar mass, m₂ is the planetary mass, a is the planet's semimajor axis, e is the eccentricity, i_sky is the orbital inclination with respect to the sky plane, ω is the argument of periapsis, and f is the true anomaly.

The caveat of the RV technique is that it renders the minimum planetary mass, msin(i_sky). Thus, it requires information about the orbital inclination with respect to the sky plane in order to better constraint its true mass, M_true = M_min/sin(i_sky). M_true is minimum when i_sky = 90°. For the RV signal analysis in this section, we have considered the best-fit values to be the true mass of the inner and outer planets (that is, i = 90° with respect to the skyplane). Also, for the injected middle planets with varying mass, their msin(i_sky) are assigned 1 M_⊕, 5 M_⊕, and 10 M_⊕, assuming that i_sky = 90°. Since the planetary mass is the function of its orbital inclination (skyplane), the amplitudes of the RV curves we generated may vary if the i_sky is less than 90°. The data presented here is a special case when the orbital inclination with respect to the sky plane is maximum for all the planets. Other orbital parameters are the same as the nominal best-fit values given in Table 1.

The synthesized RV curves for the inner planet (GJ 832c) and the middle planet for three different masses are plotted in Figure 10, top panel (a). The maximum amplitude of the RV signal for the inner planet is ∼2.0 m s⁻¹, similar to the observational value reported by Wittenmyer et al. (2014). Additional planets with bigger mass, and beyond the orbits of the inner planet, can produce a higher amplitude RV signal, however, the observation has constrained the RV signal for any new planets to be less than 2.0 m s⁻¹. For this reason, we limited the planetary true mass at ≤10 M_⊕ (for maximum i_sky) when placed at 1 au because the higher mass planet would produce an RV signal greater than 2.0 m s⁻¹. However, it is possible that a planet can have larger true mass and produce the same RV signal for the i_sky, which is significantly less than 90°.

The injected planets at 1 au generate the RV curves with varying amplitude as expected. For the masses 1 M_⊕, 5 M_⊕, and 10 M_⊕, the RV signal is 0.14 m s⁻¹ (black), 0.70 m s⁻¹ (blue), and 1.04 m s⁻¹ (green), respectively (Figure 10, (a)). The RV amplitude for 1 M_⊕ is only ∼0.14 m s⁻¹, which is much smaller than the current high accuracy RV precision of about 0.97 m s⁻¹ of the HARPS instruments (Mayor et al. 2003). A planet interior to 1 au could be less than 1 M_⊕ but no greater than 10 M_⊕. The planet bigger than 10 M_⊕ (for example, 15 M_⊕) had RVs greater than 2.0 m s⁻¹, thus we disregarded the results. Any middle planet at 1 au would have the orbital period of about ∼550 days. The highlighted RV signal for 1 M_⊕ is shown in panel (b), with smaller y-axis variation.

The other two planets with masses 15 M_⊕ and 20 M_⊕ injected separately at 2 au to obtain their synthetic RV reveal that ∼15 M_⊕ is the upper mass limit for the middle planet (Figure 10, bottom panel (c)) because its RV signal is measured less than 2.0 m s⁻¹. We chose to set the injected planet at the 2 au marker because it is the farthest stability region observed in Figure 7. The RV signals generated for 15 M_⊕ (black) and 20 M_⊕ (blue) planets are 1.50 m s⁻¹ and 2.10 m s⁻¹, respectively. Therefore, we can exclude a planet with mass bigger than 15 M_⊕. The orbital period for a planet with mass ∼15 M_⊕ is close to 1400 days. If the location of an injected planet is farther out toward the outer planet, the probable new planet has relatively higher mass than when it is closer toward the inner planet. At the same time, it is possible that the planet can have a mass bigger than 15 M_⊕ and produce the same RV signal for the i_sky less than 90°.

3.3.1. Observation Prospect of an Additional Planet

Now, based on the semi-amplitude values, K = (V_r,max − V_r,min)/2, and a single measurement precision of σ we can estimate the minimum number of observations required to detect an exoplanet (adapted from Plavchan et al. 2015), and is given by

$\begin{eqnarray}&&{N}_{\mathrm{obs}}=2{(\mathrm{SNR}\cdot \sigma /K)}^{2},\end{eqnarray} \tag{ 2 }$

where the SNR is the detection confidence. The calculated N_obs based on Equation (2) for 5σ detection are given in Table 2. The N_obs for 1 M_⊕, 5 M_⊕, and 10 M_⊕ planets at 1 au is 2500, 103, and 47, respectively; and for 15 M_⊕ at 2 au, it is 23. The N_obs increases significantly for higher σ detection.

Table 2. Radial Velocity Semi-amplitude, K = (V_r,max − V_r,min)/2, for Varying Mass and Location of an Injected Planet, and Estimated Number of Observations for Different Masses and up to 5σ Detection

Injected Planet Mass	Dist. from Star (au)	RV Semi-amplitude (m s⁻¹)	No. of Observations
1 M_⊕	1	0.14	2500
5 M_⊕	1	0.70	103
10 M_⊕	1	1.40	47
15 M_⊕	2	1.50	23

Download table as: ASCII Typeset image

4. Summary

Our studies of the lifetime maps, e_max maps, and the time evolution of the orbital elements for GJ 832c establishes the stable orbital configuration for its best-fit orbital solutions. The planet's eccentricity deviation remained within the best-fit uncertainty values during the total simulation time. Based on the e_max maps for the phase-spaces in i_pl, e_pl, and a_pl, the relative inclination of the planet is less than 39°. Also, the planet remains in stable orbits while maintaining low enough eccentricity deviations during the total integration period of 1 Gyr. The outer planet's orbital elements displayed the least deviation from their initial values during the one billion years of evolution time; however, due to its gravitational influence, the region starting from 2 au to 3.56 au remains dynamically unstable. A region similar to the solar system's Asteroid belt is likely to exist in the vicinity of 2 au, where the 3:1, 5:2, 7:3, and 2:1 resonances with the outer planet are observed.

GJ 832c maintained stable orbit for the i_pl ≤ 39° and the orbit did not vary significantly even when the middle planet with mass 15 M_⊕ was injected into the system and integrated for 1 Gyr. The middle planet also remained in a stable orbital configuration for the orbital inclination as high as ∼30° and for the semimajor axis ranging from 0.75 to 1.25 au. The e_max attained by the middle planet around the 1 au marker remained close to its initial value even after 1 Gyr orbital evolution. In general, based on the small variations in the initial values of the orbital elements, we could extrapolate our stability analysis and claim that the system is likely to maintain a stable orbital configuration for longer than the 1 Gyr timescale.

The injected middle planet could be smaller or bigger than one Earth mass. However, its upper mass limit is constrained by the RV signal of the known inner planet, 2.0 m s⁻¹ for maximum i_sky. Using this RV signal as a constraint, the synthetic RV of the middle planet is generated from our simulation data for its varying mass and i_pl for the maximum i_sky. Our results suggest that if the middle planet is located in the vicinity of the 1 au marker, it has an upper mass limit of 10 M_⊕ and generates an RV signal of 1.4 m s⁻¹. The 1 M_⊕ planet at the same location has an RV signal of 0.14 m s⁻¹ only, which is much smaller than the sensitivity of available technology. Nonetheless, the detection using the RV method is possible but with a significantly large number of observations: ∼2500 for 1 M_⊕ planet. The number of observations can be lowered to 47 if the planet has a mass of 10 M_⊕. The N_obs depends on the preferred σ detection value as well. Similarly, when the middle planet was fixed to 2 au, the upper mass limit increased to 15 M_⊕ with a synthetic RV signal of 1.5 m s⁻¹. Hence, we expect a planet with a mass less than 15 M_⊕ orbiting between the inner and the outer planets. Our RV signal calculation considers only 2° variation in i_sky. If the i_sky varies larger, the mass of the planet will vary according to msin(i_sky). The orbital periods of the planet at 1 au and 2 au are ∼500 and ∼1400 days, respectively. Our aim here is just to provide a general idea of the detection probability to the RV observation scientists. The Earth, for example, exerts 9 cm s⁻¹ wobble in the Sun. Thus, to detect any Earth-mass planets, the sensitivity of the RV measurements should be down to a few centimeters and require extreme precision radial velocities (Fischer et al. 2016).

The lower stability limit for an Earth-mass planet starts at 0.25 au and the upper limit of the star's classical habitable zone (HZ) ends at 0.28 au (from Kopparapu 2013). Hence, there is a slim window of about 0.03 au where an Earth-mass planet (or ≤15 M_⊕) could be stable as well as remain in the upper limit of the stellar HZ. However, a planet residing in an HZ does not necessarily imply that it can support life as we know it. For example, the inner planet, GJ 832c, itself orbits around the lower limit of the stellar HZ but it is not expected to be habitable; see Wittenmyer et al. (2014) for a detailed analysis of GJ 832c as a habitable-zone super Earth.

The upper limit of the planetary mass of GJ 832c is in the range of super-Earths (exoplanets bigger than Earth but smaller than Uranus (15 M_⊕) and Neptune (17 M_⊕)) and has close proximity to the star. These super-Earths can have two formation scenarios: they can form far out and migrate inward to their current location or in situ formation, which is possible if the planetary disk contains a low turbulence region (Martin & Livio 2016). Also, Izidoro et al. (2014) recently found that fast-migrating super-Earths have a modest effect on proto-planetary embryos and planetesimals leaving enough materials to form rocky, Earth-like planets. Hence, if the planet GJ 832c was a fast-migrating super Earth, that would have left enough space and matter in between the inner and outer planets to have additional bodies.

Long-term orbital stability, orbital dynamics of various phase-spaces, and the synthetic RV signal analysis suggest the possible existence of a planet ≤15 M_⊕ between the inner and outer planets in the GJ 832 system. The anticipated RV signal is much lower than the sensitivity of the RV instruments; however, a significantly large number of RV observations and the transit method, provided that the planets are along the line of sight, are viable options to get the observational verifications. Despite being at close proximity (16 ly), the star has aged enough (∼5.6 Gyr) that the residing planets may not radiate enough in order to be detected from the imaging technique. A future space telescope, such as Transiting Exoplanet Survey Satellite (2017), TESS (https://tess.gsfc.nasa.gov/), whose mission is to survey G-, K-, and M-type stars including the 1000 closest red dwarfs, is one of the best options to explore more about this system. In addition, this system is a good candidate for the Jason Webb Space Telescope (2018), JWST (http://www.jwst.nasa.gov/) to perform spectroscopic analysis of the planetary atmosphere.

We would like to thank the referee for a very comprehensive report on our paper, which allowed us to greatly improve the original version. We would also like to thank the Office of Graduate Studies at University of Texas at Arlington and their I-Engage Mentoring Program, which initiated this research project to enhance the undergraduate research program at UTA for J.G. Also, Z.E.M. acknowledges the support of this research by the Alexander von Humboldt Foundation. Special thanks to B. Quarles, M. Cuntz, and J. Noyola for their discussions, comments, and suggestions; and to Mark Sosebee at the high energy computing facility where most of our numerical simulations were carried out.

Dynamics of a Probable Earth-mass Planet in the GJ 832 System

Article metrics

Permissions

Author e-mails

Author affiliations

Dates

Abstract

1. Introduction

2. Numerical Simulations