14 Her: a likely case of planet-planet scattering

We have collected 283 archival RV data points from Keck/HIRES (Rosenthal et al., 2021; Butler et al., 2003; Wittenmyer et al., 2007; Hirsch et al., 2021) spanning over 20 years (1997-04-07 to 2020-02-26; Rosenthal et al. 2021). The High Resolution Echelle Spectrometer (HIRES; Vogt et al. 1994), mounted on the Keck I 10-m Telescope, operates between $0.3-1\,\mu$m providing spectral resolutions between R$\sim 25000-85000$ and RV precision down to 1 m/s. The RV measurements for 14 Her fluctuate between $-$73 m/s and 191 m/s, with a semi-amplitude of 100 m/s and a visible long-term slope caused by the c planet. The mean precision for our data is 1.08 m/s. The RV curve of 14 Her is shown in Figure 2.

We also found RV data in the literature from the ELODIE spectrograph at the Observatoire de Haute-Provence and the Automated Planet Finder (APF) at Lick Observatory. The ELODIE data spans epochs from 1997 to 2002, whereas the APF data also covers 5 years but between 2014 and 2019. The epoch span of the ELODIE data was sufficient to discover the b planet and fully determine its orbit (Naef et al., 2004). However, we require a longer baseline to determine the orbit of the c planet, and HIRES is the longest-running instrument with continuous RV monitoring of planet-hosting stars (Butler et al., 2017). Therefore, we decided to only use HIRES RV data to avoid any offsets or potential systematics across instruments. However, HIRES underwent an upgrade in 2004 which effectively changed its RV zero-point ( Tal-Or et al. 2019; see Section 5.3).We conservatively treat the pre- and post-upgrade data as coming from two different instruments which introduces two more degrees of freedom to the orbit fit.

We combine the radial acceleration information with tangential acceleration from absolute astrometry. The long baseline between Hipparcos and Gaia epochs and increasingly better astrometric precision from Gaia EDR3 significantly improve the uncertainty of the proper motion accelerations. We use the absolute astrometry of the Hipparcos-Gaia Catalog of Accelerations (HGCA; Brandt 2021) which has now been updated with Gaia EDR3 astrometry. The HGCA has cross-calibrated the astrometric solutions for all Hipparcos stars that are present in Gaia to a common reference frame (Brandt, 2018).

The proper motion of 14 Her across the Hipparcos and Gaia EDR3 epochs is inconsistent with a constant value by $31\sigma$. 14 Her has a mean acceleration in its proper motion direction of $(a_{\alpha},a_{\delta})=(0.812,-8.918)$ m/s/yr. Its total acceleration ($a_{tot}=8.955$ m/s/yr) is lower than that for 97% of stars in the HGCA, therefore it is a highly significant, low accelerating star.

Our derived stellar and planetary parameters are shown in Table 2. Since our RV data covers nearly five full orbits of the b planet, we are able to set strong constraints on its orbital parameters (Figure 3). The b planet has a mass of M${}_{b}=\,$${9.1}_{-1.1}^{+1.0}$$M_{\mathrm{Jup}}$ and its orbit has a semimajor axis of $a=$${2.845}_{-0.039}^{+0.038}$ AU with a moderate eccentricity of $e=$${0.3686}_{-0.0031}^{+0.0032}$and a ${32.7}_{-3.2}^{+5.3}$ degree inclination with respect to the plane of the sky. Its orbital period translates to $P=$${4.8277}_{-0.0023}^{+0.0022}$ years.

Despite only seeing a long-term trend in the RV data covering $\sim$15% of the period, we can set some constraints on the orbital parameters of the c planet. The mass of the c planet is lower than that of b (M${}_{c}=\,$${6.9}_{-1.0}^{+1.7}$ $M_{\mathrm{Jup}}$) and it is located much farther away from the star, at $a=$${27.4}_{-7.9}^{+16}$ AU. 14 Her c has a highly eccentric orbit ($e=$${0.64}_{-0.13}^{+0.12}$). With respect to the orbit of 14 Her b, the c planet has a broad distribution of inclinations ($i=$${101}_{-33}^{+31}$ degree) but hints at being misaligned (see Figure 4). Its orbital period amounts to $P=$${144}_{-58}^{+139}$ years with large uncertainties. Our $M\sin i$, semimajor axis, and eccentricity values for the b planet are in line with previous results (e.g., Hirsch et al., 2021; Rosenthal et al., 2021). However, introducing the constraint from the absolute astrometry enhances these parameters for the c planet with respect to earlier determinations using only RV data. We calculate a minimum mass of $M_{c}\,\sin i_{c}$ = $6.79^{+1.85}_{-1.03}$ $M_{\mathrm{Jup}}$ which is higher but consistent within $1\sigma$ with published values ($M_{c}\,\sin i_{c}$ = $5.8^{+1.4}_{-1.0}$ $M_{\mathrm{Jup}}$, Rosenthal et al. 2021; $M_{c}\,\sin i_{c}$ = $6.1^{+1.3}_{-0.9}$, Hirsch et al. 2021).

A surprising aspect of our MCMC results is the evidence for discrepant inclinations of the orbital planes for the two 14 Her planets. Orbit misalignment of outer giant planets can have an effect on the final configuration of super-Earths within a system, limiting their maximum mass (Childs et al., 2019), and increasing their eccentricities even to the point of migration (Lai & Pu, 2017), which in turn can have serious consequences on their potential for habitability. We evaluated the coplanarity of the orbits by calculating the 3D angle between the angular momentum vectors of the planetary orbits. The unit vectors are given by:

where $\Omega$ is the longitude of the ascending node and $i$ is the inclination of the orbit. The angle between the two angular momentum vectors is then

Calculating this expression on each chain of our MCMC sampling returns an angle of $\Theta_{bc}=$$96.3_{-36.8}^{+29.1}$ degrees. orvara places uniform priors on the orientation of each orbital plane, and this translates into a $\sin i$ prior on the relative inclination between the two planes. Comparing the distribution of mutual inclinations to the $\sin i$ prior used in orvara to give equal probability to all orbital orientations (Figure 5), we see that mutual inclinations of roughly $0-60^{\circ}$ and $130-180^{\circ}$ are strongly disfavored, indicating that the orbits are misaligned. 14 Her joins the ranks of only 3 other systems with giant planets in misaligned orbits: $\pi$ Mensae (De Rosa et al., 2020), $\nu$ Andromedae (McArthur et al., 2010), and Kepler-108 (Mills & Fabrycky, 2017). In comparison with $\pi$ Men c, whose inclination is not determined and only a plausible range is estimated, we have obtained both inclinations for the 14 Her planets and arrive at a misalignment constraint that is at least as confident as that for $\pi$ Men.

Another key measurement to constrain the orientation of the system is the spin-axis alignment of the star. Discrepant $V\sin i$ measurements for 14 Her are reported in the literature (Morris et al., 2019; Sissa et al., 2016; Valenti & Fischer, 2005; Luck, 2017), giving rise to all range of possible inclinations. However, the brightness of 14 Her makes it potentially amenable for 20 s-cadence TESS observations for asteroseismology to find the projected obliquity of the star (e.g. $\pi$ Men, Huber et al. 2021; $\alpha$ Men, Chontos et al. 2020). A third constraint on the relative angle between the star and the planets would be needed for a full orientation characterization and more clear description of the system’s past history.

The stability of an RV instrument is crucial to detect small variations that could signal the presence of low-mass and long-period companions. The HIRES spectrograph underwent a major upgrade in 2004 which enabled an RV precision of 1-3 m/s, a factor of 3 improvement (Butler et al., 2006). This upgrade changed the RV zero-point by $-1.5\pm 0.1\,$m/s, introduced a long-term drift of $<1\,$m/s, and a small nightly drift (Tal-Or et al., 2019). Rosenthal et al. (2021) modeled the pre- and post-upgrade RV data separately for their planet search without applying the Tal-Or et al. (2019) offset manually. For our results, we conservatively used the Rosenthal et al. (2021) RV data in the same way, without any modifications. We labeled the pre- and post-upgrade data points as if they were coming from two different instruments to let orvara find suitable zero-points and jitter values independently (see Table 2). Comparing the zero-points before and after the upgrade, we find a difference of $\Delta_{ZP}=$$-5.3_{-2.0}^{+2.0}$ m/s, which is larger than the one presented by (Tal-Or et al., 2019).

We tested whether manually adding the Tal-Or et al. (2019) offset to the post-upgrade data would make a difference in our results. We ran four 600,000-step chains and combined them onto a long pseudo-chain with the same settings as for our science results but with a modified RV file with the Tal-Or et al. (2019) offset and two instrument IDs to reflect the pre- and post-upgrade data. The results for both the b and c planet are essentially the same as our science results, including similar jitter values and zero-points. We also tried using the Rosenthal et al. (2021) data as-is with a single instrument ID. In this case, we found that the posterior parameters for the b planet were the same as our science results within uncertainties, given that the orbital period of the b planet is fully covered several times within the post-upgrade epoch range alone. The median of the posterior parameters of the c planet differ slightly from our science results, especially the semi-major axis ($a_{c}={31.1}_{-8.8}^{+16}$ AU), the eccentricity ($e_{c}={0.69}_{-0.11}^{+0.10}$), and the mass ($M_{c}={7.11}_{-0.82}^{+1.5}$ $M_{\mathrm{Jup}}$), although they stayed within the uncertainties of our science results. These differences imply that a small RV offset could have a significant effect in fitting long period planets. We also found more chains caught up in low probability areas compared to our science results. Therefore, it seems that results are robust as long as the pre- and post-upgrade data are treated separately.

In today’s configuration, 14 Her is a stable system. A planetary system is stable against the mutual gravity of its components if their orbital separations exceed $2\sqrt{3}$ times their mutual Hill radii (Gladman, 1993), and the large separation between the b and c planets exceeds the stability criterion by a factor of $\sim 3$. However, we postulate that today’s orbits are likely a far departure from their initial conditions.

Whether resulting from core accretion or gravitational instability, planets are expected to form in disks, causing multi-planet systems to be coplanar as a consequence (Hubickyj, 2010; Mayer, 2010). Multi-planet systems tend to have lower eccentricities than single-planet systems (Wright et al., 2009), possibly because low eccentricites are energetically favorable for long-term stability. The strongly disfavored coplanarity and high eccentricities in the case of the 14 Her system point to subsequent dynamical evolution following the birth of its planets.

In a statistical analysis of orbital parameters, Dawson & Murray-Clay (2013) found that the orbits of giant planets around metal-rich stars were more eccentric than around metal-poor stars. They postulate that metal-rich stars have solid-rich protoplanetary disks that can form more giant planets than metal-poor stars, which could ultimately engage in gravitational interactions. Multiple planets of similar mass, formed close together, and originally in coplanar, circular orbits can gravitationally excite each other’s orbits causing them to become eccentric, misaligned, and occasionally ejecting one planet out of the system in a process called planet-planet scattering (Weidenschilling & Marzari, 1996; Chatterjee et al., 2008; Marzari et al., 2002; Raymond et al., 2010). Given the high metallicity of 14 Her A, and the similar masses, large eccentricities, and misaligned orbits of 14 Her b and c, planet-planet scattering is a likely explanation for the current configuration of the system.

An external possibility is that a stellar fly-by may have triggered gravitational interactions between the planets, causing them to scatter into more eccentric orbits. Stellar fly-bys tend to disrupt a system over time scales of a few million to a few hundred million years and could lead to the ejection of one or more planets within 100 Myr (Malmberg et al., 2011; Bailer-Jones, 2015).

A more intriguing possibility is that the system initially might have had 3 nearly equal mass giant planets in relatively close, circular, coplanar orbits. After their natal disk is depleted of gas, the planets can engage in close encounters that result in excited orbital eccentricities. In the absence of gas drag, the eccentricities are dissipated through collisions, tidal circularization in the proximity of their stars, or dynamical friction by a residual population of planetesimals (Ida et al., 2013). At moderate impact parameters, a perturber can cause wide scattering rather than a collision, with recoil velocities close to the planets’ surface escape speed. At a distance of a few AU away from the star, this kind of perturbation typically leads to one planet escaping the gravitational potential of the star (Chatterjee et al., 2008; Jurić & Tremaine, 2008). In order to reach stability, the planets that survive develop large eccentricities and widely separated semimajor axes (Ford & Rasio, 2008; Carrera et al., 2019), like in the case of the 14 Her system.

The mutual gravity of the planets would have caused the scattering of the most massive one to an eccentric, closer-in orbit (i.e., b with ${9.1}_{-1.1}^{+1.0}$ $M_{\mathrm{Jup}}$ at ${2.845}_{-0.039}^{+0.038}$ AU and $e=$${0.3686}_{-0.0031}^{+0.0032}$), the intermediate one to an eccentric, far-out orbit (i.e., c with ${6.9}_{-1.0}^{+1.7}$ $M_{\mathrm{Jup}}$ at ${27.4}_{-7.9}^{+16}$ AU and $e=$${0.64}_{-0.13}^{+0.12}$), and the ejection of the least massive one out of the system. Initially resonant orbits of 3 coplanar planets can become unstable causing planet-planet scattering, leaving behind a two planet-system with a large semimajor axial ratio ($\alpha=a_{b}/a_{c}<0.3$) with mutual inclinations of $\sim 30^{\circ}$ and up to $70^{\circ}$ (Libert & Tsiganis, 2011). With a semimajor axial ratio of 0.12 and mutual inclination of $\Theta_{bc}=$$96.3_{-36.8}^{+29.1}$ degrees, the orbital parameters of today’s 14 Her system certainly fit these criteria.

The ejected planet would have had a mass lower or equal than that for the c planet ($M_{c}=$${6.9}_{-1.0}^{+1.7}$$M_{\mathrm{Jup}}$), and given the age of the primary star ($4.6^{+3.8}_{-1.3}$ Gyr), it would have had a temperature of $\lesssim 250\,$K. Isolated, planetary-mass objects at these temperatures are routinely identified as Y dwarfs, the coldest class of brown dwarfs of stellar-like origin (e.g., Cushing et al., 2011; Luhman, 2014; Marocco et al., 2019; Bardalez Gagliuffi et al., 2020). Depending on the relative occurrence of planet-planet scattering, the temperature-defined “Y dwarf” population might be of mixed origin, with both stellar-born objects and ejected planets among their ranks.

Past direct imaging campaigns have rejected the presence of stellar or substellar companions to 14 Her with confidence. In an effort to identify companions to planet-hosting FGK stars Luhman & Jayawardhana (2002) and Patience et al. (2002) imaged 14 Her as part of their target lists with Keck and Lick adaptive optics (AO), respectively. These studies ruled out the presence of companions up to a magnitude difference of $\Delta K\geq 6.4$ mag beyond $0\farcs 7$ or 12.6 AU, roughly equivalent to $\geq 0.08\,$$M_{\mathrm{\odot}}$. Carson et al. (2009) further rejected $K_{s}=18$ mag companions at a $5\farcs 0$ separation with Palomar AO. Rodigas et al. (2011) used the MMT AO system for deep imaging of 14 Her in the $L^{\prime}$-band as part of a larger survey to set direct imaging constraints on radial velocity planets.

However, no previous imaging has been able to resolve either planet due to their intrinsic faintness and large contrast with their host star. Based on our derived masses and the age of the star, we estimated effective temperatures and contrasts for the 14 Her planets for a suite of cloudless, hot start evolutionary models (Saumon & Marley, 2008; Baraffe et al., 2003). We estimate a $T_{\mathrm{eff}}$$=290-300$ K for the b planet and $T_{\mathrm{eff}}$$=260$ K for the c planet. These cold temperatures in turn imply extreme faintness in NIR bands, leading to contrasts of the order of $10^{-9}-10^{-10}$, far beyond the capabilities of current instrumentation (Table 3). Therefore, it is no surprise that these planets have eluded direct detection to date.

We also estimated the reflected light fraction for both planets by calculating the fraction of light emitted by the star that is intersected by the planet at a given distance and then reflected, based on its global atmospheric properties and phase:

where $R_{p}=1\,R_{\mathrm{Jup}}$ is the radius of the planet, $a_{p}$ is the semimajor axis of the planet, $A_{B}$ is the Bond albedo of the planet, which we approximate as Jupiter’s value ($A_{B}=0.503\pm 0.012$; Li et al. 2018), and $\Phi(\alpha)$ is the Lambert phase function at a given angle $\alpha$ (Madhusudhan & Burrows, 2012). The angle $\alpha$ is defined as the angle between the observer, planet, and star with its vertex on the planet. For simplicity, we assume an achromatic phase. We evaluated the reflected fractions when the planets were at quadrature ($\alpha=\pi/2$) (Madhusudhan & Burrows, 2012).

The $T_{\mathrm{eff}}$ of these planets rival some of the coldest known brown dwarfs (e.g., WISE J0830+2837, Bardalez Gagliuffi et al. 2020; WISE J0855$-$0714, Luhman 2014), which are brightest in the mid-infrared (Morley et al., 2012; Cushing et al., 2011). Upcoming facilities such as the Near Infrared Camera (NIRCam) aboard the James Webb Space Telescope (JWST) are ideally suited to detect objects of these temperatures at high sensitivity, although the contrast with the starlight may impact the detection of planets like 14 Her b or c (Table 3). While the b planet has a reflected contrast ratio in the order of $10^{-9}$ and could be potentially detectable with the Nancy Grace Roman Space Telescope, its angular proximity to the star ($a_{c}=158.9^{+2.1}_{-2.2}$ mas) could prove challenging for the Coronagraph Instrument. The c planet has a reflected contrast ratio in the order of $10^{-11}$, so even despite its larger angular separation from the star ($a_{c}=1529_{-442}^{+869}$ mas), it is too faint to be detected in optical wavelengths.