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

Our analysis includes 72 spectra of HAT-P-2 obtained by the CPS Team using HIRES on the W.M. Keck Observatory 10 m 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 yr (2006 September 3–2020 August 31). For comparison, de Wit et al. (2017) used 44 HIRES RVs that span 9 yr (2006 September–2015 October).

In addition to our 72 HIRES RVs, we also included 31 publicly available RVs from the HARPS-N spectrograph at the 3.6 m 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 March 12–2017 August 28. We included the published RVs from Bonomo et al. (2017) and downloaded additional publicly available RVs through the TNG Archive. ⁹

We removed RV measurements that occur during transit (Winn et al. 2007) such that we are not using RVs which are affected by the Rossiter–McLaughlin effect. The RVs are included in the Appendix in Table A1.

The astrometric data used in this work are derived from the Hipparcos-Gaia Catalog of Accelerations (Brandt 2018, 2021). This consists of three independent proper-motion measurements: the Hipparcos proper motion (μ_H), the Gaia Data Release 3 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 have a time baseline of approximately 25 yr (1991–2016).

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).

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 data sets. Our pipeline can notably perform fits of the following stable 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₀ and ω = d ω/dt · time + ω₀. 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, five 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 $\sqrt{e}\sin \omega$ and $\sqrt{e}\cos \omega$—two of which describe a linear ($\dot{\gamma }$) and quadratic ($\ddot{\gamma }$) trend with time, ¹¹ and six which correspond to the offsets (γ) and RV jitter terms for each of the RV data sets. RV jitter describes the instrumental and astrophysical stochastic signals in the data. Astrophysical signals originate from stellar variability such as inhomogeneities 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 parameterization allows mass to vary over time, which is unphysical. Parameterizing 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 five orbital parameters for a secondary companion, which sums to 18 free parameters. ¹² 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 (Vanderburg 2021), ¹³ 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 is <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₂) 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.

We use four different metrics to assess the goodness of fit of our various orbital models: (i) the loglikelihood, (ii) the reduced chi-squared statistic (${\chi }_{r}^{2}$), (iii) the Bayesian information criterion (BIC; Schwarz 1978), (iv) 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.

In order to further investigate the possibility of a long-period companion to HAT-P-2 b, we performed a joint fit of the RV and astrometric data sets. The model used to perform this joint analysis is described in detail in Venner et al. (2021). In summary, we jointly fit the two data sets 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 the ascending node (Ω) for the outer companion; proper motion of system barycentre parameters (μ_bary,RA, μ_bary,Dec); and normalization and jitter parameters for each RV data set. 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

$$\begin{eqnarray}{ \mathcal P }(P,{t}_{d},B)=\left\{\begin{array}{ll}1 & \mathrm{if}\ (P-{t}_{d})\lt B\\ (B+{t}_{d})/P & \mathrm{otherwise}\end{array}\right.,\end{eqnarray} \tag{ 1 }$$

where ${ \mathcal 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(Vanderburg 2021) 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).

For each of our four models we computed the performance metrics described in Section 3.2, and list 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 two-planet model (two Keplerians) has a significantly lower reduced ${\chi }_{r}^{2}=1.035$, lower ΔBIC = 14.90, lower ΔAIC = 12.26, and higher loglikelihood = −432.05 compared to the one -planet model (${\chi }_{r}^{2}=2.806$, ΔBIC = 165.18, ΔAIC = 175.72, loglikelihood = −518.7). Furthermore, the stable two-planet model (Keplerian + quadratic) has a similar reduced ${\chi }_{r}^{2}=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 Figures 1(a)–(f). Overall, the two planet models (either the two Keplerian or Keplerian + quadratic version) are favored in the comparison of these stable orbit models.

Zoom In Zoom Out Reset image size

Table 1. Converged MCMC Results for Each of the Stable Orbital Models and their Goodness-of-fit Values

Model	(i) Stable	(ii) Stable	(iii) Stable
	One Planet	Two Planets	Two Planets
	(Keplerian)	(Two Keplerians)	(Keplerian + Quadratic)
Goodness-of-fit metrics:
Reduced ${\chi }_{r}^{2}$	2.806	1.035	1.045
Loglikelihood	−518.77	−432.05	−432.87
ΔBIC a	165.18	14.90	21.16
ΔAIC a	175.72	12.26	15.90
MCMC values:
Mp (Mj)	${9.233}_{-0.098}^{+0.099}$	${9.033}_{-0.099}^{+0.010}$	${8.998}_{-0.098}^{+0.010}$
P (days)	${5.63346961}_{-0.00000064}^{+0.00000064}$	${5.63346961}_{-0.00000064}^{+0.00000064}$	${5.63346960}_{-0.00000064}^{+0.00000064}$
Tp–2455000 (days)	${289.46980}_{-0.00034}^{+0.00034}$	${289.48642}_{-0.00033}^{+0.00033}$	${289.48642}_{-0.00034}^{+0.00034}$
e0	${0.4990}_{-0.0086}^{+0.0085}$	${0.5003}_{-0.0089}^{+0.0090}$	${0.4960}_{-0.0088}^{+0.0088}$
ω0 (°)	${188.25}_{-0.39}^{+0.39}$	${189.62}_{-0.40}^{+0.40}$	${189.49}_{-0.41}^{+0.41}$
γHARPS−2013−2015	$-{26.8}_{-8.6}^{+8.5}$	${51}_{-17}^{+19}$	${150}_{-16}^{+16}$
γHARPS−2015−2017	$-{20.8}_{-18}^{+18}$	${35}_{-23}^{+24}$	${136}_{-22}^{+22}$
${\gamma }_{\mathrm{HIRES}}$	$-{33.8}_{-4.9}^{+4.9}$	$-{170}_{-14}^{+17}$	${77.4}_{-9.8}^{+9.9}$
Linear ($\dot{\gamma }$)			−0.117−0.011 + 0.011
Quadratic ($\ddot{\gamma }$)			${0.0000194}_{-0.0000023}^{+0.0000023}$
Mp2 (Mj)		${12.1}_{-2.2}^{+2.9}$
P2 (days)		${8721}_{-1141}^{+1407}$
Tp2–2455000 (days)		$-{9092}_{-400}^{+287}$
e02		${0.313}_{-0.069}^{+0.080}$
ω02 (°)		${44}_{-18}^{+17}$

Notes. For each of the models, we specify whether the planets were modeled as a Keplerian or a quadratic trend.

^a The ΔBIC and ΔAIC are computed by subtracting the BIC and AIC, respectively, corresponding to model (iv) two Planets (Keplerian + Quadratic), from the BIC and AIC values. The results for model (iv) are listed in Table 2.

Download table as:  ASCIITypeset image

Table 2. Converged MCMC Results for Each of the Evolving Orbital Models and their Goodness-of-fit Values

Model	(iv) Evolving	(v) Evolving	(vi) Evolving
	One Planet	Two Planets	Two Planets
	(Keplerian)	(Two Keplerians)	(Keplerian + Quadratic)
Goodness-of-fit metrics:
Reduced ${\chi }_{r}^{2}$	2.504	0.963	0.891
Loglikelihood	−502.788	−429.691	−426.918
ΔBIC b	142.47	10.181	0.00
ΔAIC b	147.74	7.546	0.00
MCMC values:
Mp (Mj)	${9.27}_{-0.10}^{+0.10}$	${9.02}_{-0.10}^{+0.10}$	${9.04}_{-0.10}^{+0.10}$
P (days)	${5.63365}_{-0.000040}^{+0.000039}$	${5.63346964}_{-0.00000062}^{+0.00000063}$	${5.633587}_{-0.000043}^{+0.000044}$
Tp (days)	${289.440}_{-0.012}^{+0.011}$	${289.48643}_{-0.00034}^{+0.00033}$	${289.48642}_{-0.00034}^{+0.00034}$
e0	${0.499}_{-0.011}^{+0.010}$	${0.502}_{-0.012}^{+0.011}$	${0.495}_{-0.011}^{+0.011}$
de/dt a (1 yr−1)	${1.75}_{-1.68}^{+1.71}\times {10}^{-3}$	$-{7.3}_{-18.0}^{+19.0}\times {10}^{-4}$	${7.7}_{-18.6}^{+18.6}\times {10}^{-4}$
ω0 (°)	${184.6}_{-1.1}^{+1.1}$	${188.80}_{-0.58}^{+0.57}$	${185.9}_{-1.2}^{+1.2}$
dω/dt a (° yr−1)	${1.12}_{-0.22}^{+0.22}$	${0.228}_{-0.11}^{+0.11}$	${0.76}_{-0.24}^{+0.24}$
γHARPS−2013−2015	$-{23.1}_{-8.7}^{+8.5}$	${45}_{-16}^{+18}$	${142}_{-16}^{+16}$
γHARPS−2015−2017	$-{24}_{-19}^{+19}$	${41.9}_{-22}^{+23}$	${134}_{-23}^{+23}$
${\gamma }_{\mathrm{HIRES}}$	$-{35.9}_{-5.0}^{+5.0}$	$-{20.8}_{-13}^{+15}$	${71}_{-10}^{+10}$
Linear ($\dot{\gamma }$)			$-{0.110}_{-0.011}^{+0.011}$
Quadratic ($\ddot{\gamma }$)			${1.78}_{-0.25}^{+0.24}\times {10}^{-5}$
Mp2 (Mj)		${11.6}_{-2.0}^{+2.8}$
P2 (days)		${8777}_{-1200}^{+1400}$
Tp2 (days)		$-{900}_{-380}^{+300}$
e02		${0.346}_{-0.072}^{+0.078}$
ω02 (°)		${43.0}_{-19.5}^{+16.0}$

Notes. For each of the models, we specify whether the planets were modeled as a Keplerian or a quadratic trend.

^a Only the orbital parameters of HAT-P-2 b are allowed to vary over time. ^b The ΔBIC and ΔAIC are computed by subtracting the BIC and AIC, respectively, corresponding to model (iv) Two Planets (Keplerian + Quadratic), from the BIC and AIC values.

Download table as:  ASCIITypeset image

Next, we can compare various evolving orbit scenarios for HAT-P-2 b. We include the one-planet (Keplerian), two-planet (two Keplerians), and two-planet (Keplerian + quadratic) scenarios, where HAT-P-2 b's e and ω can vary linearly with time. We find that the evolving one-planet model (${\chi }_{r}^{2}=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 Figures 2(a), (b). Next, we find that model (v), the evolving two-planet (two Keplerians) model, provides a slightly better fit (${\chi }_{r}^{2}=0.963$, ΔBIC = 10.181, ΔAIC = 7.546, loglikelihood = −429.691) than model (ii), the stable two-planet (two Keplerians) model. Overall, we find that model (vi), the evolving two-planet model (Keplerian + quadratic) provides the best fit with the lowest ${\chi }_{r}^{2}=0.891$, the lowest BIC = 909.453, the lowest AIC = 877.836, and the highest loglikelihood = −426.918. We plot this model in Figures 2(e), (f). We note that we also performed some additional tests to see how the rate of evolution for e and ω and their significance are impacted by how the long-term companion is modeled, and we include this in the discussion (Section 5.1.4).

Zoom In Zoom Out Reset image size

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 yr extension of the HIRES time series allows us to clearly observe nonlinearity in the RV acceleration. Previously, Knutson et al. (2014) constrained the mass of the outer companion to 8–200 M_j and the semimajor 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, 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⁻¹) 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 yr 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}_{-4.1}^{+7.1}$ yr, comparable to a 25 yr 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 nondetection allows us to exclude orbital periods significantly longer than 15,000 days at 95% confidence. Furthermore, the nondetection constrains the orbital inclination to ≈90° as this is the only orbital configuration consistent with zero net astrometric acceleration over the span of observations.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The key posterior parameters of the outer companion are as follows:

$$\begin{eqnarray}\begin{array}{rcl}{P}_{c} & = & {8500}_{-1500}^{+2600}\ \mathrm{days}\\ {K}_{c} & = & {89}_{-13}^{+39}\ {\rm{m}}\,{{\rm{s}}}^{-1}\\ {e}_{c} & = & {0.37}_{-0.12}^{+0.13}\\ {\omega }_{c} & = & 38^\circ \pm 32^\circ \\ {T}_{p,c} & = & {2454150}_{-800}^{+430}\ \mathrm{BJD}\\ {i}_{c} & = & 90^\circ \pm 16^\circ .\end{array}\end{eqnarray} \tag{ 2 }$$

The longitude of the ascending node (Ω) is unconstrained. With the orbital period and inclination resolved, we find that the mass of this second companion is ${10.7}_{-2.2}^{+5.2}$ M_j. The mass posterior is entirely substellar and corresponds to a substellar-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 substellar-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 and 2. The parameter values for HAT-P-2 c from the joint fit should be considered more robust.

5.1.1. Rates of Orbital Evolution

From our orbital fitting analyses, we find that although a two-planet, ¹⁴ 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-up observations are required to fully confirm or rule out orbital evolution in HAT-P-2 b. For this model, we find that $d\omega /{dt}\,=\,{0.76}_{-0.24}^{+0.24^\circ }$ yr⁻¹ and ${de}/{dt}\,={7.67}_{-18.5}^{+13.2}\times {10}^{-4}$ yr⁻¹, 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° yr⁻¹, de/dt = 8.9 ± 2.8 × 10⁻⁴ yr⁻¹) found in de Wit et al. (2017), they found a slightly more significant 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 Section 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.

Zoom In Zoom Out Reset image size

5.1.2. Comparison to Radial Velocity Analysis in de Wit et al. (2017)

5.1.3. Can the Change in ω Be Explained by General Relativity Alone?

We computed the rate of change expected from general relativity (GR), where

$$\begin{eqnarray}&&{\hat{\omega }}_{\mathrm{GR}}=\displaystyle \frac{3{{GM}}_{* }}{{{ac}}^{2}(1-{e}^{2})}n,\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&=\,1.185\cdot {10}^{-2}\ \mathrm{degrees}\,{\mathrm{yr}}^{-1},\end{eqnarray} \tag{ 4 }$$

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, ${\hat{\omega }}_{\mathrm{GR}}=1.185\cdot {10}^{-2}$ degrees yr⁻¹. However, in our best-fit model, we find $\hat{\omega }=0\buildrel{\circ}\over{.} 76\pm 0\buildrel{\circ}\over{.} 24$ yr⁻¹, 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 planet–star interactions may potentially be able to explain this type of rapid change.

5.1.4. The Impact of Modeling Choices for Planet c on the Rate of Evolution of e and ω

Since previous analyses (Bonomo et al. 2017; de Wit et al. 2017) modeled the long-period companion as linear or quadratic trends, we performed additional tests wherein we modeled 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 among these models (${de}/{dt}={1.75}_{-1.68}^{+1.71}\times {10}^{-3}$ yr⁻¹, $d\omega /{dt}={1.12}{_{-0.22}^{+0.22}}{^\circ}$ yr⁻¹). 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 a 5.1σ trend. The rates are the highest with ${de}/{dt}={3.28}_{-1.75}^{+1.72}\times {10}^{-3}$ yr⁻¹ and $d\omega /{dt}={1.12}{_{-0.22}^{+0.22}}^\circ$ yr⁻¹.

Table 3. Comparison of Evolving Orbit Models Where the Second Companion Is Not Modeled or Modeled as a Linear Trend, a Quadratic Trend, or a Keplerian

Model for HAT-P-2 c	None	Linear	Quadratic	Keplerian
${\chi }_{r}^{2}$	2.503	1.440	0.891	0.963
ΔBIC a	142.47	46.77	0.00	10.181
ΔAIC a	147.74	49.405	0.00	7.5460363
de/dt (1 yr−1)	${1.75}_{-1.68}^{+1.71}\times {10}^{-3}$	${3.28}_{-1.75}^{+1.72}\times {10}^{-3}$	${7.7}_{-18.6}^{+18.6}\times {10}^{-4}$	$-{7.3}_{-18.0}^{+19.0}\times {10}^{-4}$
dω/dt (° yr−1)	${1.12}_{-0.22}^{+0.22}$	${1.12}_{-0.22}^{+0.22}$	${0.76}_{-0.24}^{+0.24}$	${0.228}_{-0.11}^{+0.11}$
σde/dt	1.03	1.90	0.41	0.32
σd ω/dt	5.20	5.11	3.23	2.11

Note.

^a The ΔBIC and ΔAIC are computed by subtracting the BIC and AIC, respectively, corresponding to model (iv) Two Planets (Keplerian + Quadratic) from the BIC and AIC values. Within this table, that is equivalent to subtracting the "Quadratic" model BIC and AIC values.

Download table as:  ASCIITypeset image

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}\times {10}^{-4}$ yr⁻¹, $d\omega /{dt}={0.76}{_{-0.24}^{+0.24}}^\circ$ yr⁻¹). 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\omega /{dt}={0.228}{_{-0.11}^{+0.11}}^\circ$ yr⁻¹) and the direction of de/dt changes (${de}/{dt}=-{7.3}_{-18.0}^{+19.0}\times {10}^{-4}$ yr⁻¹). However, this trend in e is so insignificant, indicating that this 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}_{-1500}^{+2600}\ \mathrm{days}$), and so are many of the other orbital parameters. Thus, further monitoring of the 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.

5.1.5. Spitzer Transits and Secondary Eclipses

5.1.6. Radial Velocity 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 (2023 April 1–2023 April 5). 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 Figures 6(b), (c). In Figure A1, 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 three to four RV measurements on HIRES/Keck I. This could be done anytime in the observing window (2023 April 1–2023 September 28) and would extend the RV baseline from 14 to 17 yr. 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 total RV measurements spanning three orbital periods (i.e., three RVs per orbit). The particular orbital phases are indicated in green in Figures 6(b), (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.

Zoom In Zoom Out Reset image size

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, 2021). Though no acceleration is detected in the Hipparcos-Gaia astrometry, this nondetection is sufficiently informative to constrain the orbit of the second companion. HAT-P-2 c is a ${10.7}_{-2.2}^{+5.2}$ M_j substellar-mass object with an orbital period of ${P}_{c}={8500}_{-1500}^{+2600}$ days, placing it among the longest-period substellar companions that have ever been discovered from RVs.

The constraints of the astrometric nondetection require HAT-P-2 c to have an edge-on orbital inclination (i_c = 90° ± 16°), 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 nondetection. HAT-P-2 c's orbital inclination is consistent with HAT-P-2 b, which 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.

We use the loglikelihood, reduced chi-squared statistic, BIC, and AIC 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 ($\mathrm{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 ($\mathrm{ln}(L)$) for practical purposes. In our case, the loglikelihood is defined as

$$\begin{eqnarray}&&\mathrm{ln}(L)=-\displaystyle \sum _{i}0.5\ast \displaystyle \frac{{\left({x}_{i}-\hat{{x}_{i}}\right)}^{2}}{{\sigma }^{2}}+\mathrm{ln}(\sigma ),\end{eqnarray} \tag{ A1 }$$

where x_i is the ith observation in x and σ is the estimated variance.

The reduced chi-squared statistic (${\chi }_{v}^{2}$) is one of the most widely used metrics to assess the goodness of fit. A lower ${\chi }_{v}^{2}$ indicates a better fit to the data. The ${\chi }_{v}^{2}$ is defined as

$$\begin{eqnarray}&&{\chi }_{v}^{2}=\displaystyle \frac{1}{N-1}\displaystyle \sum _{i}\displaystyle \frac{{\left({x}_{i}-\hat{{x}_{i}}\right)}^{2}}{{\sigma }^{2}},\end{eqnarray} \tag{ A2 }$$

where N is the number of observations, k is the number of free parameters, x_i is the ith observation in x, $\hat{{x}_{i}}$ is the model prediction for the ith observation, and σ is the estimated variance. In general, a ${\chi }_{v}^{2}\gt 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 ${\chi }_{v}^{2}\ll 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 ${\chi }_{v}^{2}$ close to 1 corresponds to a proper fit and a lower value indicates a better fit to the data.

The 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

$$\begin{eqnarray}&&\mathrm{BIC}=k\mathrm{ln}(n)-2\mathrm{ln}(\hat{L}),\end{eqnarray} \tag{ A3 }$$

where $\hat{L}$ is the maximized value of the likelihood function, n is the number of observations, and k is the number of free parameters in the model.

The AIC (Akaike 1974) also uses the likelihood function and is founded in information theory, where it 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 problems of overfitting and underfitting the data. AIC is defined as

$$\begin{eqnarray}&&\mathrm{AIC}=k-2\mathrm{ln}(\hat{L}),\end{eqnarray} \tag{ A4 }$$

where the model with the lowest AIC value is the preferred model.