14 Her: a likely case of planet-planet scattering

In this Letter, we measure the full orbital architecture of the two-planet system around the nearby K0 dwarf 14 Herculis. 14 Her (HD 145675, HIP 79248) is a middle-aged ($4.6^{+3.8}_{-1.3}$ Gyr) K0 star with two eccentric giant planets identified in the literature from radial velocity (RV) variability and long-term trends. Using archival RV data from Keck/HIRES in concert with \textit{Gaia-Hipparcos} acceleration in the proper motion vector for the star, we have disentangled the mass and inclination of the b planet to ${9.1}_{-1.1}^{+1.0}$ $M_\mathrm{Jup}$ and ${32.7}_{-3.2}^{+5.3}$ degrees. Despite only partial phase coverage for the c planet's orbit, we are able to constrain its mass and orbital parameters as well to ${6.9}_{-1.0}^{+1.7}$ $M_\mathrm{Jup}$ and ${101}_{-33}^{+31}$ degrees. We find that coplanarity of the b and c orbits is strongly disfavored. Combined with the age of the system and the comparable masses of its planets, this suggests that planet-planet scattering may be responsible for the current configuration of the system.


INTRODUCTION
The orbital parameters of a planetary system are sculpted by its formation processes and evolutionary history in a quest for a stable configuration. Circular, coplanar planetary orbits are a natural consequence of formation in a disk, whether by core accretion or gravitational instability (e. g., Lissauer 1993;Boss 2001). This effect is observed in the low eccentricities and low mutual inclinations of the Solar System planets, and multi-exoplanet systems (e. g., Limbach & Turner 2015). ALMA observations of protoplanetary disks with gaps also suggest coplanarity likely as a result 2 Bardalez Gagliuffi et al.
of newborn planets clearing gas and dust from the disk (e.g., HL Tau Dipierro et al. 2015).
The dynamical evolution of an early planetary system determines the final configuration of its components (Chatterjee et al. 2008;Carrera et al. 2019;. Fly-by events can excite orbital eccentricities and even eject planets on timescales proportional to the impact parameter of the fly-by (e.g., Malmberg et al. 2011). Fly-bys can also trigger planet-planet scattering, decreasing the semimajor axes of some planets (typically the more massive ones) while leaving others in much wider orbits (Scharf & Menou 2009;Veras et al. 2009). The orbital architectures of mature systems reflect both their formation conditions and their dynamical evolution.
In this Letter, we present the full orbital architecture of the 14 Her system. Over 20 years of RV follow-up show a clear signature of a giant planet (e.g. Butler et al. 2003;Naef et al. 2004;Goździewski et al. 2006), and a long-term trend for a second one (e. g., Hirsch et al. 2021;Rosenthal et al. 2021;Wittenmyer et al. 2007). With absolute astrometry from Hipparcos and Gaia, we are able to break the M sin i degeneracy for both planets and calculate their orbital parameters, as well as identify strong evidence for a high mutual inclination between the orbits. The orbits of the 14 Her planets hint at a turbulent past marked by dynamical interactions between the two massive planets that led to their current, peculiar configurations.
We structure the Letter as follows. In Section 2, we present the existing stellar characteristics from the literature. In Section 3, we describe our RV and absolute astrometric data. Section 4 describes the orvara tool for orbit fitting and our analysis procedure. In Section 5, we report the orbital parameters derived for the two planets of this system. In Section 6, we explore the implications of our results for the formation and evolution of the 14 Her system. We present our conclusions in Section 7.

SYSTEM CHARACTERIZATION
14 Her (α = 16:10:24.315, δ = +43:49:03.50) is a middle-aged K0 dwarf located 17.9416 ± 0.0072 pc away with an estimated T eff = 5282 K (Gaia Collaboration et Lindegren et al. 2021). Due to its brightness and proximity to Earth, 14 Her was one of the first stars monitored with RV to search for exoplanets (e. g., Perrier et al. 2003;Naef et al. 2003Naef et al. , 2004. Even though significant RV variations had been found in the ELODIE RV data by 1996, the first formal discovery publication for 14 Her b was not presented until 2003 (Butler et al. 2003). 14 Her c was not identified until 2007 (Wittenmyer et al. 2007) due to its long orbital period.
Similarly to other planet-hosting stars discovered around the same time (e.g., ρ 1 Cnc), 14 Her is "super"-metal-rich with a [Fe/H] = 0.50 ± 0.05 (Taylor 1996;Gonzalez et al. 1999 We infer the fundamental parameters of 14 Her using the Bayesian activity-age dating method of Brandt et al. (2014), the luminosity, effective temperature, and angular diameter relations of Casagrande et al. (2010), and the PARSEC isochrones (Bressan et al. 2012). We find an age of 4.6 +3.8 −1.3 Gyr based on the star's low chromospheric activity and slow rotation period. We use the V T − J color from Tycho-2 (Høg et al. 2000) and 2MASS (Cutri et al. 2003) and adopt a metallicity of [Fe/H] = 0.43±0.07 to infer a luminosity of 0.67 ± 0.02 L using the relations given in Casagrande et al. (2010). Our adopted metallicity range spans 2/3 of the measurements in the PASTEL catalog ), all of which are from high-dispersion, high signal-to-noise spectroscopy. We do require a slight extrapolation of the Casagrande et al. (2010) relations, which are only validated to [Fe/H] = 0.4. The effective temperature relations of Casagrande et al. (2010), based on the V T − K s color, give T eff = 5310 ± 30 K after including both measurement errors and the 18 K rms scatter about their calibrated relation. These combine with the measured Gaia parallax to give a radius of 0.97 ± 0.02 R . This compares well with the J-band-based angular radius of 0.99 ± 0.02 R using the V T − K s color (Table 6, Casagrande et al. 2010). Finally, we use our inferred age distribution, a uniform [Fe/H] distribution, and the Salpeter initial mass function as priors, and construct a posterior mass distribution using the PARSEC isochrones together with our inferred metallicity and luminosity. We finally obtain a mass for 14 Her A of 0.98 ± 0.04 M .
Previous estimations of the rotation period of this star place it at 22.38 days (Sissa et al. 2016) and 48.5 ± 1.137 days (Watson et al. 2010), both of these coming from activity-period relations based on logR HK measurements. While 14 Her has been observed in TESS sectors 24 and 25, we are limited by the observing window of a given sector of 27.4 days, and we do not attempt to combine TESS sectors since is not trivial. From an archival ASAS-SN light curve, spanning just over 5 years (2013-02-13 to 2018-09-06), we have obtained a significant period of 29.5 days with a Lomb-Scargle periodogram (Figure 1), which is consistent with the age of the star. 3. DATA

Keck/HIRES radial velocity
We have collected 283 archival RV data points from Keck/HIRES Butler et al. 2003;Wittenmyer et al. 2007;Hirsch et al. 2021) spanning over 20 years (1997-04-07 to 2020Rosenthal 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 µm providing spectral resolutions between R∼ 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 planethosting 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.

Hipparcos and Gaia absolute stellar astrometry
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σ. 14 Her has a mean acceleration in its proper motion direction of (a α , a δ ) = (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.

THREE-BODY ORBIT FIT WITH ORVARA
In order to fully determine the orbits of the 14 Her planets, we use the RV and absolute astrometry data in concert. The Orbits from Radial Velocity, Absolute, and/or Relative Astrometry Python package (orvara; Brandt et al. 2021b) works by fitting Keplerian orbits to a combination of absolute astrometry from Hipparcos-Gaia, RV, and relative astrometry if available. For 3-body orbits, as is the case for 14 Her, the star's motion is approximated as a superposition of one Keplerian orbit due to each companion (see Brandt et al. 2021a for a demonstration and validation of the approach).
We use orvara to infer masses and orbital parameters for both 14 Her planets simultaneously. We use a parallel-tempered MCMC sampler (Foreman-Mackey et al. 2013;Vousden et al. 2021) to robustly explore the parameter space with several copies of the system, randomly initialized at 30 temperatures. For each temperature, we use 100 walkers with 600,000 steps per walker, discarding the first 30,000 steps as burn-in and saving every 50 th step which are then used for inference. For our final values, we combined four such chains into a long pseudo-chain. We set a Gaussian prior on the primary mass of 0.98 ± 0.04 M and a log-flat prior on the RV jitter to range between 10 −5 − 10 m/s, letting the jitter of the pre-and post-2004 HIRES upgrade vary independently. We test the convergence of the MCMC algorithm by running several 600,000 step chains and confirming a stable plateau in log likelihood for both planet solutions.
The misaligned planets of 14 Her 7 5. RESULTS

Orbital parameters of 14 Her b and c
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.0 −1.1 M Jup and its orbit has a semimajor axis of a =2.845 +0.038 −0.039 AU with a moderate eccentricity of e =0.3686 +0.0032 −0.0031 and a 32.7 +5.3 −3.2 degree inclination with respect to the plane of the sky. Its orbital period translates to P =4.8277 +0.0022 −0.0023 years. Despite only seeing a long-term trend in the RV data covering ∼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.7 −1.0 M Jup ) and it is located much farther away from the star, at a =27.4 +16 −7.9 AU. 14 Her c has a highly eccentric orbit (e =0.64 +0.12 −0.13 ). With respect to the orbit of 14 Her b, the c planet has a broad distribution of inclinations (i =101 +31 −33 degree) but hints at being misaligned (see Figure 4). Its orbital period amounts to P =144 +139 −58 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

Relative orientation of 14 Her b and c
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 Ω 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 8 Bardalez Gagliuffi et al. Calculating this expression on each chain of our MCMC sampling returns an angle of Θ bc =96.3 +29.1 −36.8 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 ( In comparison with π 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 π Men. Another key measurement to constrain the orientation of the system is the spin-axis alignment of the star.  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. π Men, Huber et al. 2021; α 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.

Radial velocity offset in HIRES between pre-and post-upgrade
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  suitable zero-points and jitter values independently (see Table 2). Comparing the zero-points before and after the upgrade, we find a difference of ∆ 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 postupgrade 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 +16 −8.8 AU), the eccentricity (e c = 0.69 +0.10 −0.11 ), and the mass (M c = 7.11 +1.5 −0.82 M 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 √ 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 ∼ 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 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.0 −1.1 M Jup at 2.845 +0.038 −0.039 AU and e =0.3686 +0.0032 −0.0031 ), the intermediate one to an eccentric, far-out orbit (i.e., c with 6.9 +1.7 −1.0 M Jup at 27.4 +16 −7.9 AU and e =0.64 +0.12 −0.13 ), 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 (α = a b /a c < 0.3) with mutual inclinations of ∼ 30 • and up to 70 • (Libert & Tsiganis 2011). With a semimajor axial ratio of 0.12 and mutual inclination of Θ bc =96.3 +29.1 −36.8 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.7 −1.0 M Jup ), and given the age of the primary star (4.6 +3.8 −1.3 Gyr), it would have had a temperature of 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.

Potential as a future extreme-AO target
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 planethosting 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 ∆K ≥ 6.4 mag beyond 0. 7 or 12.6 AU, roughly equivalent to ≥ 0.08 M . Carson et al. (2009) further rejected K s = 18 mag companions at a 5. 0 separation with Palomar AO. Rodigas et al. (2011) used the MMT AO system for deep imaging of 14 Her in the L -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 eff = 290 − 300 K for the b planet and T 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 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 ± 0.012; Li et al. 2018), and Φ(α) is the Lambert phase function at a given angle α (Madhusudhan & Burrows 2012). The angle α 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 (α = π/2) (Madhusudhan & Burrows 2012).
The T 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 · · · · · · · · · · · · · · · · · · 1.16E-11 · · · · · · Baraffe et al. 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 +869 −442 mas), it is too faint to be detected in optical wavelengths.

CONCLUSIONS
In this paper we have characterized the orbital parameters and dynamical evolution of the 14 Her planetary system. Using orvara, which combines RV and astrometric accelerations, we have obtained a dynamical mass of 9.1 +1.0 −1.1 M Jup and an inclination of 32.7 +5.3 −3.2 degrees for the b planet, hence disentangling the M sin i degeneracy for this object for the first time. We also set dynamical mass and orbital constraints on the c planet, albeit with larger uncertainties. We have also characterized the fundamental parameters of the star in order to study this system as an ensemble.
Our results describe a middle-aged K0 star with two massive planets in highly eccentric, misaligned orbits. The mutual orientation between the b and c orbits is Θ bc = 96.3 +29.1 −36.8 degrees. Coplanarity is disfavored for this system, a fact that combined with the large eccentricities, suggests a disruptive planet-planet scattering event leading to the current architecture. An N-body dynamical simulation could strengthen the hypothesis that a third 7 M Jup planet was ejected from the system.
Based on the age of the star and the dynamical planetary masses derived, we infer the effective temperature of the planets from hot start evolutionary models to be