Revisiting Orbital Evolution in HAT-P-2 b and Confirmation of HAT-P-2 c

One possible formation mechanism for Hot Jupiters is that high-eccentricity gas giants experience tidal interactions with their host star that cause them to lose orbital energy and migrate inwards. We study these types of tidal interactions in an eccentric Hot Jupiter called HAT-P-2 b, which is a system where a long-period companion has been suggested, and hints of orbital evolution (de Wit et al. 2017) were detected. Using five additional years of radial velocity (RV) measurements, we further investigate these phenomena. We investigated the long-period companion by jointly fitting RVs and Hipparcos-Gaia astrometry and confirmed this long-period companion, significantly narrowed down the range of possible periods ($P_2 = 8500_{-1500}^{+2600}$ days), and determined that it must be a substellar object ($10.7_{-2.2}^{+5.2}$ $M_j$). We also developed a modular pipeline to simultaneously model rapid orbital evolution and the long-period companion. We find that the rate and significance of evolution are highly dependent on the long-period companion modeling choices. In some cases the orbital rates of change reached $de/dt = {3.28}_{-1.72}^{+1.75} \cdot 10^{-3}$/year, $d\omega/dt = 1.12 \pm 0.22 ^{\circ}$/year which corresponds to a $\sim 321$ year apsidal precession period. In other cases, the data is consistent with $de/dt = 7.67 \pm 18.6 \cdot 10^{-4}$/year, $d\omega/dt = 0.76\pm 0.24 ^{\circ}$/year. The most rapid changes found are significantly larger than the expected relativistic precession rate and could be caused by transient tidal planet-star interactions. To definitively determine the magnitude and significance of potential orbital evolution in HAT-P-2 b, we recommend further monitoring with RVs and precise transit and eclipse timings.


INTRODUCTION
The role of planet-star interactions in the migration of Hot Jupiters is not well understood. The original discovery of Jupiter-sized planets on short-period orbits (Mayor & Queloz 1995) challenged our theories of solar system formation and it was proposed that Hot Jupiters must have migrated inward after forming beyond the ice line (Lin et al. 1996). One possible migration mechanism proposed is that high-eccentricity gas giants experience tidal interactions with their host star that cause them to lose orbital energy and migrate to a close-in orbit (Rasio & Ford 1996;Weidenschilling & Marzari 1996;Lin & Ida 1997).
High-eccentricity tidal migration theories often focus on the tides within the planet. However, equilibrium tides (Hut 1981;Eggleton et al. 1998) in the star are accompanied by dynamical tides (Goodman & Dickson 1998;Fuller & Lai 2012) consisting of resonantly excited g-modes, which may result in energy and angular momentum exchanges between the star and orbit that can be several orders of magnitude larger (Witte & Savonije 2002;Ma & Fuller 2021). The resulting changes in the orbit could impact radial velocity and photometric ob-servations of a planetary system. Hints of rapid orbital evolution were detected in the HAT-P-2 system where other signs of strong planet-star interactions have been observed (de Wit et al. 2017).
HAT-P-2 is an F8 type star and hosts an 8.70 +0.19 −0.20 M j mass planet on an eccentric (e = 0.51023 ± 0.00042), short period (P = 5.6334675 ± 0.0000013 days) orbit. Due to this combination of characteristics, it has been identified as an ideal target for studying planet-star interactions (Jordán & Bakos 2008;Fabrycky & Winn 2009;Levrard et al. 2009;Cowan & Agol 2011;Lewis et al. 2013Lewis et al. , 2014Salz et al. 2016). Now, we are revisiting this interesting system because the California Planet Search (CPS) team observed HAT-P-2 for an additional five years with the HIRES spectrograph, which extends the RV baseline from nine years to fourteen years. In addition, we also merged RVs from multiple instruments and increased the number of RVs from 54 to 103 compared to de Wit et al. (2017). With these additional RV measurements, we were able to further investigate whether the hints of rapid orbital evolution appear to persist and can be confirmed. We developed our own modular orbital evolution pipeline that allows for an apples-to-apples comparison between various orbital models.
In addition, we performed an in-depth assessment of the possibility of an outer long-period companion as has been suggested by Lewis et al. (2013); Knutson et al. (2014); Bonomo et al. (2017); Ment et al. (2018). We combined our RV measurements with astrometric observations and confirmed a persistent signal that caused by a substellar long-period companion HAT-P-2 c.
Our paper is organized as follows. In Section 2, we describe the RV and astrometric data included in this analysis, and in Section 3, we describe RV orbital fitting pipeline and how we performed a joint astrometric-RV fit. In Section 4, we describe our results and to evalulate our models. In Section 5, we describe the interpretation of these results. Lastly, in Section 6 and 7, we provide directions for future work and we conclude.

Precise Radial Velocity Measurements
Our analysis includes 72 spectra of HAT-P-2 obtained by the California Planet Search (CPS) Team using the High Resolution Echelle Spectrograph (HIRES) on the W.M. Keck Observatory 10m telescope Keck-I (Vogt et al. 1994). The standard CPS data reduction pipeline was used, which uses an iodine cell as a wavelength reference as further described in Howard et al. (2010). These 72 RVs span 14 years (2006-09-03 through 2020-08-31).
In addition to our 72 HIRES RVs, we also included 31 publicly available RVs from the HARPS-N spectrograph at the 3.6m Telescopio Nazionale Galileo on La Palma in the Canary Islands. HARPS-N is a temperature and pressure stabilized echelle spectrograph (Cosentino et al. 2012). HARPS-N spans the wavelength range from 383 to 693 nm and has a resolving power of λ/∆λ = 115,000. These RVs span 2013-03-12 through 2017-09-28. We included the published RVs from Bonomo et al. (2017) and downloaded additional publicly available RVs through the TNG Archive 1 .
We removed RV measurements that occur during transit (Winn et al. 2007) such that we are not using RVs that are affected by the Rossiter-McLaughlin effect. The RVs are included in the appendix in Table A.4.

Astrometric Measurements
The astrometric data used in this work are derived from the Hipparcos-Gaia Catalogue of Accelerations (HGCA; Brandt 2018Brandt , 2021. This consists of three independent proper motion measurements: the Hipparcos proper motion (µ H ), the Gaia DR3 proper motion (µ G ), and the mean Hipparcos-Gaia proper motion (µ HG ). See Brandt (2018) and Brandt et al. (2019) for further details. The astrometric data has a time baseline of approximately 25 years .
RVs provide us with the minimum mass M p sin(i) of an object, where i is the inclination of the system. However, when combining RVs with astrometric measurements, we can constrain the inclination of the system and derive the actual mass (M p ). In this way, astrometric measurements can help constrain the nature of the outer companion of HAT-P-2 (stellar, brown dwarf, planet).

Radial Velocity Orbital fitting pipeline
We designed a modular RV orbital fitting pipeline that allows for comparison of various dynamical scenarios to determine the origins of the RV variation seen in de Wit et al. (2017). We use radvel's kepler.rv drive function (Fulton et al. 2018) to solve Kepler's equation in our pipeline and build upon this software by allowing orbital parameters to change over time. Our modular pipeline is parallelized to speed up orbital fitting and can easily be applied to other RV datasets. Our pipeline can notably perform fits of the following sta-ble orbital configurations: (i) a stable, one-planet Keplerian fit, (ii) a stable, two-planet Keplerian fit, (iii) a stable, two-planet fit where the long-period companion is modeled as a quadratic trend. In addition, our pipeline can model the eccentricity (e) and argument of periastron (ω) as linearly evolving parameters such that e = de/dt · time + e 0 and ω = dω/dt · time + ω 0 . We refer to models where e and ω can linearly evolve with time as first-order, evolving models. In particular, our pipeline can perform the following evolving orbital configurations: (iv) a first-order evolving, one-planet Keplerian fit, (v) a first-order evolving, two-planet Keplerian fit, (vi) and a first-order evolving, two-planet fit where the long-period companion is modeled as a quadratic trend.
In case (i), the model can be described by 13 free parameters, 5 of which describe the orbit (period (P ), planet mass (M p ), time of periastrion passage (t p ), eccentricity (e) and argument of periastron (ω) which are parameterized as √ e sin ω and √ e cos ω), 2 of which describe a linear (γ) and quadratic (γ) trend with time 2 and 6 which correspond to the offsets (γ) and RV jitter terms for each of the RV datasets. RV jitter describes the instrumental and astrophysical stochastic signals in the data. Astrophysical signals originate from stellar variability such as inhomogeneties on the stellar surface, pressure-mode oscillations, and granulation. Since each instrument has its own noise characteristics, it is common to fit for a jitter term for each instrument. We note that RV models are commonly parametrized using semi-amplitude (K) rather than mass. However, when e and ω are allowed to vary over time for an evolving planet model, this parametrization allows mass to vary over time, which is unphysical. Parametrizing the RV model in terms of mass prevents this and instead only allows K to vary as a function of e and ω. In case (ii), we gain an additional 5 orbital parameters for a secondary companion, which sums to 18 free parameters 3 . In case (iii), the secondary companion is modeled as a second-order background trend so this adds a linear and quadratic term and thereby two additional free parameters. In case (iv), we add first-order polynomials for e and/or ω and this adds one free parameter per polynomial. In case (v), we add an additional five parameters to fit for the orbital parameters of a secondary companion. Lastly, in case (vi), we instead model the secondary companion as a linear and quadratic trend so this adds two free parameters.
Once a given orbital configuration is chosen, we sample our parameter space using edmcmc 4 , which is an implementation of Markov Chain Monte Carlo (MCMC) that incorporates differential evolution (DE;Ter Braak 2006). DE is a genetic algorithm that solves the problem of choosing an appropriate scale and orientation for the jumping distribution, which significantly speeds up the fitting process.
After performing the fit, the pipeline produces several diagnostic plots and performance metrics (described in Section 3.2) to assess goodness of fit. Convergence of the parameters is determined using the Gelman-Rubin statistic (Gelman & Rubin 1992) where our Gelman-Rubin threshold < 1.02. We computed 10, 000 MCMC chains and discarded the first 25% burn-in samples to produce our posterior samples. We implemented Gaussian priors to the P and t p parameters based on transit and occultation times from previous studies (Ivshina & Winn 2022). For our models that include a Keplerian long-period companion, we also apply Gaussian priors to the companion's period (P 2 ) and time of periastron passage (t p2 ) based on the astrometry analysis described in Section 3.3 since the RVs do not sufficiently span its orbital period to constrain these parameters alone. To derive our final parameter values, we computed the median and 68.3% confidence intervals. The results for our comparison of models (i)-(vi) are described in Section 4.

Performance metrics
We use four different metrics to assess the goodness of fit of our various orbital models: (1) the loglikelihood, (2) the reduced chi-squared statistic (χ 2 r ), (3) Bayesian information criterion (BIC; Schwarz 1978), (4) and the Akaike information criterion (AIC; Akaike 1974). Each of these metrics have different advantages and emphasize different features of the fitted models. Generally, each metric penalizes the number of free parameters in a different way and the AIC tends to be the least punitive of more complex models. In our analysis, these different metrics all provide consistent results. Each metric is described in additional detail in the appendix.

Joint fit of astrometry and Radial Velocities
In order to further investigate the possibility of a longperiod companion to HAT-P-2 b, we performed a joint fit of the RV and astrometric data sets. The model used Figure 1. Comparison of our stable, non-evolving orbital models. From top to bottom, we include the stable one-planet model with no background trends (a, b), the stable two-planet model with the long-period companion modeled as a Keplerian (c,d), the stable two-planet model with the long-period companion modeled as a quadratic trend (e,f). For more details on these orbital models, see Section 3.1. In each left hand panel (a, c, e), we plot the orbital fit of HAT-P-2 b minus a secondary companion (or no companion if this was not included in the model). The data is color-coded by instrument (turquoise = HARPS-N 2013-2015, red = HARPS-N 2015-2017, black = HIRES). The median model is plotted in dark purple and a 100 draws from the posteriors are plotted in light purple. For each right hand panel (b, d, f), we plot the secondary companion (or no trend/companion if this was not included in the model) minus the median HAT-P-2 b model. We also include the goodness-of-fit metrics corresponding to each orbital model in the yellow box. Red. χ 2 r is an abbreviation for reduced χ 2 r . to perform this joint analysis is described in detail in Venner et al. (2021). In summary, we jointly fit the two datasets with a three-body Keplerian model with 22 variable parameters (the stellar mass and parallax, which were assigned Gaussian priors of 1.33 ± 0.03 M ⊙ and 7.781 ± 0.012 mas respectively; P , K, e, ω, t p for both companions; orbital inclination (i) and longitude of ascending node (Ω) for the outer companion; proper motion of system barycentre parameters (µ bary,RA , µ bary,Dec ); and normalization and jitter parameters for each RV dataset). In comparison to the models described in Section 3.1, this is equivalent to a stable two-planet Keplerian fit.
Cursory examination of the model reveals a degeneracy between P and e for the outer companion, arising due to incomplete coverage of the orbit. If the outer planet has a shorter orbital period, it is more likely to be observed during its periastron passage than for longer orbital periods. We account for this in the same way as Blunt et al. (2019) and Venner et al. (2022) by adopting an informed prior on the orbital period such that: where P is the probability of observing periastron passage, P is the orbital period, t d is the duration of periastron passage, and B is the baseline of observations. This prior is analogous to the period priors used for long-period transiting exoplanets (e.g. Vanderburg et al. 2016;Kipping 2018). Imposing this prior slightly penalizes orbital solutions where the periastron passage of a long orbit happens to occur during the comparatively short span of observations. This discourages our model from fine-tuned long-period orbital solutions for the outer companion, without excluding them entirely. Blunt et al. (2019) tested a range of values for t d and found that their orbital fits were indistinguishable, and thus chose t d = 0. We do the same in this analysis.
We again use edmcmc to explore the parameter space of the model. We produce posterior samples by discarding the first 25% of the chain as burn-in and save every hundredth step for each walker. We extract the median and 68.3% confidence intervals for the parameters from the samples and compute the same confidence intervals for physical parameters that can be derived (i.e. companion mass).

Comparing Radial Velocity Orbital Fits
For each of our four models, we computed the performance metrics described in Section 3.2 and listed them in Tables 1 and 2. First, we can compare the stable, one-planet model with the stable model that allows for a long-period companion modeled as either a Keplerian or a quadratic trend. We find that the stable, twoplanet model (2 Keplerians) has a significantly lower reduced χ 2 r = 1.035, lower ∆BIC= 14.90, lower ∆AIC= 12.26, and higher loglikelihood= −432.05 compared to the one-planet model (χ 2 r = 2.806, ∆BIC= 165.18, ∆AIC= 175.72, loglikelihood= −518.7). Furthermore, the stable, two-planet model (Keplerian + Quadratic) has a similar reduced χ 2 r = 1.045, similar ∆BIC= 21.16, ∆AIC= 15.90, and similar loglikelihood= −432.87 compared to the other two-planet model. To visually examine these orbital models, we plot the median fit along with 100 random draws for all three of these models in Figure 1a-f. Overall, the two-planet models (either the 2 Keplerian or Keplerian + Quadratic version) are favoured in the comparison of these stable orbit models.
Next, we can compare various evolving orbit scenarios for HAT-P-2 b. We include the one-planet (Keplerian), two-planet (2 Keplerians), and two-planet (Keplerian + Quadratic) scenario where the HAT-P-2b's e and ω can vary linearly with time. We find that the evolving, one-planet model (χ 2 r = 2.504, ∆BIC= 142.47, ∆AIC= 147.74, loglikelihood = −502.788) provides a marginally better fit to the observations than a stable one-planet model. This is not a significant difference however. The evolving one-planet model is plotted in Figure 2a, b. Next, we find that model (v), the evolving 2 planets (2 Keplerians) model provides a slightly better fit (χ 2 r = 0.963, ∆BIC= 10.181, ∆AIC= 7.546, loglikelihood = −429.691) than model (ii), the stable 2 planets (2 Keplerians) model. Overall, we find that model (vi), the evolving, two-planet model (Keplerian + Quadratic) provides the best fit with the lowest χ 2 r = 0.891, the lowest BIC= 909.453, the lowest AIC= 877.836, and the highest loglikelihood= −426.918. We plot this model in Figure 2e,f. We note that we also performed some additional tests to see how the rate of evolution for e, ω and their significance are impacted by how the long-term companion is modeled and include this in the discussion (Section 5.1.4).

Astrometric-Radial Velocity joint fit
The results of the previous section strongly suggest the presence of an outer companion in the system, as has been found by past authors (Lewis et al. 2013;Knutson et al. 2014;Bonomo et al. 2017;Ment et al. 2018). Furthermore, our 5-year extension of the HIRES timeseries allows us to clearly observe non-linearity in the RV acceleration. Previously, Knutson et al. (2014)  and the semi-major axis to 4 − 31 AU based on the nondetection of a stellar companion in direct imaging. The strong orbital motion suggests that the companion has a relatively short period, and hence a mass on the lower end of this range. HAT-P-2 has been observed by both Hipparcos and Gaia, allowing us to make use of Hipparcos-Gaia astrometry (Brandt 2018(Brandt , 2021 to place additional constraints on the outer companion. We expect that even the lowest-mass solutions should produce a detectable (≳0.2 mas yr −1 ) astrometric signal based on the RV signal, but surprisingly we find that the Hipparcos-Gaia astrometry of HAT-P-2 is consistent with zero net acceleration. A possible explanation for this is that the orbital period is close to the 25-year astrometric timespan. Specifically, if the orbital inclination is close to edge-on and the outer companion is at conjunction during both the Hipparcos and Gaia observations, then there will be no net acceleration between the Gaia proper motion and Hipparcos-Gaia mean proper motion, as observed. This is borne out by our joint fit to the RVs and astrometry, the results of which we plot in Figure 3. As anticipated, we find an orbital period of the outer companion of 23.2 +7.1 −4.1 yr, comparable to 25 year timespan of the astrometric data. To demonstrate the importance of the astrometry for this result, in Figure 4 we compare the period distribution from our joint fit with that of an RV-only model. The astrometric non-detection allows us to exclude orbital periods significantly longer than 15000 days at 95% confidence. Furthermore, the non-detection constrains the orbital inclination to ≈90 degrees as this is the only orbital configuration consistent with zero net astrometric acceleration over the span of observations. The key posterior parameters of the outer companion are as follows: −13 m/s e c = 0.37 +0.13 −0.12 ω c = 38 ± 32 degrees T p,c = 2454150 +430 −800 BJD i c = 90 ± 16 degrees. (2) The longitude of node (Ω) is unconstrained. With the orbital period and inclination resolved, we find that the mass of this second companion is 10.7 +5.2 −2.2 M j . The mass posterior is entirely substellar and corresponds to a planetary-mass object, comparable in mass to HAT-P-2 b, which straddles the deuterium-burning limit at 13 M j . We hence designate this companion as HAT-P-2 c, confirming the presence of a second planetary-mass companion in the HAT-P-2 system. We note that the parameters of HAT-P-2 c derived from this joint astrometric-RV fit agree within error with the RV fits in Tables 1, 2. The parameter values for HAT-P-2 c from the joint fit should be considered more robust.

Rates of orbital evolution
From our orbital fitting analyses, we find that although a two-planet 5 , evolving orbit model that allows e and ω to vary linearly with time fits the RV measurements best, the orbital parameters of HAT-P-2 c are poorly constrained in these RV models and follow- Figure 3. The orbit of HAT-P-2 c from a joint fit to the RVs (left) and astrometry (right). Colors are as in Figure 1. The astrometry is consistent with zero net acceleration; the limits of this non-detection allow us to place useful constraints on the orbital period and inclination (Pc = 8500 +2600 −1500 days, ic = 90 ± 16 degrees). The resulting mass of HAT-P-2 c is 10.7 +5.2 −2.2 Mj, establishing it as a planetary-mass object. We indicate the next observing window (April 1st, 2023 -September 28, 2023) with green shading. Further RV observations would help to improve constraints on the companion's orbit. up observations are required to fully confirm or rule out orbital evolution in HAT-P-2 b. For this model, we find that dω/dt = 0.76 +0.24• −0.24 /year and de/dt = 7.67 +13.2 −18.5 ·10 −4 /year which correspond to 3.23σ and 0.4σ detections respectively. In Figure 5, we illustrate how these changes in ω and e correspond to changes in the shape of the orbit over time. Although these rates of change agree within error with the rates (dω/dt = 0.91 ± 0.31 • /year, de/dt = 8.9 ± 2.8·10 −4 /year) found in de Wit et al. (2017), they found a slightly more signif-icant change in e. Additional RV measurements in the next few years should allow for a further refinement of de/dt and dω/dt as described in 5.1.6. Specific choices in how HAT-P-2 c's long-term signal was modeled may also offer an explanation for the slight differences in evolution rates found in these different analyses as described in detail in Sections 5.1.2 and 5.1.4. et al. (2017) In de Wit et al. (2017), the long-term trend corresponding the outer companion and the orbital evolution were modeled differently than in this analysis. In particular, the data available to de Wit et al. (2017) only sampled the linearly decreasing component of the longterm trend and therefore de Wit et al. (2017) modeled the outer companion as a linear trend rather than a quadratic or Keplerian. Now, with an additional five years of RV data, we were able to model the outer companion as a quadratic or Keplerian instead. For completeness, we also included a test where we approximate the outer companion as a linear trend as described in Section 5.1.4.

Comparison to Radial Velocity analysis in de Wit
In addition to modeling the outer companion more accurately in our analysis, we also were able to model the orbital evolution in a more robust way with the availability of more RVs. de Wit et al. (2017) checked for changes in e and ω by i) splitting the RV data into two approxi-  . Variation in orbital shape over time for HAT-P-2 as e and ω increase linearly over time. The long-period companion is modeled as a quadratic trend here. The trends in e and ω were fit using the entire dataset, but for visualization, the data is split into four equal-sized chunks (∼ 25 observations per chunk) here. The median model corresponding to that time period is computed using the median e and ω for the corresponding time period.
mately equal halves (data prior and after BJD 2454604) and ii) fitting each dataset independently for an average value of e and omega. de Wit et al. (2017) then computed an estimate for the trends in e and omega based on the difference in values for each of these halves of the data. In this analysis, we instead fit the entire RV dataset simultaneously and define e = de/dt · time + e 0 , ω = dω/dt · time + ω 0 such that we fit for de/dt and dω/dt directly. This approach should require that the trend is more persistent across the entire RV baseline and provide more robust estimates of the uncertainties on the trends as well.

5.1.3.
Can the change in ω be explained by general relativity alone?
We computed the rate of change expected from general relativity (GR) wherê where n is the Keplerian mean motion, a is the semimajor axis, and M * is the mass of the star. We find that for HAT-P-2 b,ω GR = 1.185e − 2 degrees/year. However, in our best fit model, we findω = 0.76 ± 0.24 degrees/year, which is nearly 65 times the rate of change expected from GR. Thus, this type of rate of change cannot be explained by GR alone. Tidal planetstar interactions may potentially be able to explain these types of rapid changes. 5.1.4. The impact of modeling choices for planet c on the rate of evolution of e and ω Since previous analyses (de Wit et al. 2017;Bonomo et al. 2017) modeled the long-period companion as linear or quadratic trends, we performed additional tests where we model the second companion in the same way and investigated its impact on the derived rate of orbital evolution. In Table 3, we list the derived rates of orbital evolution and their significance of detection when not modeling the outer companion and when modeling the outer companion as a linear trend, a quadratic trend, or as a Keplerian. We find that when we do not model HAT-P-2 c, the significance of the trend in ω (5.2σ) is highest and the significance of the trend in e (1.03σ) is second-highest among these models. The derived rates of evolution are also among the highest amongst these models (de/dt = 1.75 +1.71 −1.68 ·10 −3 / year, dω/dt = 1.12 +0.22 −0.22 • /year). When we model HAT-p-2 c as a linear trend, the significance of the trend in e (1.9σ) increases and the significance of dω/dt is marginally lower but still The two best-fit models (ii, vi) are plotted side by side to compare their differences. In green, we indicate the three orbital phases we propose to sample three times to detect the orbital evolution or lack thereof. (c) The residual RVs from subtracting the expected RVs from model (vi) from model (iii). a 5.1σ trend. The rates are the highest with de/dt = 3.28 +1.72 −1.75 ·10 −3 / year and dω/dt = 1.12 +0.22 −0.22 • /year. For the quadratic model, the significance of both trends decreases (σ de/dt = 0.41, σ dω/dt = 3.23) and the rates also decrease (de/dt = 7.7 +18.6 −18.6 ·10 −4 / year, dω/dt = 0.76 +0.24 −0.24 • /year). For the Keplerian model, we find that significance in de/dt slightly increases (σ de/dt = 0.32) and the significance in dω/dt also decreases (σ dω/dt = 2.11). The rate of evolution continues to decrease for dω/dt (dω/dt = 0.228 +0.11 −0.11 • / year) and direction of de/dt changes (de/dt = −7.3 +19.0 −18.0 ·10 −4 / year). However, this trend in e is so insignificant, indicating that that parameter is unconstrained by the data and so we do not attribute a physical interpretation to this change in sign. Clearly, the exact modeling choices of HAT-P-2 c strongly affect the rates and significance of the evolution rates derived for HAT-P-2 b. The Keplerian model is the most physically motivated, but the period is still relatively unconstrained (P c = 8500 +2600 −1500 days) and so are many of the other orbital parameters. Thus, further monitoring of HAT-P-2 system to allow us to model the long-term RV changes induced by HAT-P-2 c most accurately are necessary and should allow us to constrain the rates and possible evolution of HAT-P-2 b as well. In particular, we propose follow-up RVs to constrain the orbital parameters of HAT-P-2 c and follow-up transit and eclipse observations to solve precisely for the current values of e and ω, which when combined with archival Spitzer data should allow us to definitively constrain the rate of evolution.

Spitzer Transits and Secondary Eclipses
To get additional constraints on e and ω and look for trends, We re-processed and analyzed 4.5 µm Spitzer transit and secondary eclipse measurements using the pipeline described in Berardo et al. (2019). We computed the best-fit systematic model using the pixel-level decorrelation (PLD) technique (Deming et al. 2015). We analyzed the primary and secondary eclipses by splitting the data into three chunks (2011, 2013, and 2015 observations) and performed three independent fits to obtain a measurements of e and ω for each time period and look for changes in their values. For the 2011 observations, we find e = 0.5104 ± 0.0025 and ω = 189.10 ± 2.19 • . For 2015 observations, we find e = 0.5128 ± 0.0014 and ω = 190.95 ± 1.23 • . Since the 2013 observations only included eclipses and no transits, we did not have the transit eclipse pair that is neccesary to constrain e for the 2013 window. Overall the Spitzer transit and eclipse measurements provide some constraints on the possible rates of change for e and ω, but they only contain 3-4 transits/eclipses in each subset. With this little amount of data available for each of these subsets, we do not get enough precision to accurately measure drifts in e. We plan to propose for JWST observations such that we can get the necessary precision to further constrain the rates of change of e and ω.

RV Sensitivity Analysis
In order to further characterize the HAT-P-2 system such that we can increase the significance of detection of the orbital evolution (or lack thereof), we performed a sensitivity analysis for the next observing window. In Figure 6, we plot the expected RVs for the first orbit of the spring-fall 2023 observing season (April 1st, 2023 -April 5th, 2023). The expected RVs for a stable orbit will diverge significantly from an evolving orbit model and we should be able to clearly distinguish the two scenarios as seen in Figure 6b, c. In Figure A.7, we plot the entire next observing window, which illustrates that we should also be able to distinguish clearly between the one or two-planet models and further constrain the period of a potential HAT-P-2 c. Thus, constraining the period of the long-period companion would require 3-4 RV measurements on HIRES/ Keck I. This could be done anytime in the observing window (April 1, 2023 -September 28, 2023) and would extend the RV baseline from 14 to 17 years. However, also detecting the orbital evolution would require additional RVs and would place constraints on the timing of the observations. In particular, to get a 6σ detection of the orbital evolution, we require nine new total RV measurements spanning three orbital periods (i.e. three RVs per orbit). The particular orbital phases are indicated in green in Figure 6b,c. Although taking the RV observations at these orbital phases would be required, the observations could be taken for any orbit during the observing window.

HAT-P-2 c
The presence of a long-period outer companion to HAT-P-2 has been suggested in previous studies on the basis of a significant RV acceleration (Lewis et al. 2013;Knutson et al. 2014;Bonomo et al. 2017;Ment et al. 2018). We have constrained the orbit of this companion for the first time by jointly fitting the RVs with astrometry from the Hipparcos-Gaia Catalog of Accelerations (Brandt 2018(Brandt , 2021. Though no acceleration is detected in the Hipparcos-Gaia astrometry, this non-detection is sufficiently informative to constrain the orbit of the second companion. HAT-P-2 c is a 10.7 +5.2 −2.2 M j planetary mass object with an orbital period of P c = 8500 +2600 −1500 days, placing it among the longest-period substellar companions that have ever been discovered from RVs. The constraints of the astrometric non-detection requires HAT-P-2 c to have an edge-on orbital inclination (i c = 90 ± 16 degrees), as any other configuration would produce a detectable astrometric signal. This is comparable to the result of Errico et al. (2022), who determined the true mass of the long-period planet HD 83443 c from an astrometric non-detection. HAT-P-2 c's orbital inclination is consistent with HAT-P-2 b, which is in turn has a low projected obliquity from the host star's rotational axis (Winn et al. 2007;Albrecht et al. 2012). This suggests that the orbital and rotational planes in the HAT-P-2 system are compatible with alignment, agreeing well with the conclusions of Becker et al. (2017), who found that outer planets in Hot Jupiter systems must generally have low mutual inclinations to reproduce the low obliquities observed in these systems. The absence of strong misalignments in the HAT-P-2 system appears inconsistent with a chaotic dynamical history, which may be surprising considering the remarkably high orbital eccentricity of HAT-P-2 b. The formation history of this unusual planetary system will be an interesting topic to explore in future studies.

FUTURE WORK
We have proposed for RV follow-up observations to further constrain the period and orbital parameters of HAT-P-2 c, which will help us determine the magnitude and significance of orbital evolution for HAT-P-2 b (or lack thereof). In particular, we proposed to observe during specific phases of the orbit which have the information content necessary to disentangle between the one to two-planet and evolving or non-evolving scenarios as further described in Section 5.1.6.
In addition, we plan to propose for a JWST partial phase curve of system with three goals in mind: 1. Obtain ultra-precise transit and eclipse timing measurements: these measurements will enable us to get extremely tight constraints on the current values of e and ω, and thereby constrain their evolution.
2. Measure the transient heating of HAT-P-2 b: Due to its high eccentricity and short period orbit, HAT-P-2 b is an ideal target for studying how atmospheres redistributes the time-varying heat flux from host stars. Analogous to Lewis et al. (2014), we would propose measure the transient heating and use 3D atmospheric circulation models to further our understanding of the atmospheric response of exoplanets on highly eccentric orbits.
3. Follow-up on the planet-induced oscillations observed on the star noted in de Wit et al. (2017) and check for wavelength-dependence of these oscillations. This will also allow us to construct a 2D heat map of HAT-P-2 b's surface, which was not possible in previous papers due to the complications added by the oscillations.

CONCLUSION
We developed an modular orbital evolution fitting pipeline that allowed us to perform apples-to-apples comparison between various dynamical scenarios that may explain observed radial velocity variation in the HAT-P-2 system. Using this pipeline, we confirmed that the hints of rapid change in argument of periastron (ω) and eccentricity (e) noted in de Wit et al. (2017) appear to persist with an additional 5 years of RV data. However, the significance and exact rates of this evolution is highly dependent on the modeling of the longperiod companion and requires additional follow-up observations to constrain. The largest changes in ω are significantly larger than what would be expected solely from general relativistic precession and could be instead explained by tidal planet-star interactions. Since these precession rates are not possible to explain with general relativity alone, further investigation to model the tidal planet-star interactions with stellar evolution and pulsation codes (Paxton et al. 2011;Townsend et al. 2014, Townsend et al. (GYRE;) is warranted. We also ran a joint fit with Hipparcos-Gaia astrometry and the RV measurements. This joint fit allowed us to put precise constraints on the outer companion in the system, which we find to be substellar: HAT-P-2 c. Furthermore, this work sets the stage for applying our orbital evolution fitting pipeline to a larger sample of stars that host eccentric Hot Jupiters, with the ultimate goal of modeling Hot Jupiter precession on a statistical scale.

CODE AVAILABILITY
The DE-MCMC Orbital evolution fitting pipeline developed for this manuscript is publicly available on Zenodo (de Beurs 2023) and a living version is available on github 6 . The code was run on 7 cores of an M1 Macbook Air (16 GB RAM).

APPENDIX
In this appendix, we include the radial velocity measurements included in this analysis (Table A.4). In addition, we provide a supplemental figure of the expected RV measurements for each of the possible models during the next observing window (Figure A.7). Lastly, we include more detailed descriptions of the performance metrics used in evaluating each of our models (Section A).

A. PERFORMANCE METRICS
We use the loglikelihood, reduced chi-squared statistic, Bayesian information criterion, and Akaike information criterion to evaluate the goodness of fit for our models. For completeness and pedagogical purposes, we include more detailed descriptions of each of these metrics here.
The loglikelihood (ln(L)) is the natural log of the joint probability of the observations as a function of the parameters of the model. The likelihood is often written as L(θ X) to emphasize that it is a function of the parameters θ given the data X. To find the most optimal θ given X, we sample a range of θ to maximize the likelihood. In practice, we do this by maximizing (ln(L)) for practical purposes. In our case, the loglikelihood is defined as: where x i is the ith observation in x and σ is the estimated variance. The reduced chi-squared statistic (χ 2 v ) is one of the most widely used metrics to assess the goodness of fit. A lower χ 2 v indicates a better fit to the data.The χ 2 v is defined as: where N is the number of observations, k is the number of free parameters, x i is the ith observation in x,x i is the model prediction for the ith observation, and σ is the estimated variance. In general, a χ 2 v > 1 corresponds to a model that does fully capture the complexity of the data and that error variance has been underestimated. In a model where χ 2 v ≪ 1, the model is overfitting and thus either (a) the error variance is overestimated or (b) the model is not properly fitting the noise in the observations. Generally, a χ 2 v close to 1 corresponds to a proper fit and a lower value indicates a better fit to the data.
The Bayesian information criterion (BIC; Schwarz 1978) is based on the likelihood function and adds a penalty term to the number of parameters in the models to reduce overfitting. Models with a lower BIC are generally preferred. The BIC is defined as BIC = k ln(n) − 2 ln(L) whereL is the maximized value of the likelihood function, n are the number of observations, and k are the number of free parameters in the model. The Akaike information criterion (AIC; Akaike 1974) also uses the likelihood function and is founded in information theory where AIC estimates the relative amount of information lost by a given model. The better model is the model where less information is lost. AIC allows one to strike a balance between optimizing the goodness of fit and the simplicity of the model. In this way, AIC allows one to deal with the problem of overfitting and underfitting the data. AIC is defined as AIC = k − 2 ln(L) where the model with the lowest AIC value is the preferred model.