HAT-P-31b,c: A TRANSITING, ECCENTRIC, HOT JUPITER AND A LONG-PERIOD, MASSIVE THIRD BODY^*

In previous HATNet papers, we have used a simplified model for the transit light curve of the HATNet data. For HAT-P-31, no precise photometry exists and thus we fit the HATNet data using a more sophisticated quadratic limb-darkening Mandel & Agol (2002) algorithm with limb-darkening coefficients interpolated from the tables by Claret (2004). One caveat is that the instrumental blending factor, B_inst, is unknown as discussed earlier in Section 2. We point out that experience with previous HATNet planets suggests B_inst is within 2σ of unity for light curves processed using reconstructive TFA in all cases and thus can be accounted for by conservatively doubling the uncertainties on p and R_P. Further support for a B_inst factor not greatly different from unity come from the fact HAT-P-31 is fairly isolated and there are no neighbors in 2MASS or a DSS image which contribute significant flux to the HATNet aperture.

Due to the low-precision photometry, the stellar density cannot be determined to high precision using the method of Seager & Mallén-Ornelas (2003) and in fact spectroscopic estimates were found to be more precise. However, we can reverse this well-known trick by implementing a Bayesian prior in our fitting process for the stellar density.

We use the spectroscopically determined stellar density from Section 3.1 as a prior in our fits. Since the period of the transiting planet is well constrained for even low signal-to-noise photometry, reasonably precise estimates for P_b and ρ_* are possible. With these two constrained, (a_b/R_*) is therefore also constrained. The transit light curve is essentially characterized by four parameters, τ_b, δ_b, b_b, and (a_b/R_*) and thus one of these parameters is constrained by the combination of the transit times and the stellar density prior alone.

The photometry which is fitted in the global modeling is corrected for instrumental systematics through the EPD and reconstructive TFA correction procedures prior to performing the fit (see Section 2.1 and Bakos et al. 2010 for details).

For the RV fits, we found a single planet fit gave a very poor fit to the observations (χ² = 204.2 for 40 RV points) and quickly appreciated some kind of second signal was present. Exploring different models, such as Trojan offsets, polynomial time trends, and outer companions, (see Table 6 for comparison), we found that an eccentric transiting planet with a quadratic trend in the RVs was the preferred model. The pivot point (t_pivot) of the polynomial models, including the drift and quadratic trends, was selected to be the weighted mean of the RV time stamps.

The most likely physical explanation for a quadratic trend is a third body in the system, described by a Keplerian model. Indeed, the Keplerian model provides an improved χ² for a circular orbit and then again for an eccentric orbit but the extra degrees of freedom penalize our model selection criterion. We also found that these models were highly unconstrained and convergence in the associated fits was unsatisfactory. An illustration of the lack of convergence is shown in Figure 5. Therefore, we will adopt the quadratic model in our final reported parameters in Table 7.

Zoom In Zoom Out Reset image size

Table 7. Global Fit Results for HAT-P-31

Parametera	Value
Fitted parameters
Pb (days)	5.005425+0.000091−0.000092
τb (BJDTDB − 2,450,000)	4320.8866+0.0051−0.0053
$[\tilde{T}_1]_b$ (s)	16300+1000−1000
p2b (%)	0.65+0.18−0.12
bb	0.57+0.23−0.31
OOT	1.00061+0.00016−0.00016
Kb (ms−1)	232.5+1.1−1.1
ebsin ωb	−0.2442+0.0043−0.0043
ebcos ωb	0.0185+0.0080−0.0079
γKeck (ms−1)	−29.0+1.4−1.4
γSubaru (ms−1)	18.8+3.1−3.1
γFIES (ms−1)	92.7+5.6−5.6
$\dot{\gamma }$ (ms−1 day−1)	0.141+0.025−0.025
$\ddot{\gamma }$ (ms−1 day−2)	0.00226+0.00021−0.00021
SME derived quantities
u1b,e	0.2078*
u2b$^,$e	0.3550*
ρ*c$^,$e (g cm−3)	0.69+0.34−0.26
HAT-P-31b derived properties
Ψb	0.4737+0.0065−0.0064
eb	0.2450+0.0045−0.0045
ωb (°)	274.3+1.8−1.8
log (gb (cgs))	3.61+0.15−0.32
(ab/R*)	8.9+1.4−2.3
ib (°)	87.1+1.8−2.7
[T1, 4]b (s)	18500+1700−1200
[T2, 3]b (s)	14200+1300−2400
pb	0.080+0.022−0.015
Mb (MJ)	2.171+0.105−0.077
Rb (RJ)	1.07+0.24−0.16
Corr(Rb,Mb)	0.795
ρb (g cm−3)	2.18+1.24−0.93
ab (AU)	0.055+0.015−0.015
[Teq]b (K)	1450+230−110
Θe	0.190+0.036−0.056
Fperi (109 erg s−1 cm−2)	1.69+1.39−0.44
Fap (109 erg s−1 cm−2)	0.62+0.51−0.16
〈F〉f (109 erg s−1 cm−2)	0.99+0.82−0.26
HAT-P-31c derived quantities
τc + 0.25Pc (BJDTDB − 2,450,000)	5254.9+7.4−6.8
KcP−2c (ms−1 day−2)	7.64+0.71−0.70
Pc (years)	⩾2.8
Kc (ms−1)	⩾60
Mc (MJ)	⩾3.4

Notes. ^a Quoted values are medians of MCMC trials with errors given by 1σ quantiles. "b" subscripts refer to planet HAT-P-31b and "c" subscripts refer to HAT-P-31c. ^bFixed parameter. ^cParameter is treated as a prior. ^dParameter determined using SME and YY isochrones. ^eThe Safronov number is given by $\Theta = \frac{1}{2}(V_{\rm esc}/V_{\rm orb})^2 = (a/R_{p})(M_{p}/ M_\star)$ (see Hansen & Barman 2007). ^fIncoming flux per unit surface area, averaged over the orbit.

Download table as:  ASCIITypeset image

The quadratic model may be used to infer some physical parameters of the third planet. To make some meaningful progress, we will assume the outer planet is on a circular orbit. One may compare the quadratic and Keplerian model descriptions via

$$\begin{eqnarray} \mathrm{RV}_{\mathrm{quad}}^c &= \gamma ^{\prime } + \dot{\gamma } (t - t_{\mathrm{pivot}}) + 0.5 \ddot{\gamma } (t - t_{\mathrm{pivot}})^2 \end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray} \mathrm{RV}_{\mathrm{Kep}}^c &= \gamma ^{\prime } - K_c \sin \Bigg (\frac{2\pi (t-\tau _c)}{P_c}\Bigg). \end{eqnarray} \tag{ 2 }$$

By differentiating both expressions and solving for the time when the signals are minimized, one may write

$$\begin{eqnarray} \tau _c + \frac{P_c}{4} &= \frac{ -\dot{\gamma } + \ddot{\gamma }t_{\mathrm{pivot}} }{\ddot{\gamma }}. \end{eqnarray} \tag{ 3 }$$

Differentiating both RV models with respect to time twice and evaluating at the moment when both signals are minimized yields

$$\begin{eqnarray} \frac{K_c}{P_c^2} &= \frac{\ddot{\gamma }}{4\pi ^2}. \end{eqnarray} \tag{ 4 }$$

Equations (3) and (4) may therefore be used to determine some information about HAT-P-31c.

To evaluate the statistical significance of HAT-P-31c, we performed an F-test between the one-planet and two-planet models. In both cases, HAT-P-31b is assumed to maintain non-zero orbital eccentricity. Assuming HAT-P-31c is on a circular orbit, the false-alarm probability (FAP) from an F-test is 3.0 × 10⁻¹², or 7.0σ. Assuming HAT-P-31c is on an eccentric orbit requires two more degrees of freedom and thus reduces the FAP to 1.3 × 10⁻¹⁰, or 6.4σ. From a statistical perspective then, the presence of HAT-P-31c is highly secure. We point out this determination of course assumes purely Gaussian uncertainties and no outlier measurements.

We utilize a Metropolis–Hastings Markov Chain Monte Carlo (MCMC) algorithm to globally fit the data, including the stellar density prior (our routine is described in Kipping & Bakos (2011)). To ensure the parameter space is fully explored, we used five independent MCMC fits which stop once 1.25 × 10⁵ trials have been accepted and burn out the first 20%. This leaves us with a total of 5 × 10⁵ points for the posterior distributions. At the end of the fit, a more aggressive downhill simplex χ² minimization is used, for which the final solution is used for Figures 2 and 3.

There were 14 free parameters in the global fit: {τ_b, p²_b, $[\tilde{T}_1]_b$, b_b, P_b, e_bsin ω_b, e_bcos ω_b, γ_Keck, γ_FIES, γ_Subaru, K_b, $\dot{\gamma }$, $\ddot{\gamma }$, OOT}, which we elaborate on here. τ_b is the time of transit minimum (Kipping 2011), frequently dubbed by the misnomer "mid-transit time." p²_b is the ratio-of-radii squared and P_b is the orbital period. $[\tilde{T}_1]_b$ is the "one-term" approximate equation (Kipping 2010) for the transit duration between the instant when the center of the planet crosses the stellar limb to exiting under the same condition. b_b is the impact parameter, defined as the sky-projected planet–star separation in units of the stellar radius at the instant of inferior conjunction. e_b is the orbital eccentricity and ω_b is the associated position of pericenter. γ terms relate to the instrumental offsets for the RVs. Similarly, OOT is the out-of-transit flux for the HATNet photometry. Finally, K_b is the RV semi-amplitude, and $\dot{\gamma }$ (drift) and $\ddot{\gamma }$ (curl) are the first and second time derivatives of γ.

Final quoted results are the median of the marginalized posterior for each fitted parameter with 34.15% quantiles either side for the 1σ uncertainties (see Table 7). The uncertainties on p, p², and R_P have been conservatively doubled for reasons described in Section 2. Histograms of the posterior distributions for the fitted parameters are provided in Figure 6. We find that the stellar jitter of this star is at or below the measurement uncertainties (∼2 m s⁻¹).

Zoom In Zoom Out Reset image size

Table 7 provides estimates for some minimum limits on various parameters of interest relating to HAT-P-31c. These limits are determined by the known constraints on the minimum P_c. We determined this value by forcing a circular orbit Keplerian fit for planet c through the data, stepping through a range of periods from 1 year to 5 years in 1 day steps. The minimum limit on P_c is defined as when Δχ² = 1, relative to the quadratic trend fit, occurring at 2.8 years.

Here we discuss our procedure to test the dynamical stability of two possible orbital configurations. It should be noted that the eccentricity of HAT-P-31b is highly secure but the eccentricity of HAT-P-31c remains unclear. For this reason, we repeat our simulations assuming both a circular and eccentric orbit for HAT-P-31c, beginning with the former. We utilize the Systemic Console (Meschiari et al. 2009) for this purpose assuming a coplanar configuration. Employing the Gragg–Bulirsch–Stoer integrator, orbital evolution was computed for 250,000 years for the HAT-P-31 system (see Figure 7).

Zoom In Zoom Out Reset image size

We first consider the circular case. The orbital period and mass of HAT-P-31c are non-convergent parameters and so we can only provide an orbital solution which gives a good fit to the data, but is not necessarily unique. To this end, we proceeded to input HAT-P-31c with P_c = 4.86 years and M_c = 13.0 M_J, corresponding to the solution presented in Table 6. This test revealed minor evolution over the course of our simulations, indicating a stable and essentially static configuration.

To test the eccentric fit, we again used the lowest χ² solution presented earlier in Table 6, corresponding to P_c = 4.82 years and e_c = 0.285. We found that the system was also stable over 250,000 years (see Figure 7). However, the simulations do show the eccentricity of planet b varying sinusoidally over a timescale of ∼125,000 years with an amplitude of ∼0.01. The eccentricity of planet c also varies in anti-phase but with a much smaller amplitude. The eccentric evolution shows faster apsidal precession for planet b, but this is unlikely to be observable through changes in the transit duration. We estimate the duration will change by 0.2 s over 10 years (corresponding to Δω_b = 0027) using the expressions of Kipping (2010).

The orbital period and semimajor axes of both bodies were stable over the 250,000 years of integration considered here.

We tried adding a habitable-zone Earth-mass planet on a circular orbit into the system and testing stability. One may argue that the probable history of this system involved the inward migration of HAT-P-31b and that this migration through the inner protoplanetary disk would essentially eliminate the possibility of an Earth-like planet forming in the habitable-zone. However, Fogg & Nelson (2007) have shown that this not necessarily true. In their simulations, it is found that >60% of the solid disk survives, including planetesimals and protoplanets, by being scattered by the giant planet into external orbits where dynamical friction is strong and terrestrial planet formation is able to resume. In one simulation, a planet of 2 M_⊕ formed in the habitable-zone after a hot-Jupiter passed through and its orbit stabilized at 0.1 AU.

For a planet to receive the same insolation as the Earth, we estimate P = 604 days. For our circular orbit solution of HAT-P-31c, the habitable-zone Earth-mass planet is stable for over 100,000 years. For our eccentric orbit solution, the Earth-like planet is summarily ejected in less than 1000 years.