THE HAT-P-13 EXOPLANETARY SYSTEM: EVIDENCE FOR SPIN–ORBIT ALIGNMENT AND A THIRD COMPANION

Precise radial velocity (RV) measurements have revealed more than 30 multiple-planet systems (Wright 2010). However, in only a few cases have transits been detected for any of the planets in those systems. Those cases are potentially valuable because the transit observables—the times of conjunction, orbital inclination, and projected spin–orbit angle, among others—provide a much more complete description of a planetary system, which may in turn give clues about its formation and evolution. The Corot-7 system has two orbiting super-Earths, one of which transits (Léger et al. 2009; Queloz et al. 2009). The HAT-P-7 system has a transiting hot Jupiter in a polar or retrograde orbit, as well as a longer-period companion that could be a planet or a star (Pál et al. 2008; Winn et al. 2009; Narita et al. 2010). The HAT-P-13 system, the subject of this paper, features a G4 dwarf star with two previously known orbiting companions (Bakos et al. 2009). The inner companion (HAT-P-13b, or simply "b" hereafter) is a transiting hot Jupiter in a 2.9 day orbit. The outer companion ("c") has an eccentric 1.2 yr orbit and a minimum mass (M_csin i_c) of about 15 Jupiter masses, although its true mass (M_c) and orbital inclination (i_c) are unknown. In particular, transits of companion c have neither been observed nor ruled out.

Batygin et al. (2009) and Mardling (2010) showed that it may be possible to use the observed state of the system to determine planet b's Love number k₂, a parameter that depends on the planet's interior density distribution. This would be of great interest, as few other methods exist for investigating the interior structure of exoplanets. The method is based on the theoretical expectation that tidal evolution has aligned the apsides of the orbits of b and c. This method has not yet yielded meaningful constraints on k₂, partly because of the large uncertainty in the eccentricity of b's orbit. Another relevant parameter is the mutual inclination between the orbits, which is not known at all.

RV observations are usually powerless to determine mutual inclinations, unless the planets are in a mean-motion resonance (see, e.g., Correia et al. 2010). However, for a transiting planet it is possible to assess the alignment between the orbit and the stellar equator through the Rossiter–McLaughlin (RM) effect. A system with mutually inclined planetary orbits might also be expected to have large angles between the orbits and the stellar equator. In particular, Mardling (2010) presented a formation scenario for HAT-P-13 involving gravitational scattering by a putative third companion, which could have caused large mutual inclinations and a large stellar obliquity.

In this paper, we present new RV measurements of HAT-P-13 bearing on all these issues. The new data are presented in Section 2. Our analysis is presented in Section 3, and includes evidence for a third companion "d" (Section 3.1), refined estimates of the eccentricity and apsidal orientation of b's orbit (Section 3.2), modeling of the RM effect (Section 3.3), and updated predictions for the next inferior conjunction (possible transit window) of companion c (Section 3.4). In Section 4, we discuss the implications for further dynamical investigations of HAT-P-13.

We observed HAT-P-13 with the High Resolution Spectrograph (HIRES; Vogt et al. 1994) on the Keck I 10 m telescope, using the same instrument settings and observing protocols that were used by Bakos et al. (2009) and are used by the California Planet Search (Howard et al. 2009). In particular, we used the iodine gas absorption cell to calibrate the instrumental point-spread function and the wavelength scale. The total number of new spectra is 75, which are added to the 30 spectra presented by Bakos et al. (2009). Of the new spectra, 40 were obtained on the night of 2009 December 27–28, spanning a transit of HAT-P-13b, and were gathered to measure the RM effect. The other 35 were obtained on arbitrary nights. They extend the timespan of the data set by approximately 1 yr, and thereby help to refine the orbital parameters.

The RV of each spectrum was measured with respect to an iodine-free template spectrum, using the algorithm of Butler et al. (1996) with subsequent improvements. All of the spectra obtained by Bakos et al. (2009) were re-reduced, for consistency. Measurement errors were estimated from the scatter in the fits to individual spectral segments spanning 2 Å. The RVs are given in Table 1, and plotted as a function of time in Figures 1 and 2. The model curves appearing in those figures are explained in Section 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Our model for the RV data took the form

$$\begin{equation} V_{\rm calc}(t) = V_{\rm b}(t) + V_{\rm c}(t) + V_{\rm RM}(t) + \gamma + \dot{\gamma }(t-t_0), \end{equation} \tag{ 1 }$$

where V_calc is the calculated RV, V_b and V_c are the RVs of non-interacting Keplerian orbits, V_RM is the transit-specific "anomalous velocity" due to the RM effect (Section 3.3), and $\lbrace \gamma,\dot{\gamma },t_0\rbrace$ are constants. The first constant, γ, specifies the velocity offset between the system barycenter and the arbitrary template spectrum that was used to calculate RVs. The second constant, $\dot{\gamma }$, allows for a constant radial acceleration, and was included because models with $\dot{\gamma }=0$ did not fit the data (Section 3.1). We interpret $\dot{\gamma }$ as the acceleration produced by a newly discovered long-period companion "d." The third constant, t₀, is an arbitrary reference time that was taken to be the time of the first RV datum (BJD 2,454,548.80650).

Our RM model was based on that of Winn et al. (2005), in which simulated spectra are used to calibrate the relation between the phase of the transit and the measured RV. For this case, we used the relation

$$\begin{equation} \Delta V(t) = -(v\sin i_\star) \delta (t) \left[ 0.9833 - 0.0356\left(\frac{v_p(t)}{v\sin i_\star } \right)^2 \right], \end{equation} \tag{ 2 }$$

where vsin i_⋆ is the sky-projected stellar rotation speed, δ is the fractional loss of light, and v_p is the RV of the portion of the stellar photosphere that is hidden by the planet. To calculate v_p, we assumed that the stellar photosphere rotates uniformly with an angle λ between the sky projections of the spin vector and the orbital angular momentum vector (see, e.g., Ohta et al. 2005; Gaudi & Winn 2007).

Since δ depends on the planet-to-star radius ratio R_p/R_⋆, orbital inclination i, and impact parameter b_tra, all of which are more tightly bounded by observations of photometric transits than by the RM effect, we simultaneously fitted a composite i'-band transit light curve based on the photometric data of Bakos et al. (2009). For the photometric model we assumed a quadratic limb-darkening law and used the analytic formula of Mandel & Agol (2002), as implemented by Pál (2008). Because the photometric data are not precise enough to constrain both of the limb-darkening coefficients u₁ and u₂, we fixed u₂ ≡ 0.3251, the value obtained by interpolating the Claret (2004) tables, and allowed u₁ to vary freely.¹⁰ For the RM model, we used a linear law with a fixed coefficient of 0.72, as appropriate for the V band, the approximate spectral range from which the RV signal is drawn.

All together there were 18 adjustable parameters, of which 12 were controlled almost exclusively by the RV data and 6 by the photometric data. The data set had 105 RVs and 107 flux data points. Thus, the total number of degrees of freedom was 194, of which 93 pertained to RVs and 101 to photometry.

We determined the best-fitting parameter values and their 68.3% confidence limits with a Monte Carlo Markov Chain (MCMC) algorithm that we have described elsewhere (see, e.g., Winn et al. 2007). Uniform priors were adopted for all parameters except for b's time of transit and orbital period, for which we adopted Gaussian priors based on the ephemeris of Bakos et al. (2009). We doubled the quoted errors in the ephemeris, out of concern that systematic errors or transit-timing variations have affected the results. The likelihood was taken to be exp(−χ²/2) with

$$\begin{eqnarray} \chi ^2 &=& \sum _{i=1}^{105} \left[ \frac{v_i(\textnormal {obs}) - v_i(\textnormal {calc})}{\sigma _i} \right]^2\nonumber\\ &&+ \sum _{j=1}^{107} \left[ \frac{f_j(\textnormal {obs}) - f_j(\textnormal {calc})}{\sigma _j} \right]^2 = \chi ^2_v + \chi ^2_f, \end{eqnarray} \tag{ 3 }$$

where v_i(obs) and σ_i are the RV data and associated uncertainties, v_i(calc) are the calculated RVs, f_j(obs) and σ_j are the flux data and associated uncertainties, and f_j(calc) are the calculated fluxes.

Each flux uncertainty σ_i was taken to be the scatter in the ≈16 data points contributing to each 4 minute time bin. Each RV uncertainty σ_j was taken to be the quadrature sum of the internally estimated measurement error and a "jitter" term of 3.4 m s⁻¹. The jitter term was set by the requirement χ²_v = 93, the relevant number of degrees of freedom, and is consistent with the empirical jitter estimates of Wright (2005) for stars similar to HAT-P-13. In the best-fitting model, χ²_f = 109.5 and χ²_v = 93.0, the rms photometric residual was 470 ppm and the rms RV residual was 3.6 m s⁻¹.

Table 2 gives the results for the parameter values. Figure 3 shows the RVs as a function of the orbital phases of b and c, expressed in days. In the left panel, the RVs are plotted as a function of the orbital phase of b, after subtracting the calculated contributions to the RV from companions c and d. (The contribution due to d is a linear function of time.) Likewise, the right panel of Figure 3 shows the RVs as a function of the orbital phase of c, after subtracting the calculated contributions from b and d. Figure 4 shows the data over a restricted time range centered on b's transit. The top panel shows the light curve. The bottom panel shows the data after subtracting the calculated orbital RV, thereby isolating the RM effect.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2. Model Parameters for HAT-P-13

Parameter	Value
Star
Mass, M⋆ (M☉)	1.22+0.05−0.10
Radius, R⋆ (R☉)	1.559 ± 0.080
Projected stellar rotation rate, vsin i⋆ (km s−1)	1.66 ± 0.37
Planet b
Mass, Mb (MJup)	0.851 ± 0.038
Radius, Rb (RJup)	1.272 ± 0.065
Orbital period, Pb (days)	2.916250 ± 0.000015
Planet-to-star radius ratio, Rp/R⋆	0.08389 ± 0.00081
Star-to-orbit radius ratio, R⋆/a	0.1697 ± 0.0072
Orbital inclination, i (deg)	83.40 ± 0.68
Impact parameter, btra	0.679 ± 0.043
Time of midtransit, Ttra,b (BJD)	2454779.92976 ± 0.00075
Orbital eccentricity, eb	0.0133 ± 0.0041
Argument of pericenter, ωb (deg)	210+27−36
ebcos ωb	−0.0099 ± 0.0036
ebsin ωb	−0.0060 ± 0.0069
Velocity semiamplitude, Kb (m s−1)	106.04 ± 0.73
Projected spin–orbit angle, λ (deg)	1.9 ± 8.6
Companion c
Minimum mass, Mcsin ic (MJup)	14.28 ± 0.28
Orbital period, Pc (days)	446.27 ± 0.22
Time of inferior conjunction, Tcon,c (BJD)	2455312.80 ± 0.74
Orbital eccentricity, ec	0.6616 ± 0.0054
Argument of pericenter, ωc (deg)	175.29 ± 0.35
eccos ωc	−0.6594 ± 0.0056
ecsin ωc	0.0543 ± 0.0038
Velocity semiamplitude, Kc (m s−1)	440 ± 11
Other system parameters
Angle between apsides, ωb − ωc (deg)	36+27−36
Velocity offset, γ (m s−1)	−100.3 ± 2.0
$\dot{\gamma }$ (m s−1 yr−1)	17.51 ± 0.90

Download table as:  ASCIITypeset image

3.1. Evidence for a Third Companion

The extra acceleration, $\dot{\gamma }$, was included in the RV model because a model consisting of only two Keplerian orbits gave an unacceptable fit to the data. With $\dot{\gamma } = 0$, the RV-specific portions of the data and model had χ²_v = 458.6 and 94 degrees of freedom (χ²_v/N_dof,v = 4.9). The pattern of residuals is displayed in Figure 1. There are large and time-correlated residuals that are not easily attributed to stellar jitter or underestimated measurement errors.

In contrast, when $\dot{\gamma }$ was allowed to vary freely, the best-fitting model had $\dot{\gamma } = 17.5$ m s⁻¹ yr⁻¹, and χ²_v = 93 with 93 degrees of freedom. The exact match between χ²_v and N_dof,v is not significant in itself, as it follows from our choice of 3.4 m s⁻¹ for the jitter term. However, it is significant that an acceptable fit was found for a choice of jitter term that is in line with observations of similar stars. Even more significant is that the correlated pattern of residuals vanished. As shown in Figure 2, the residuals scattered randomly around zero.

The failure of the two-Keplerian model, and the success of a model with an additional radial acceleration, is evidence for a third companion to HAT-P-13 ("d") with a long orbital period. With the limited information available, though, the properties of d are largely unknown. Assuming its orbit to be nearly circular, and its mass to be much smaller than that of the star, we may set $\dot{\gamma } \sim GM_d\sin i_d/a_d^2$ to give an order-of-magnitude constraint

$$\begin{equation} \left(\frac{M_d\sin i_d}{M_{\rm Jup}}\right) \left(\frac{a_d}{{\rm 10 AU}} \right)^{-2} \sim 9.8, \end{equation} \tag{ 4 }$$

where a_d is the orbital distance. By this standard, the newly discovered object could be a 2.5 M_Jup planet at 5 AU, or a 10 M_Jup planet at 10 AU, or a 90 M_Jup (0.09 M_☉) star at 30 AU, etc. The properties of d could be substantially different depending on its eccentricity, argument of pericenter, and time of conjunction. Orbits closer than ∼5 AU would be subject to additional constraints by the requirement of dynamical stability.

More information about d could be gleaned from any significant curvature in the RV signal, beyond the effects of the two Keplerian orbits and a linear trend. We experimented with models that include a "jerk" parameter, $\ddot{\gamma }$, finding this parameter to be highly covariant with the mass, orbital period, and eccentricity of c. More elaborate models and detailed constraints on companion d will only be justified after another few years of observing, when the properties of companion c will have been well established. The uncertainties given in Table 2 must therefore be understood as subject to the assumption that d is producing no significant RV curvature.

3.2. Orbital Eccentricities

Figure 5 shows the results for the orbital eccentricities. Planet c's orbit is strongly eccentric, with e_c = 0.6616 ± 0.0054. Planet b's orbit is nearly circular, with e_b = 0.0133 ± 0.0041. To assess the significance of the detection of eccentricity, it is simpler to examine the components of the eccentricity vector because they obey Gaussian distributions, while e obeys a Rayleigh distribution (see, e.g., Shen & Turner 2008). We found e_bcos ω_b = −0.0099 ± 0.0036 (i.e., nonzero with 2.8σ significance) and e_csin ω_c = −0.0060 ± 0.0069. The eccentricity of b's orbit is right on the edge of detectability.¹¹

Zoom In Zoom Out Reset image size

Because of the low significance of this detection, it is impossible to draw a firm conclusion about whether its orbit is aligned with that of companion c. Our result is Δω ≡ ω_b − ω_c = 36⁺²⁷₋₃₆ deg. The red dashed lines in Figure 5 show the 3σ allowed region for ω_c. The lines intersect the allowed region for planet b, as shown in the right panel. Most of the uncertainty in Δω arises from the poorly determined orientation of b's orbit.¹²

The best way to check on the eccentricity of b's orbit would be to observe an occultation with the Spitzer Space Telescope. The timing of occultations would allow e_bcos ω_b to be determined with an accuracy of about 0.001, several times better than the RV result. However, even after such an observation, considerable uncertainty would remain in Δω because the accuracy in e_bsin ω_b would not be much improved.

3.3. The Rossiter–McLaughlin Effect

The RV data obtained during transits exhibit a prograde RM effect: an anomalous redshift for the first half of the transit, followed by an anomalous blueshift for the second half. The fit to the data is shown in Figure 4, and the resulting constraints on λ and vsin i_⋆ are shown in Figure 6. The finding of λ = 1.9 ± 8.6 deg implies a close alignment between the rotational angular momentum of the star, and the orbital angular momentum of the planet, at least as projected on the sky. Our result for the projected stellar rotation velocity, vsin i_⋆ = 1.66 ± 0.37 km s⁻¹, is about 1σ smaller than the result of 2.9 ± 1.0 km s⁻¹ reported by Bakos et al. (2009).

Zoom In Zoom Out Reset image size

3.4. Inferior Conjunction of Planet C

It is not yet known whether the inclination of c's orbit is close enough to 90° for transits to occur. Observations of transits would reveal the mass and radius of the companion, allow a more precise characterization of its orbit, and place constraints on the mutual inclination of orbits b and c.

Using our results we predicted the times of inferior conjunctions of planet c, which is when transits would occur. The accuracy of the predicted time is limited by correlations with the uncertainties in c's velocity semiamplitude and eccentricity (see Figure 7). Table 3 gives the results. The quoted uncertainties represent 1σ (68.3%) confidence levels. It would be prudent to keep the star under continuous photometric surveillance for the entire 3σ time range, at least. The maximum transit duration is 14.9 hr.

Zoom In Zoom Out Reset image size

Table 3. Predicted Times of Inferior Conjunction for HAT-P-13c

Year	Month	Date	Hour (UT)	Julian Date	Uncertainty (days)
2010	Apr	26	7.3	2,455,312.80	0.74
2011	Jul	16	13.9	2,455,759.08	0.85
2012	Oct	4	20.4	2,456,205.35	1.00
2013	Dec	25	3.0	2,456,651.62	1.17
2015	Mar	16	9.6	2,457,097.90	1.36

Download table as:  ASCIITypeset image

HAT-P-13 was already a noteworthy system, as the first known case of a star with a transiting planet and a second close companion. We have presented evidence for a third companion in the form of a long-term radial acceleration of the star. The properties of the newly discovered long-period companion will remain poorly known until additional RV data are gathered over a significant fraction of its orbital period. Our analysis of the RM effect shows that planet b's orbital axis is aligned with the stellar rotation axis, as projected on the sky. Our new data also agree with the previous finding that the orbit of planet b is slightly eccentric.

The latter two findings are relevant to the second reason why HAT-P-13 is noteworthy: its orbital configuration may represent an example of two-planet tidal evolution. In this scenario, first envisioned by Wu & Goldreich (2002) and investigated further by Mardling (2007), tidal circularization of the inner planet's orbit is delayed due to gravitational interactions with the outer planet. The interactions drive the system into a state of apsidal alignment, where it remains as both orbits are slowly circularized. As it turned out, the specific planetary system that inspired Wu & Goldreich (2002) was irrelevant to their theory, because the "outer planet" was found to be a spurious detection (Butler et al. 2002).

Batygin et al. (2009) welcomed HAT-P-13 as a genuine system that followed the path predicted by Wu & Goldreich (2002), and with the additional virtue that the inner planet is transiting. For this interpretation to be valid, the apsides of b and c must be aligned, whereas we have found the angle between the apsides to be 36⁺²⁷₋₃₆ deg, differing from zero by 1σ. We do not consider this result to be significant enough to draw a firm conclusion, especially in light of the uncertainties due to the ad hoc stellar jitter term and our simplified treatment of the influence of companion d. Further RV monitoring and observations of occultations are needed to make progress.

Batygin et al. (2009) also showed that the existence of transits would allow for an empirical estimate of the tidal Love number k₂ of planet b, as mentioned in Section 1. The requirement that the apsidal precession rates of b and c are equal leads to a condition on k₂, because b's precession rate is significantly affected by its tidal bulge. Subsequent work by Mardling (2010) showed that for a unique determination of k₂ it is necessary for the mutual inclination Δi between orbits b and c to be small. If instead the orbits are mutually inclined, then tidal evolution drives the system into a state in which e_b and Δi undergo oscillations: a cycle in parameter space, instead of a fixed point. Furthermore, Mardling (2010) argued that a large mutual inclination should be considered plausible, or even likely, given c's high eccentricity. She proposed that b and c began with nearly circular and coplanar orbits, but c's orbit was strongly perturbed by an interaction with a hypothetical outer planet. Those same perturbations would likely have tilted c's orbit.

The relation, if any, between the newly discovered HAT-P-13d and Mardling's hypothetical outer planet is unclear. In her scenario, the outer planet is ejected from the system. This seems important to the scenario, as otherwise d would continue interacting with c, and interfere with the tidal evolution of b and c. Thus, unless d's pericenter was somehow raised to avoid further encounters with c, it does not seem likely to have played the role envisioned by Mardling (2010). Of course the scenario could still be correct even if the third companion d was not the scattering agent; a fourth (ejected) companion may have been responsible.

Our study of the RM effect pertains to the angle ψ_⋆,b between planet b's orbit and the stellar equator, and has no direct bearing on the angle Δi between the orbital planes of b and c. However, there is an indirect connection, through the nodal precession that would be caused by mutually inclined orbits. As shown by Mardling (2010), planet b is far enough from the star that its orbital precession rate is likely to be dominated by the torque from c, rather than the quadrupole moment J₂ of the star. The critical orbital distance inside which the stellar torque is dominant is ∼(2J₂a²_cM_c/M_⋆)^0.2 (Burns 1986), which is 0.020 AU assuming J₂ = 2 × 10⁻⁷ as for the Sun. This is smaller than the actual orbital distance of 0.043 AU. Hence if Δi were large, then b's orbit would nodally precess around c's orbital axis, which would cause periodic variations in ψ_⋆,b. Therefore, at any given moment in the system's history, we would be unlikely to observe a small value of ψ_⋆,b unless Δi were small. However, it is impossible to draw firm conclusions about Δi because of the dependence on initial conditions, the possible effects of companion d, and the fact that only the sky-projected angle λ is measured rather than the true obliquity ψ_⋆,b.

It may be possible to estimate Δi based on transit-timing variations of planet b (Nesvorný & Beaugé 2010; Bakos et al. 2009). An even more direct estimate of Δi could be achieved if transits of c were detected. The existence of transits would show that i_c is nearly 90°, as is i_b. This would suggest Δi is small, although it would still be possible that the orbits are misaligned and their line of nodes happens to lie along the line of sight. The most definitive result would be obtained by observing the RM effect during transits of c, and comparing c's value of λ with that of planet b. In effect, the rotation axis of the star would be used as a reference line from which the orientation of each orbit is measured (Fabrycky 2009). This gives additional motivation to observe HAT-P-13 throughout the upcoming conjunctions of companion c.

We thank Dan Fabrycky for helpful conversations, especially about the dynamical implications of our results. We also thank Debra Fischer and John Brewer for investigating the spectroscopic determination of the stellar rotation rate. We are grateful to Scott Gaudi and Greg Laughlin for comments on the manuscript.

J.N.W. gratefully acknowledges support from the NASA Origins program through award NNX09AD36G and the MIT Class of 1942. A.W.H. acknowledges a Townes Postdoctoral Fellowship from the Space Sciences Laboratory at UC Berkeley. G.A.B. was supported by NASA grant NNX08AF23G and an NSF Astronomy & Astrophysics Postdoctoral Fellowship (AST-0702843). G.T. acknowledges partial support from NASA grant NNX09AF59G. S.A. acknowledges the support of the Netherlands Organisation for Scientific Research (NWO). N.N. was supported by a Japan Society for Promotion of Science (JSPS) Fellowship for Research (PD: 20-8141). J.N.W. and N.N. were also supported in part by the National Science Foundation under Grant No. NSF PHY05-51164 (KITP program "The Theory and Observation of Exoplanets" at UCSB).

The data presented herein were obtained at the W. M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California, and the National Aeronautics and Space Administration, and was made possible by the generous financial support of the W. M. Keck Foundation. We extend special thanks to those of Hawaiian ancestry on whose sacred mountain of Mauna Kea we are privileged to be guests. Without their generous hospitality, the Keck observations presented herein would not have been possible.

Facilities: Keck:I (HIRES) - KECK I Telescope