Dynamical Packing in the Habitable Zone: The Case of Beta CVn

The boundaries of the HZ have been cataloged for a variety of known exoplanetary systems (Kane & Gelino 2012), including the candidate exoplanet systems discovered by the Kepler mission (Kane et al. 2016; Hill et al. 2018). The stellar distance plays a key role in determining the stellar properties that influence the HZ boundaries, highlighting the importance of accurate distance estimates (Johns et al. 2018; Kane 2018). Here we calculate the width of the HZ as a function of the stellar mass, since mass is one of the more critical intrinsic stellar properties that determines the overall stellar evolution pathway.

To investigate the width of the HZ region as a function of stellar mass for main-sequence stars, we utilize the MESA Isochrones & Stellar Tracks (MIST), described in detail by Choi et al. (2016) and Dotter (2016). We adopted a solar metallicity and used an age of 3 Gyr in order to encompass isochrones for stars more massive than the Sun, including a mass range of 0.137–1.199 M_⊙. For each stellar mass, we used the isochrone luminosity and effective temperature to calculate the HZ boundaries as formulated by Kopparapu et al. (2013, 2014). We adopt both the conservative and optimistic HZ boundaries, described in detail by Kane et al. (2016).

The extracted isochrone stellar parameters of luminosity (solid line) and effective temperature (dashed line) along with our calculations of the HZ width (dotted–dashed line) are represented in Figure 1. The values for all three have been normalized to solar values, which is why they intersect at a normalized quantity of unity for a solar mass star. We further calculated the Hill radius for an Earth-mass planet located in the middle of the HZ for each stellar mass. These Hill radii were also normalized, this time to the Hill radius of Earth, and are represented in Figure 1 as a dotted line. The rise in Hill radius with stellar mass is slightly less than the rise in HZ width. This implies that dynamical constraints on the number of planets in the HZ may increase with increasing stellar mass.

Zoom In Zoom Out Reset image size

To test the dynamical limits of terrestrial planets within the HZ, we conducted an extensive suite of N-body integrations using the Mercury Integrator Package, described in detail by Chambers (1999). The simulations adopted the hybrid symplectic/Bulirsch–Stoer integrator with a Jacobi coordinate system, generally providing more accurate results for multi-planet systems (Wisdom & Holman 1991; Wisdom 2006). For each of the stellar masses extracted from the isochrones described in Section 2.1, we performed the following steps.

• 1.  
  The boundaries of the optimistic HZ were calculated based on the stellar parameters.
• 2.  
  Earth-mass planets were placed in circular orbits within the HZ, with one at the inner edge, one at the outer edge, and the others evenly spaced in between (Weiss et al. 2018). This was carried out for a range of 3–7 planets for each stellar mass.
• 3.  
  The orbital period at the inner and outer edges of the HZ were calculated. The orbital period at the inner edge was used to determine the time step of the N-body integrations. These were set to 1/20 of the inner edge orbital period to ensure perturbative reliability of the simulation, consistent with the recommendations of Duncan et al. (1998). The total simulation time was set to 10⁸ times the orbital period at the outer edge of the HZ.
• 4.  
  The dynamical simulation was carried out ∼5 times with randomized starting locations (mean anomalies) for each of the planets to robustly capture the dynamical stability of the orbital architecture.
• 5.  
  The assessment of stability was determined by requiring that all planets survive the duration of the simulation. Non-survival means that one or more of the planets has been captured by the gravitational well of the host star or ejected from the system.

Following the above steps resulted in several thousand simulations to fully explore the dynamical stability limits of packing evenly spaced terrestrial planets within the HZ for a range of stellar masses.

The results of the simulation suite described above are shown in Figure 2, where the colored shaded regions indicate the approximate spectral types that correspond to the stellar masses. For the three and four planet cases, we found that the dynamical integrity of the systems are retained for the entire range of stellar masses tested. We show the five planet case (solid line) as a demonstration of where instability begins to arise at the extreme low-mass end of the stellar mass range. The six planet case (dashed line) exhibits significant instability within the HZ for stellar masses less than ∼0.7 M_⊙. The trend of relative instability for lower mass stars is consistent with the relationship between the HZ width and the size of the Hill sphere discussed in Section 2.1. Significant collapse of the dynamical integrity for all stellar masses starts to occur for the seven planet case (dotted line), with no stable scenarios for masses less than ∼0.25 M_⊙. Furthermore, the stellar mass range of 0.8–1.1 M_⊙ (G stars) is unlikely to harbor seven evenly spaced planets within the HZ. Calculation of the corresponding orbital periods of the planets for this stellar mass range reveals a cause of mean motion resonances (MMRs). Specifically, the third planet in the sequence of seven planets (from the inner HZ edge to the outer HZ edge) lies close to the 5:4 MMR with planet four, the 3:2 MMR with planet five, and 2:1 MMR with planet seven. These MMRs result in strong perturbative forces on planet three that are able to successfully compromise the system stability and rapidly lead to close encounters and subsequent ejections. Thus, in some cases, MMRs can lead to a potential additional advantage of K dwarfs over G dwarfs for biosignature and direct imaging studies (Arney 2019).

Zoom In Zoom Out Reset image size

For example, the case of our solar system has an optimistic HZ that extends from 0.74 to 1.73 au. Our simulations show that six terrestrial planets are able to remain stable in this region but the seven HZ planet architecture falls within the regime of MMR perturbations and catastrophic instability. This is consistent with the results of Smith & Lissauer (2009) and Obertas et al. (2017) whose simulations were for the specific solar mass case. These simulations neglect the effect of giant planets beyond the HZ that can also induce perturbative effects that can disrupt even the six planet scenarios.

The number of known exoplanetary systems is now sufficient to develop statistical models that can predict orbital distributions based on collisional formation scenarios (Tremaine 2015). However, these distributions are greatly influenced by giant planets, particularly during their early periods of migration (Morbidelli & Raymond 2016). The presence of the giant planets in the solar system has undoubtedly played a significant role in the formation, evolution, and final architecture of the inner solar system (Tsiganis et al. 2005; Batygin & Laughlin 2015; Horner et al. 2020). Such important dynamical effects of giant planets have also shaped the architecture of other planetary systems (Matsumura et al. 2013; Morbidelli & Raymond 2016). Indeed, both the solar system and other planetary systems have been found to be close to instability boundaries that would compromise the dynamical integrity of those systems (Laskar 1996; Lecar et al. 2001; Murray & Holman 2001; Barnes & Quinn 2004). For example, a planet within the asteroid belt of the solar system would significantly alter the orbital dynamics of the inner solar system by exchanging angular momentum with the outer solar system (Lissauer et al. 2001).

With this in mind, the HZ packing scenarios described in Section 2.2 have been considered in isolation but are clearly dependent on the presence of other planets outside of the HZ, especially giant planets. Jupiter has a well-known perturbing influence at resonance locations, especially the induced chaos near the Kirkwood gaps (Wisdom 1983). The case of seven planets within the HZ, described in Section 2.2, suffers from instability due to MMR effects internal to that region. External sources of MMR, particularly 5:2, 3:1, would result in similar sources of instability which would have a cascading transfer of angular momentum, similar to the asteroid belt planet scenario described by Lissauer et al. (2001). Therefore, the most likely cases of maximum terrestrial planets within the HZ are those solar-type stars for which giant planets may be excluded.

The proximity, spectral type, and apparent brightness of Beta CVn have resulted in numerous observational studies of the star. The star has been found to be photometrically and spectroscopically stable, resulting in its frequent use as an RV standard (Konacki 2005). A study of Sun-like stars by Radick et al. (2018) presented 18 yr of Beta CVn photometry, acquired with the Automated Photoelectric Telescopes (Henry 1999). We extracted the Beta CVn photometry from their data collection and conducted a Fourier analysis of the complete time series. The highest peak in the power spectrum from this analysis occurs at 182.5 days, and is an observational artifact from the observing cadence. The second highest peak occurs at 110 days, which could be related to intrinsic stellar variability, although is too long for an expected rotation period.

Given the quiet nature of the star, the question arises as to the stellar equatorial plane relative to the plane of the sky since this will influence the detectability of planetary signatures, assuming those orbits are coplanar with the stellar equator. In order to optimize blind transit searches, Herrero et al. (2012) provided a catalog of stars whose inclination of the equatorial plane is likely > 80° relative to the plane of the sky. They adopt a rotational velocity of $v\sin i=2.9$ km s⁻¹ and activity index of $\mathrm{log}({R}_{\mathrm{HK}}^{{\prime} })=-4.885$ for Beta CVn, and their methodology predicts a high probability (50% larger than the average of their sample) of the star having an inclination $\gt 80^\circ$. The predicted probability is 10% higher than the mean probability for the host stars in their sample that are known to host transiting planets. This in turn increases the probability that any planets present are in close to edge-on orbits, increasing the expected RV signal.

The proximity of Beta CVn has also enabled several direct measurements of the stellar radius using interferometric techniques (van Belle & von Braun 2009; Boyajian et al. 2012). We adopt the stellar parameters provided by Boyajian et al. (2012), which include a stellar radius of R_⋆ = 1.123 ± 0.028 R_⊙, a luminosity of L_⋆ = 1.151 ± 0.018 L_⊙, and an effective temperature of T_eff = 5653 ± 72 K. Based on these stellar parameters, we calculate the conservative HZ as 1.03–1.82 au and the optimistic HZ as 0.81–1.92 au. The significance of the HZ distances for the star are described in more detail in Section 3.3.

Even though giant planets are more common around solar analogs than later-type stars (Zechmeister et al. 2013; Wittenmyer et al. 2016), they are still relatively rare, with an average occurrence rate of 6.7% beyond orbital periods of 300 days (Wittenmyer et al. 2020). Consequently, the solar system may not be representative of orbital architectures, and it is useful to directly compare with another nearby solar-type star, such as Beta CVn. To determine the limits on such giant companions, we present here 235 precision RV measurements of Beta CVn acquired over a period of 20 yr. Observations were carried out using the HIRES echelle spectrograph on the Keck I telescope (Vogt et al. 1994) and the Levy spectrometer on the APF (Radovan et al. 2014; Vogt et al. 2014). A subset of the RV data is shown in Table 1, where all the data have been calibrated to the same zero-point. Table 1 also includes a column for the instrument that was used to acquire the data. Note that the HIRES subscripts refer to data that were acquired prior to the 2004 upgrade to the instrument (k) and those that were acquired after (j). The data from all three instruments are shown plotted in Figure 3.

Zoom In Zoom Out Reset image size

We performed a Fourier analysis of the RV data to investigate possible periodic signals. The periodogram resulting from this analysis is shown in Figure 4. The horizontal dotted line indicates a false-alarm probability (FAP) threshold of 0.001 (0.1%). The vertical dashed lines show the location of common aliases caused by Earth's rotation (diurnal) and orbital (annual) motions, as well as that caused by the lunar orbit. The highest peak occurs at a period of 169.2 days, which is still well below the FAP threshold. We conclude that there are no significant periodic signals present in the data.

Zoom In Zoom Out Reset image size

The RV data described here were used to perform an injection-recovery test to quantify the completeness of the observations for excluding planetary signatures. This method injects planetary signatures of various masses (${M}_{p}\sin i$) and semimajor axes (a) into the data using the same observation epochs and noise properties of the data. Circular orbits are assumed and the Keplerian orbital fits are performed using the RadVel package (Fulton et al. 2018). The methodology for this process, including the criteria for planetary signature recovery, is described in detail by Howard & Fulton (2016). To perform the injection-recovery test in terms of ${M}_{p}\sin i$, a stellar mass must be assumed. There are a variety of estimates for the mass of Beta CVn in the literature, including 1.05 ± 0.14 M_⊙ by Valenti & Fischer (2005) and 0.852 ± 0.023 M_⊙ by Boyajian et al. (2012). We adopt the latter of these to be consistent with the stellar parameters used to calculate the HZ in Section 3.1.

The results of the injection-recovery test are shown in Figure 5. The blue dots represent injected planets that were recovered and the red dots represent those that were not recovered. The shaded contours indicate the probability of a successful detection, with the color scale shown on the right vertical axis. The results thus demonstrate that the RV data are sufficient to rule out the presence of 100 ${M}_{\oplus }$ planets out to ∼5 au and 300 M_⊕ (Jupiter mass) planets out to ∼10 au. These data are thus sufficient to rule out the presence of most giant planets around Beta CVn inside an orbital radius of 10 au. Note however that the data are not sufficient to rule out terrestrial planets in the system.

Zoom In Zoom Out Reset image size

As described in Sections 2.2 and 2.3, the relative lack of giant planets in the Beta CVn system may allow for the presence of numerous terrestrial planets in a HZ dynamical packing scenario. Based on Figure 2 and the stellar mass adopted in Section 3.2, the HZ of the star would be unlikely to contain seven terrestrial planets even without a giant planet due to internal MMR perturbation effects. Beta CVn may therefore be a suitable system for which to test HZ packing scenarios.

As calculated in Section 3.1, the optimistic HZ of the system extends from 0.81 to 1.92 au. Continuing to use the stellar parameters of Boyajian et al. (2012), we calculate several predicted exoplanet signatures at the inner and outer edge of the optimistic HZ, shown in Table 2. These include the orbital period P, the RV semi-amplitude K, the transit depth d, and the transit probability p (Kane & von Braun 2008). As expected, the K values are well below that required of the RV data presented in this work, but may be feasible with developing facilities (Fischer et al. 2016). Additionally, the work of Herrero et al. (2012) showed that Beta CVn has a relatively high probability that the inclination of the equatorial plane is > 80° relative to the plane of the sky (see Section 3.1), and so the transit probabilities are likely underestimated.

However, a further option to explore the possible planetary architecture of Beta CVn, particularly given the proximity of the star, is direct imaging (Kane et al. 2018; Dulz et al. 2020). At a stellar distance of 8.44 pcs, the 0.81–1.92 au optimistic HZ translates into angular separations of 96–227 mas. Consistent with the final reports of the Habitable Exoplanet Observatory (HabEx) mission (Gaudi et al. 2020) and the Large UV/Optical/Infrared Surveyor (LUVOIR) mission (The LUVOIR Team 2019), we assume (a) an inner working angle (IWA) of ∼58 mas, and (b) a contrast limit of ∼4 × 10⁻¹¹, which corresponds to the Earth–Sun contrast ratio at the outer edge of the Sun's conservative HZ (1.7 au).

Given the above assumptions, Figure 6 shows the stars having the largest fraction of their HZs falling within the detection limits of the instrument (scaled by the Hill radius), and the maximum V-band brightness of Earth analogs within the detectable HZ space, assuming a Bond albedo of 0.3. The color of the stars shown in the figure represents the effective temperature (spectral type) of the host stars, with M red, K orange, G yellow, and F cyan; similar to Figure 2. The most promising targets are shown in the right of the plot where both the predicted HZ Earth analog brightness and the visibility of the HZ are largest. For example, Beta CVn has a fractional HZ visibility of ∼75% and an expected Earth analog brightness of V ∼ 29. Note that the stellar sample shown in Figure 6 is dominated by K and G dwarfs for which the occupancy of the HZ by such Earth analogs is maximized, except in cases where internal MMRs truncate the number of allowed orbits (see Figure 2). Thus, direct imaging with facilities that meet the criteria laid down in the HabEx and LUVOIR reports may naturally favor targets that are predisposed toward larger than expected HZ planet inventories.

Zoom In Zoom Out Reset image size