# The HAT-P-13 Exoplanetary System: Evidence for Spin-Orbit Alignment and a Third Companion

## Abstract

We present new radial-velocity measurements of HAT-P-13, a star with two previously known companions: a transiting giant planet “b” with an orbital period of 3 days, and a more massive object “c” on a 1.2 yr, highly eccentric orbit. For this system, dynamical considerations would lead to constraints on planet b’s interior structure, if it could be shown that the orbits are coplanar and apsidally locked. By modeling the Rossiter-McLaughlin effect, we show that planet b’s orbital angular momentum vector and the stellar spin vector are well-aligned on the sky ( deg). The refined orbital solution favors a slightly eccentric orbit for planet b (), although it is not clear whether it is apsidally locked with c’s orbit ( deg). We find a long-term trend in the star’s radial velocity and interpret it as evidence for an additional body “d”, which may be another planet or a low-mass star. Predictions are given for the next few inferior conjunctions of c, when transits may happen.

## 1 Introduction

Precise radial-velocity 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 () of about 15 Jupiter masses, although its true mass () and orbital inclination () are unknown. In particular, transits of companion c have neither been observed nor ruled out.

Batygin, Bodenheimer, & Laughlin (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 , 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 , 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.

Radial-velocity 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 radial-velocity measurements of HAT-P-13 bearing on all these issues. The new data are presented in § 2. Our analysis is presented § 3, and includes evidence for a third companion “d” (§ 3.1), refined estimates of the eccentricity and apsidal orientation of b’s orbit (§ 3.2), modeling of the RM effect (§ 3.3), and updated predictions for the next inferior conjunction (possible transit window) of companion c (§ 3.4). In § 4 we discuss the implications for further dynamical investigations of HAT-P-13.

## 2 Observations

We observed HAT-P-13 with the High Resolution Spectrograph (HIRES; Vogt et al. 1994) on the Keck I 10m telescope, using the same instrument settings and observing protocols that were used by Bakos et al. (2009) and that 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 Dec 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 radial velocity (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-2. The model curves appearing in those figures are explained in § 3.

## 3 Analysis

Our model for the radial-velocity data took the form

(1) |

where is the calculated RV, and are the radial velocities of non-interacting Keplerian orbits, is the transit-specific “anomalous velocity” due to the Rossiter-McLaughlin effect (§ 3.3), and 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, , allows for a constant radial acceleration, and was included because models with did not fit the data (§ 3.1). We interpret as the acceleration produced by a newly-discovered long-period companion “d”. The third constant, , 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 radial velocity. For this case we used the relation

(2) |

where is the sky-projected stellar rotation speed, is the fractional loss of light, and is the radial velocity of the portion of the stellar photosphere that is hidden by the planet. To calculate 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 or Gaudi & Winn 2007).

Since depends on the planet-to-star radius ratio
, orbital inclination , and impact parameter , all of which are more tightly bounded by observations of
photometric transits than by the RM effect, we simultaneously fitted a
composite -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 and , we fixed , the value obtained by interpolating the Claret (2004) tables,
and allowed to vary freely.^{14}

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 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 with

(3) |

where and are the RV data and associated uncertainties, are the calculated RVs; and are the flux data and associated uncertainties; and are the calculated fluxes.

Each flux uncertainty was taken to be the scatter in the 16 data points contributing to each 4 min time bin. Each RV uncertainty 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 , 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, and , 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.

### 3.1 Evidence for a third companion

The extra acceleration, , was included in the RV model because a model consisting of only two Keplerian orbits gave an unacceptable fit to the data. With , the RV-specific portions of the data and model had and 94 degrees of freedom (. 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 was allowed to vary freely, the best-fitting model had m s yr, and with 93 degrees of freedom. The exact match between and 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 to give an order-of-magnitude constraint

(4) |

where is the orbital distance. By this standard, the newly discovered object could be a 2.5 planet at 5 AU, or a 10 planet at 10 AU, or a 90 (0.09 ) 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, , 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 .
Planet b’s orbit is nearly circular, with . 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 obeys a Rayleigh
distribution (see, e.g., Shen & Turner 2008). We found
(i.e., nonzero with 2.8
significance) and . The
eccentricity of b’s orbit is right on the edge of
detectability.^{15}

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 degrees. The red dashed lines in Figure 5
show the 3 allowed region for . 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.^{16}

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 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 would not be much improved.

### 3.3 The Rossiter-McLaughlin effect

The RV data obtained during transits exhibit a prograde Rossiter-McLaughlin 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 are shown in Figure 6. The finding of 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, km s, is about 1 smaller than the result of km s reported by Bakos et al. (2009).

### 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 time range, at least. The maximum transit duration is 14.9 hr.

## 4 Discussion

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 Rossiter-McLaughlin 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 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 of planet b, as mentioned in § 1. The requirement that the apsidal precession rates of b and c are equal leads to a condition on , because b’s precession rate is significantly affected by its tidal bulge. Subsequent work by Mardling (2010) showed that for a unique determination of it is necessary for the mutual inclination 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 and 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 Rossiter-McLaughlin effect pertains to the angle between planet b’s orbit and the stellar equator, and has no direct bearing on the angle 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 of the star. The critical orbital distance inside which the stellar torque is dominant is (Burns 1986), which is 0.020 AU assuming as for the Sun. This is smaller than the actual orbital distance of 0.043 AU. Hence if were large, then b’s orbit would nodally precess around c’s orbital axis, which would cause periodic variations in . Therefore, at any given moment in the system’s history, we would be unlikely to observe a small value of unless were small. However, it is impossible to draw firm conclusions about 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 .

It may be possible to estimate based on transit timing variations of planet b (Nesvorný & Beaugé 2010; Bakos et al. 2009). An even more direct estimate of could be achieved if transits of c were detected. The existence of transits would show that is nearly , as is . This would suggest 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 Rossiter-McLaughlin 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.

BJD | RV [m s] | Error [m s] |
---|---|---|

Note. – The RV was measured relative to an arbitrary template spectrum; only the differences are significant. The uncertainty given in Column 3 is the internal error only and does not account for “stellar jitter.”

Parameter | Value |
---|---|

Star | |

Mass, [] | |

Radius, [] | |

Projected stellar rotation rate, [km s] | |

Planet b | |

Mass, [] | |

Radius, [] | |

Orbital period, [d] | |

Planet-to-star radius ratio, | |

Star-to-orbit radius ratio, | |

Orbital inclination, [deg] | |

Impact parameter, | |

Time of midtransit, [BJD] | |

Orbital eccentricity, | |

Argument of pericenter, [deg] | |

Velocity semiamplitude, [m s] | |

Projected spin-orbit angle, [deg] | |

Companion c | |

Minimum mass, [] | |

Orbital period, [d] | |

Time of inferior conjunction, [BJD] | |

Orbital eccentricity, | |

Argument of pericenter, [deg] | |

Velocity semiamplitude, [m s] | |

Other system parameters | |

Angle between apsides, [deg] | |

Velocity offset, [m s] | |

[m s yr] |

Year | Month | Date | Hour [UT] | Julian Date | Uncertainty [d] |
---|---|---|---|---|---|

2010 | April | 26 | 7.3 | 2455312.80 | 0.74 |

2011 | July | 16 | 13.9 | 2455759.08 | 0.85 |

2012 | October | 4 | 20.4 | 2456205.35 | 1.00 |

2013 | December | 25 | 3.0 | 2456651.62 | 1.17 |

2015 | March | 16 | 9.6 | 2457097.90 | 1.36 |

### Footnotes

- affiliation: Department of Physics, and Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology, Cambridge, MA 02139
- affiliation: Kavli Institute for Theoretical Physics, UCSB, Santa Barbara, CA 93106, USA
- affiliation: Department of Astrophysics, California Institute of Technology, MC 249-17, Pasadena, CA 91125
- affiliation: Department of Astronomy, University of California, Mail Code 3411, Berkeley, CA 94720
- affiliation: Townes Postdoctoral Fellow, Space Sciences Laboratory, University of California, Berkeley, CA 94720
- affiliation: Department of Astronomy, University of California, Mail Code 3411, Berkeley, CA 94720
- affiliation: Harvard-Smithsonian Center for Astrophysics, 60 Garden St., Cambridge, MA 02138
- affiliation: NSF Fellow
- affiliation: Harvard-Smithsonian Center for Astrophysics, 60 Garden St., Cambridge, MA 02138
- affiliation: Harvard-Smithsonian Center for Astrophysics, 60 Garden St., Cambridge, MA 02138
- affiliation: Department of Physics, and Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology, Cambridge, MA 02139
- affiliation: Kavli Institute for Theoretical Physics, UCSB, Santa Barbara, CA 93106, USA
- affiliation: National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan
- The result, , was consistent with the tabulated value of 0.3068.
- For this reason the results are also sensitive to the choice of priors for the fitting parameters. The results described in this section and given in Table 2 are based on uniform priors for and . If instead uniform priors are adopted for and , then we find .
- If we assume the orbits are apsidally locked, and repeat the fitting procedure with the requirement , we find .

### References

- Bakos, G. Á., et al. 2009, ApJ, 707, 446
- Batygin, K., Bodenheimer, P., & Laughlin, G. 2009, ApJ, 704, L49
- Burns, J. 1986, in Satellites, eds. J. Burns and M. Matthews (University of Arizona Press), p. 117
- Butler, R. P., Marcy, G. W., Williams, E., McCarthy, C., Dosanjh, P., & Vogt, S. S. 1996, PASP, 108, 500
- Butler, R. P., et al. 2002, ApJ, 578, 565
- Claret, A. 2004, A&A, 428, 1001
- Correia, A. C. M., et al. 2010, A&A, 511, A21
- Fabrycky, D. C. 2009, in Transiting Planets, eds. F. Pont, D. Sasselov, M. Holman, Proc. IAU Symposium 253, p. 173-179 (Cambridge Univ. Press)
- Gaudi, B. S., & Winn, J. N. 2007, ApJ, 655, 550
- Howard, A. W., et al. 2009, ApJ, 696, 75
- Léger, A., et al. 2009, A&A, 506, 287
- Mandel, K., & Agol, E. 2002, ApJ, 580, L171
- Mardling, R. A. 2007, MNRAS, 382, 1768
- Mardling, R. A. 2010, MNRAS, in press [arXiv:1001.4079]
- Narita, N., et al. 2010, PASJ, in press [arXiv:1004.2458]
- Nesvorný, D., & Beaugé, C. 2010, ApJ, 709, L44
- Ohta, Y., Taruya, A., & Suto, Y. 2005, ApJ, 622, 1118
- Pál, A. 2008, MNRAS, 390, 281
- Pál, A., et al. 2008, ApJ, 680, 1450
- Queloz, D., et al. 2009, A&A, 506, 303
- Shen, Y., & Turner, E. L. 2008, ApJ, 685, 553
- Vogt, S. S., et al. 1994, Proc. SPIE, 2198, 362
- Winn, J. N., et al. 2005, ApJ, 631, 1215
- Winn, J. N., Holman, M. J., & Fuentes, C. I. 2007, AJ, 133, 11
- Winn, J. N., Johnson, J. A., Albrecht, S., Howard, A. W., Marcy, G. W., Crossfield, I. J., & Holman, M. J. 2009, ApJ, 703, L99
- Wright, J. T. 2005, PASP, 117, 657
- Wright, J. T. 2010, EAS Publications Series, 42, 3
- Wu, Y., & Goldreich, P. 2002, ApJ, 564, 1024