HAT-P-16b: A 4 M_J PLANET TRANSITING A BRIGHT STAR ON AN ECCENTRIC ORBIT^*,^†

The transits of HAT-P-16b were detected with the HAT-6 and HAT-7 telescopes in Arizona and the HAT-8 and HAT-9 telescopes in Hawaii. The star GSC 2792-01700 (also known as 2MASS 00381756+4227470; α = 00^h38^m1756, δ = +42°27'471; J2000; V = 10.812; Droege et al. 2006) lies in the intersection of three separate HATNet fields internally labeled as 123, 163, and 164. Field 123 was observed with an I filter and a 2K×2K CCD, while fields 163 and 164 were observed through R filters with 4K×4K CCDs. All three fields were observed with a 5 minute exposure time and at a 5.5 minute cadence in the period between 2007 July and 2008 September during which over 12,000 exposures were gathered for the three fields. Each image contains between 27,000 and 76,000 stars down to a magnitude of I ∼ 13 for field 123 and R ∼ 15 for fields 163 and 164, yielding a photometric precision for the brightest stars in the field of ∼2 mmag for field 123 and ∼5 mmag for fields 163 and 164.

Frame calibration, astrometry, and aperture photometry were done in an identical way to recent HATNet TEP discoveries, as described in Bakos et al. (2010) and Pál (2009). The resulting light curves were decorrelated (cleaned of trends) using the External Parameter Decorrelation (EPD) technique in "constant" mode (see Bakos et al. 2010) and the Trend Filtering Algorithm (TFA; see Kovács et al. 2005). The light curves were searched for periodic box-like signals using the Box Least-Squares method (BLS; see Kovács et al. 2002). A significant signal was detected in the light curve of GSC 2792-01700, with a depth of ∼10.2 mmag, and a period of P = 2.7760 days. The dip had a relative duration (first to last contact) of q ≈ 0.0460 ± 0.0005, corresponding to a total duration of Pq ≈ 3.062 ± 0.031 hr (see Figure 1).

Zoom In Zoom Out Reset image size

Transiting planet candidates found by ground-based, wide-field photometric surveys must undergo a rigorous vetting process to eliminate the many astrophysical systems mimicking transiting planets (called false positives), the rate of which has proved to be much higher than the occurrence of true planets (10–20 times higher). Low signal-to-noise ratio (S/N) high-resolution reconnaissance spectra are used to extract stellar parameters such as effective temperature, gravity, metallicity, and rotational and radial velocities (RVs) to rule out these false positives. Examples of false positives which are discarded are eclipsing binaries and triple systems. The latter can be either hierarchical or chance alignment systems where the light of the eclipsing pair of stars is diluted by the light of a third brighter star. Rapidly rotating and/or hot host stars whose spectrum is unsuitable for high-precision velocity work are also discarded.

We acquired seven reconnaissance spectra with the CfA Digital Speedometer (DS; Latham 1992) mounted on the FLWO 1.5 m Tillinghast Reflector between 2008 December and 2009 January. The extracted modest precision RVs gave a mean absolute RV = −16.83 km s⁻¹ with an rms of 0.51 km s⁻¹, which is consistent with no detectable RV variation. The stellar parameters derived from the spectra (Torres et al. 2002), including the effective temperature T_eff⋆ = 6000 ± 100 K, surface gravity log g_⋆ = 4.0 ± 0.25 (log cgs) and projected rotational velocity vsin i = 3.8 ± 1.0 km s⁻¹, correspond to an F8 dwarf.

The high significance of the transit detection by HATNet, together with the stellar spectral type and small RV variations, encouraged us to gather high-resolution, high S/N spectra to determine the orbit of the system. We have taken 21 spectra between 2009 August and October using the FIbre-fed Échelle Spectrograph (FIES) at the 2.5 m Nordic Optical Telescope (NOT) at La Palma, Spain (Djupvik & Andersen 2009). We used the medium- and the high-resolution fibers (13 projected diameter) with resolving powers of λ/Δλ ≈ 46,000 and 67,000, respectively, giving a wavelength coverage of ${\sim} 3600\hbox{--}7400\ \rm {\AA}$. We used the wavelength range from approximately ${\sim} 4000\hbox{--}5500\ \rm {\AA}$ to determine the RVs. The exposure time was approximately 20 minutes yielding an S/N from 45 to 85 pixel⁻¹ in the wavelength range used.

We also used the High-Resolution Echelle Spectrometer (HIRES) instrument (Vogt et al. 1994) on the Keck I telescope located on Mauna Kea, Hawaii. We obtained six exposures with an iodine cell, plus one iodine-free template with Keck I/HIRES. The S/N per pixel ranged from 90 to 120 and the exposure time was approximately 10 minutes. The observations were made on six nights between 2009 July 3 and 2009 October 29. The width of the spectrometer slit used on HIRES was 086, resulting in a resolving power of λ/Δλ ≈ 55,000, with a wavelength coverage of ∼3800–8000 Å. The iodine gas absorption cell was used to superimpose a dense forest of I₂ lines on the stellar spectrum and establish an accurate wavelength fiducial (see Marcy & Butler 1992) and only the part of the spectrum where the iodine lines are present was used in the determination of the RVs. Relative RVs in the solar system barycentric frame were derived as described by Butler et al. (1996), incorporating full modeling of the spatial and temporal variations of the instrumental profile.

The final RV data and their errors (for both instruments) are listed in Table 1. The folded data are plotted in Figure 2. The systemic gamma velocities (that are the result of the global modeling, as laid out in Section 3) have been subtracted to ensure a common zero point. The best orbital fit (see Section 3) is superimposed in the figure.

Zoom In Zoom Out Reset image size

In the same figure, we also show the relative S index, which is a measure of the chromospheric activity of the star derived from the flux in the cores of the Ca ii H and K lines. This index was computed following the prescription given by Vaughan et al. (1978), after matching each spectrum to a reference spectrum using a transformation that includes a wavelength shift and a flux scaling that is a polynomial as a function of wavelength. The transformation was determined on the regions of the spectra that are not used in computing this indicator. We also measured the R'_HK index for the system to be R'_HK = −4.862, as described by Noyes et al. (1984). Note that our relative S index has not been calibrated to the scale of Vaughan et al. (1978). We do not detect any significant variation of the index correlated with orbital phase; such a correlation might have indicated that the RV variations could be due to stellar activity, casting doubt on the planetary nature of the candidate.

Since this is the first time the NOT has been used to determine the orbit of a HATNet planet, we describe briefly the extraction of the RVs. We have built a custom pipeline to rectify and cross-correlate spectra from Échelle spectrographs with the goal of providing very precise RVs. The first step in the extraction process is to remove the bias level and crop the raw images. To effectively remove cosmic rays, the observation is divided into three separate exposures, which enables us to combine the raw images using median filtering, removing virtually all cosmic rays. We use a flat field to trace the Échelle orders and we rectify the spectra by using the "optimal extraction" algorithm (Hewett et al. 1985; Horne 1986). The blaze function is determined by fitting cubic splines to a combined high S/N flat field exposure and is saved separately in order to preserve the original flux in the stellar exposure. By dividing the normalized blaze function into the rectified flat field spectrum, we can determine the pixel to pixel variations and correct for these. The scattered light in the two-dimensional raw image is determined and removed by masking out the signal in the Échelle orders and fitting the inter-order scattered light flux with a two-dimensional polynomial.

Thorium argon (ThAr) calibration images are obtained through the science fiber before and after the stellar observation. The two calibration images are combined to form the basis for the fiducial wavelength calibration. We have determined that the best wavelength solution is achieved by choosing an exposure time that saturates the strong argon lines but enhances the forest of weaker thorium lines.

FIES is not a vacuum spectrograph and is only temperature controlled to 0.1°C. Consequently, the RV errors are dominated by our limited ability to correct for shifts due to pressure, humidity, and temperature variations. In order to successfully remove these large variations (>1.5 km s⁻¹), it is critical that the ThAr light travels through the same light path as the stellar light and thus acts as an effective proxy to remove these variations. We have therefore chosen to interleave the stellar observation between two ThAr exposures instead of using the simultaneous ThAr technique, which may not exactly describe the induced shifts in the science fiber. The centers of the ThAr lines are found by fitting a Gaussian function to the line profiles and a two-dimensional fifth-order Legendre polynomial is used to describe the wavelength solution.

Once the spectra have been extracted, a multi-order cross-correlation is performed order by order. First, the spectra are interpolated to a common oversampled log wavelength scale with the same monotonic wavelength increment in all orders. A high and low pass filter is applied in the Fourier domain and the ends of the spectra are apodized with a cosine bell function. The orders are cross-correlated using a fast Fourier transform and the cross-correlation function (CCF) for each order is co-added. This automatically weights each order by its flux, giving more weight to the orders with more photons. This CCF peak is fitted with a Gaussian function to determine its center. The internal precision is estimated by finding the RV for each order and the RV error is thus $\sigma =\mathrm{\mbox{rms}}(v)/\sqrt{N}$, where v is the RV of the individual orders and N is the number of orders. For the final RVs, a template spectrum is constructed by shifting and co-adding all the observed spectra, and the individual spectra are cross-correlated against this co-added template spectrum to minimize the contribution of noise in the template. We note that FIES has also been used to obtain a spectroscopic orbit for one of the first Kepler planets, namely Kepler-7b (Latham et al. 2010).

The observations were gathered while the moon was up, but the moon was well separated from the target on the sky and we had mostly clear conditions. Furthermore, the relative RV of the moon was well separated from the target, so if any lunar contamination did occur, it would have a negligible effect on the RVs.

To confirm the transit signal and obtain high-precision light curves for modeling the system, we conducted photometric follow-up observations with the KeplerCam CCD on the FLWO 1.2 m telescope. We observed four transit events of HAT-P-16 on the nights of 2009 September 10 and 21, October 19 and 30 (Figure 3). On the four nights, 303, 181, 293, and 471 frames were acquired with a cadence of 41, 52, 41, and 32 s (30, 40, 30, and 20 s of exposure time) in the Sloan i band, respectively.

Zoom In Zoom Out Reset image size

The reduction of the images, including frame calibration, astrometry, and photometry, was performed as described in Bakos et al. (2010). We also performed EPD and TFA to remove trends simultaneously with the light curve modeling (for more details, see Section 3). The final light curves are shown in the upper plots of Figure 3, superimposed with our best-fit transit light curve models (see also Section 3); the photometry is provided in Table 2.

The stellar atmospheric parameters were derived from the template spectrum obtained with the Keck/HIRES instrument. We used the Spectroscopy Made Easy (SME) analysis package of Valenti & Piskunov (1996), along with the atomic-line database of Valenti & Fischer (2005). This yielded the following initial values and uncertainties (which we have conservatively increased to include our estimates of the systematic errors): effective temperature T_eff⋆ = 6175 ± 80 K, stellar surface gravity log g_⋆ = 4.44 ± 0.10 (cgs), metallicity [Fe/H] = 0.15 ± 0.06 dex, and projected rotational velocity vsin i = 4.4 ± 0.5 km s⁻¹.

We could now use the effective temperature and the surface gravity found by SME to determine other stellar parameters, such as the mass, M_⋆, and radius, R_⋆, using model isochrones. However, there may be strong correlations between effective temperature, gravity, and metallicity in the spectroscopic determination of these parameters. Also, the effect of log g_⋆ on the spectral line shapes is typically subtle, and as a result log g_⋆ is generally a rather poor luminosity indicator. Instead, we used the a/R_⋆, the normalized semimajor axis (related to ρ_⋆, the mean stellar density), which can be derived directly from the light curves (Sozzetti et al. 2007).

The stellar parameters were determined simultaneously with the modeling of the light curves and RVs, as described next. We began by adopting the values of T_eff⋆, [Fe/H], and log g_⋆ from the SME analysis to fix the quadratic limb-darkening coefficients from the tabulations by Claret (2004), which are needed to model the light curves (a_i, b_i). This modeling yields the probability distribution of a/R_⋆ via a Monte Carlo approach, which is fully described in Section 3.3. We then used the a/R_⋆ distribution together with Gaussian distributions for T_eff⋆ and [Fe/H], with 1σ uncertainties as reported previously, to estimate M_⋆ and R_⋆ by comparison with the stellar evolution models of Yi et al. (2001). This was done for each of the 15,000 simulations. Parameter combinations corresponding to unphysical locations in the H-R diagram (4% of the trials) were ignored and replaced with another randomly drawn parameter set. The resulting stellar parameters and their uncertainties were determined from the a posteriori distributions obtained in this way.

In particular, the resulting surface gravity of log g_⋆ = 4.34 ± 0.03 is somewhat different from that derived in the initial SME analysis, which is not surprising in view of the possible strong correlations among T_eff⋆, [Fe/H], and the surface gravity. Therefore, in a second iteration, we adopted this value of log g_⋆ and held it fixed in a new SME analysis, adjusting only T_eff⋆, [Fe/H], and vsin i. This gave T_eff⋆ = 6158 ± 80 K, $\rm [Fe/H]=+0.17\pm 0.08$, and vsin i = 3.5 ± 0.5 km s⁻¹, which we adopt as the final atmospheric values for the star. Finally, we repeated the calculation of the stellar mass and radius with these values and the stellar evolution models and obtained M_⋆ = 1.218 ± 0.039 M_☉, R_⋆ = 1.237 ± 0.054 R_☉, and L_⋆ = 1.97 ± 0.22 L_☉, along with other parameters summarized in Table 3.

The model isochrones from Yi et al. (2001) for metallicity [Fe/H] = +0.17 are plotted in Figure 4, with the final choice of effective temperature T_eff⋆ and a/R_⋆ marked, and encircled by the 1σ and 2σ confidence ellipsoids.

Zoom In Zoom Out Reset image size

The stellar evolution modeling provides color indices that may be compared against the measured values as a sanity check. The best available measurements are the near-infrared magnitudes from the Two Micron All Sky Survey (2MASS) catalog (Skrutskie et al. 2006), J_2MASS = 9.850 ± 0.022, H_2MASS = 9.623 ± 0.022, and K_2MASS = 9.553 ± 0.020, which we have converted to the photometric system of the models (ESO) using the transformations by Carpenter (2001). The resulting measured color index is J − K = 0.319 ± 0.032. This is within 1σ of the predicted value from the isochrones of J − K = 0.33 ± 0.02. The distance to the object may be computed from the absolute K magnitude from the models (M_K = 2.74 ± 0.10) and the 2MASS K_s magnitude, which has the advantage of being less affected by extinction than optical magnitudes. The result is 235 ± 10 pc, where the uncertainty excludes possible systematics in the model isochrones that are difficult to quantify.

We performed a line bisector analysis of the observed spectra following Torres et al. (2007) to explore the possibility that the main cause of the RV variations is actually distortions in the spectral line profiles due to stellar activity or a nearby unresolved eclipsing binary. The bisector analysis was performed on the Keck spectra as described in Section 5 of Bakos et al. (2007a); the spectra from the NOT were analyzed in a similar manner.

The resulting BS variations are plotted in Figure 2 and have a low amplitude compared to the orbital semi-amplitude and show no correlation with the RVs. We have set the mean of the bisector variations from the three sets of observations to zero for a better comparison between the different instruments. We conclude that there is no significant bisector variation and therefore modeling the RV variation as a result of the reflex motion of the host star of a close-in giant planet is justified.

Our model for the follow-up light curves used analytic formulae based on Mandel & Agol (2002) for the eclipse of a star by a planet, where the stellar flux is described by quadratic limb darkening. The limb-darkening coefficients are based on the SME results (Section 3.1) and the tables provided by Claret (2004) for the i band. The transit shape was parameterized by the normalized planetary radius R_p/R_⋆, the square of the impact parameter b², and the reciprocal of the half mid-ingress to mid-egress duration of the transit ζ/R_⋆. We chose these parameters because of their simple geometric meanings and the fact that they show negligible correlations (see Bakos et al. 2010). Our model for the HATNet data was the simplified "P1P3" version of the Mandel & Agol (2002) analytic functions (an expansion by Legendre polynomials), for the reasons described in Bakos et al. (2010). Following the formalism presented by Pál (2009), the RV curve was modeled as an eccentric Keplerian orbit with semi-amplitude K and Lagrangian orbital elements (k, h) = e × (cos ω, sin ω).

We assumed that there is a strict periodicity in the individual transit times. In practice, we assigned the transit number $N_{\rm{tr}}=0$ to the first high quality follow-up light curve gathered on 2009 September 10. The adjusted parameters in the fit were the first transit center observed with HATNet, T_c,−284, and the last high quality transit center observed with the FLWO 1.2 m telescope, T_c,+18. The transit center times for the intermediate transits were interpolated using these two epochs and the N_tr transit number of the actual event (Bakos et al. 2007b). The model for the RV data contains the ephemeris information through the T_c,−284 and T_c,+18 variables. Altogether, the eight parameters describing the physical model were T_c,−284, T_c,+18, R_p/R_⋆, b², ζ/R_⋆, K, k = ecos ω, and h = esin ω. Nine additional parameters were related to the instrumental configuration. These are the three instrumental blend factors B_inst,i of HATNet (one for each of the three fields), which account for possible dilution of the transit in the HATNet light curve due to the wide (20'' wide FWHM) point-spread function and possible crowding, the three HATNet out-of-transit magnitudes, M_0,HN,i, and three relative RV zero-points γ_rel,j (one each for the Keck, high-resolution FIES, and medium-resolution FIES observations).

We extended our physical model with an instrumental model that describes the systematic variations of the data. This was done in a similar fashion to the analysis presented in Bakos et al. (2010). The HATNet photometry has already been EPD- and TFA-corrected before the global modeling, so we only considered systematic correction of the follow-up light curves. We chose the "ELTG" method, i.e., EPD was performed in "local" mode with EPD coefficients defined for each night, and TFA was performed in "global" mode using the same set of stars and TFA coefficients for all nights. The underlying physical model was based on the Mandel & Agol (2002) analytic formulae, as described earlier. The five EPD parameters were the hour angle (characterizing a monotonic trend that linearly changes over time), the square of the hour angle, and the stellar profile parameters (equivalent to FWHM, elongation, and position angle). The functional form of the above parameters contained six coefficients, including the auxiliary out-of-transit magnitude of the individual events. The EPD parameters were independent for all four nights, implying 24 additional coefficients in the global fit. For the global TFA analysis, we chose 20 template stars that had good quality measurements for all nights and on all frames, implying an additional 20 parameters in the fit. Thus, the total number of fitted parameters was 17 (physical model) + 24 (local EPD) + 20 (global TFA) = 61.

The joint fit was performed as described in Bakos et al. (2010). We minimized the χ² in the parameter space by using a hybrid algorithm, combining the downhill simplex method (AMOEBA; see Press et al. 1992) with the classical linear least-squares algorithm. Uncertainties for the parameters were derived using the Markov Chain Monte Carlo method (see Ford 2006; Pál 2009).

The eccentricity of the system appeared as significantly non zero (k = −0.030 ± 0.003, h = −0.021 ± 0.006). The best-fit results for the relevant physical parameters are summarized in Table 4. We also list the RV "jitter," which is a quantity added in quadrature to the calculated RV measurement uncertainties in order to have χ²/dof = 1 from the RV data for the global fit. The "jitter" may be of astrophysical origin (pulsations, activity) or it may be in part due to instrumental effects.

Table 4. Orbital and Planetary Parameters

Parameter	Value
Light curve parameters
P (days).........	2.775960 ± 0.000003
Tc (BJD).........	2455027.59293 ± 0.00031
T14 (days)a.........	0.1276 ± 0.0013
T12 = T34 (days)a.........	0.0150 ± 0.0014
a/R⋆.........	7.17 ± 0.28
ζ/R⋆.........	17.73 ± 0.10
Rp/R⋆.........	0.1071 ± 0.0014
b2.........	0.193+0.063−0.069
b ≡ acos i/R⋆.........	0.439+0.065−0.098
i (deg).........	86.6 ± 0.7
RV parameters
K (m s−1).........	531.1 ± 2.8
kRVb.........	-0.030 ± 0.003
hRVb.........	-0.021 ± 0.006
e.........	0.036 ± 0.004
ω.........	214 ± 8°
RV jitter (m s−1).........	8.0
RV rms from fit (m s−1)......	10.0
Secondary eclipse parameters
Ts (BJD).........	2455028.929 ± 0.005
Ts,14.........	0.1234 ± 0.0020
Ts,12.........	0.0142 ± 0.0013
Planetary parameters
Mp (MJ).........	4.193 ± 0.094
Rp (RJ).........	1.289 ± 0.066
C(Mp, Rp)c.........	0.57
ρp ($\rm g\,cm^{-3}$).........	2.42 ± 0.35
a (AU).........	0.0413 ± 0.0004
log gp (cgs).........	3.80 ± 0.04
Teq (K).........	1626 ± 40
Θd.........	0.220 ± 0.011
〈F〉 (109$\rm erg\,s^{-1}\,cm^{-2}$)e......	1.58 ± 0.16

Notes. ^aT₁₄: total transit duration, time between first to last contact; T₁₂ = T₃₄: ingress/egress time, time between first and second or third and fourth contact. ^bLagrangian orbital parameters derived from the global modeling and primarily determined by the RV data. ^cCorrelation coefficient between the planetary mass M_p and radius R_p. ^dThe 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). ^eIncoming flux per unit surface area.

Download table as:  ASCIITypeset image

In addition, some auxiliary parameters (not listed in the table) are: T_c,−284 = 2454297.5154 ± 0.0009 (BJD), T_c,+18 = 2455135.8554 ± 0.0003 (BJD), γ_rel = 4.5 ± 3.8 m s⁻¹, γ_rel,FIEShi = −434.9 ± 4.2 m s⁻¹(FIES, high resolution), γ_rel,FIESmed = −442.9 ± 2.8 m s⁻¹(FIES, medium resolution), B_instr,123 = 0.77 ± 0.08, B_instr,163 = 0.93 ± 0.04, and B_instr,164 = 0.96 ± 0.04. Note that these gamma velocities do not correspond to the true center-of-mass RV of the system, but are only relative offsets. The true systemic velocity of the system is RV = −16.83 ± 0.19 km s⁻¹, as described in Section 2.2. The planetary parameters and their uncertainties can be derived by the direct combination of the a posteriori distributions of the light curve, RV, and stellar parameters. We found that the planet is fairly massive with M_p = 4.193 ± 0.094 M_J, and compact with radius of R_p = 1.289 ± 0.066 R_J, yielding a mean density of ρ_p = 2.42 ± 0.35 g cm⁻³. The final planetary parameters are summarized at the bottom of Table 4.