ORBITAL PHASE VARIATIONS OF THE ECCENTRIC GIANT PLANET HAT-P-2b

We determine the background level in each image from the region outside of a 10 pixel radius from the central pixel. This minimizes contributions from the wings of the star's point-spread function (PSF) while still retaining a substantial statistical sample. We iteratively trim 3σ outliers from the background pixels then fit a Gaussian to a histogram of the remaining pixel values to estimate a sky value. The background flux is 0.6% and 0.2% of the total flux in the science aperture for 3.6 and 4.5 μm observations, respectively. We correct for transient hot pixels by flagging pixels more than 4.5σ away from the median flux at a given pixel position across each set of 64 images then replacing flagged pixels by their corresponding position median value.

Several different methods exist to determine the stellar centroid on the Spitzer array such as flux-weighted centroiding (e.g., Knutson et al. 2008), Gaussian centroiding, and least asymmetry methods (e.g., Stevenson et al. 2010; Agol et al. 2010). We compare photometry calculated using both Gaussian and flux-weighted centroiding estimates and find that for the HAT-P-2 data flux-weighted centroiding produces the smallest scatter in the final light curve solutions. This is in contrast with the work by Stevenson et al. (2010) and Agol et al. (2010), which advocate Gaussian fits and least asymmetry methods to determine the stellar centroid for observations lasting less than 10 hr. We find that flux-weighted centroids give more stable position estimates over long periods, especially if the stellar centroid crosses a pixel boundary during the duration of the observation. For these data, we calculate the flux-weighted centroid for each background subtracted image using a range of aperture sizes from 2.0 to 5.0 pixels in 0.5 pixel increments. We find that the stellar centroid aperture sizes that best reduce the scatter in the final time series are 4.5 and 3.5 pixels for the 3.6 μm and 4.5 μm observations, respectively.

We estimate the stellar flux from each background subtracted image using circular aperture photometry. A range of aperture sizes from 2.0 to 5.0 pixels in 0.25 pixel increments were tested to determine the optimal aperture size for each observation. Additionally, we tested time-varying aperture sizes based on the noise pixel parameter described in Appendix A. We find that we obtain the lowest standard deviation of the residuals from our best-fit solution at 3.6 μm using a variable aperture size given by the square-root of the measured noise pixel value for each image (median aperture size of 2.4 pixels). For the 4.5 μm observations we find that a fixed 2.25 pixel aperture gives the lowest scatter in the final solution. We remove outliers from our final photometric data sets by discarding points more than 4.5σ away from a moving boxcar median 16 points wide. We find that using a narrower median filter or larger value for the outlier cutoff will miss some of the significant outliers. We also find that using a wider median filter or smaller value for the outlier cutoff will often selectively trim the points at the top and bottom of the "saw-tooth" pattern seen in the resulting photometry for the 3.6 and 4.5 μm observations (Figures 1 and 2), which is the result of intrapixel sensitivity variations discussed in Appendix B.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For our 5.8 μm data set, we determine the background level in each image using the same methodology as presented for the 3.6 and 4.5 μm data sets (Section 2.1). The background flux is 0.6% of the total flux in the science aperture. We tested both flux-weighted and Gaussian fits to the image methods of determining the stellar centroid. For the Gaussian fits to the image we first fit the image allowing both the x and y width of the Gaussian to vary. We then refit the same data set using a symmetric fixed width for the Gaussian that is the mean of the previous x and y widths. We find that Gaussian fits to the image give a lower standard deviation of the residuals in the final fits compared with flux-weighted centroids. This preference for the Gaussian centroiding method is likely the result of the short-time scale of these observations compared with the 3.6 and 4.5 μm observations and the less than 0.1 pixel change in the stellar centroid position over the course of the observation.

We estimate the stellar flux from each background subtracted image using circular aperture photometry for a range of apertures sizes from 2.0 to 5.0 pixels in 0.25 pixel increments. We find that an aperture size of 2.5 pixels gives the lowest standard deviation of the residuals in the final fits. We remove outliers in our data set more than 3.0σ away from a moving boxcar median 50 points wide. Figure 3 show our resulting photometry for the 5.8 μm observations. We test for possible intrapixel sensitivity variations using the methodology discussed in Appendix B, but find that including intrapixel sensitivity corrections actually degraded the precision of the final solution.

Zoom In Zoom Out Reset image size

Unlike the 3.6, 4.5, and 5.8 μm data, the 8.0 μm data were observed in the full-array (256 × 256 pixels) mode of the IRAC instrument. Our photometry determines the flux from the star in a circular aperture of variable radius, as well as the background thermal emission surrounding the star. We calculate two values of the background, using different spatial regions of the frame. First, we calculate the background that applies to the entire frame. Second, we isolate an annulus between 6 and 35 pixels from the star. In each region (entire frame, or annulus) we calculate the histogram of intensity values from the pixels within the region, and we fit a Gaussian to that histogram. We use the centroid of that Gaussian as the adopted background value. This method has the advantage of being insensitive to "hot" or otherwise discrepant values. We find that using the background values calculated from the region near the star, as opposed to the entire frame, produces a lower scatter in the residuals in our final solution.

We locate the center of the star by fitting a two-dimensional Gaussian to a 3 × 3 pixel median filtered version of the frame; we find that this method produces slightly better photometric precision than using a center-of-light algorithm. We estimate the stellar flux from each background subtracted image using circular aperture photometry for a range of apertures sizes from 2.0 to 5.0 pixels in 0.25 pixel increments. We find that an aperture size of 3.5 pixels gives the lowest standard deviation of the residuals in the final fits. We remove outliers in our data set more than 3.0σ away from a moving boxcar median 50 points wide. Figure 4 show our resulting photometry for the 8.0 μm observations. We test for possible intrapixel sensitivity variations using the methodology discussed in Appendix B, but find that including intrapixel sensitivity corrections actually degraded the precision of the final solution.

Zoom In Zoom Out Reset image size

Our data exhibit a ramp-like increase, or decrease, in the flux values at the start of each observation and after each break in the 3.6 and 4.5 μm data for downlink. This ramp in the flux values has been previously noted for IRAC 8.0 μm observations (e.g., Knutson et al. 2007, 2009b; Agol et al. 2010) as well as 3.6 and 4.5 μm observations (e.g., Beerer et al. 2011; Todorov et al. 2012). While the precise nature of this ramp is not well-understood, it can be at least partially attributed to thermal settling of the telescope at a new pointing, which contributes an additional drift in position. However, we see that the ramp often persists beyond this initial drift in position, and we therefore speculate that there is an additional component analogous to the effect observed at long wavelengths due to charge-trapping in the array. We tested several functional forms to describe this ramp, including quadratic, logarithmic, and linear functions of time, but find that functional form that gives us the best correction is given by the formulation presented in Agol et al. (2010):

$$\begin{equation} F^{\prime }/F=1\pm a_1 e^{-t/a_2} \pm a_3 e^{-t/a_4} \end{equation} \tag{ 1 }$$

where F' is the measured flux, F is the stellar flux incident on the array, t is the time from the start of the observation in days, and a₁–a₄ are free parameters in the fit.

We find that in the 3.6 and 4.5 μm data the asymptotic shape of the ramp converges on timescales less than an hour for most of the data segments. Instead of including parameters in our fits for this ramp behavior, we instead elect to simply trim the first hour at the start of each observation and subsequent downlinks. This reduces the complexity of our fits while avoiding possible correlations between the shape of the underlying phase curve of HAT-P-2b and the ramp function. We find that the ramp persists beyond the one hour mark at the start of the 3.6 μm observations, which affects the shape of the first observed secondary eclipse and therefore include a ramp correction for that data segment. For the 5.8 μm data, we find that the flux asymptotically decays from the start of the observation (Figure 3). This decrease in the 5.8 μm flux is also apparent in the derived background flux values and has been previously noted by studies such as Stevenson et al. (2012). As seen in Figure 4, the 8.0 μm data exhibit the well know asymptotic increase in flux (e.g., Knutson et al. 2007, 2009b; Agol et al. 2010).

In order to select the best functional form for the ramp correction in each data set we use the Bayesian Information Criterion (BIC), defined as

$$\begin{equation} {\rm BIC}=\chi ^2+k\ln (n) \end{equation} \tag{ 2 }$$

where k is the degrees of freedom in the fit and n is the total number of points in the fit (Liddle 2007). This allows us to determine if we are "over-fitting" the data by including additional parameters to describe the ramp correction in our fits. We find that for the 3.6 μm data set using only a single exponential gives us the lowest BIC value. For the 8.0 μm data, using all the terms (a₁–a₄) in Equation (1) with no trimming of the data near the start of the observation gives us the lowest BIC value. For the 5.8 μm data, we find using a single decaying exponential gives the lowest BIC if we trim data within 30 minutes of the start of the observation. Trimming significantly more or less than this amount ether gives a higher BIC or reduces the amount of out of eclipse data to less than the eclipse duration.

We model our transit and eclipse events using the equations of Mandel & Agol (2002) modified to account for the orbital eccentricity, e, and argument of periapse, ω, of the HAT-P-2 system. For an eccentric system the normalized separation of the planet and star centers, z, is given by

$$\begin{equation} z=\frac{r(t)}{R_{\star }}\sqrt{1-(\sin {i}\times \sin {(\omega +f(t))})^2} \end{equation} \tag{ 3 }$$

where r(t) is the radial planet–star distance as a function of time, R_* is the stellar radius, i is the orbital inclination of the planet, and f(t) is the true anomaly as a function of time. The r(t)/R_⋆ term in Equation (3) is calculated as

$$\begin{equation} \frac{r(t)}{R_{\star }}=\frac{a}{R_{\star }}\frac{1-e^2}{1+e\cos {(f(t))}}. \end{equation} \tag{ 4 }$$

The true anomaly angle (f) that appears in both Equations (3) and (4) is determined using Kepler's equation (Murray & Dermott 1999) and is a function of e, ω, and the orbital period of the planet, P. Because of the degeneracies that exist between e, ω, and P in determining f(t) we elect not to use P as a free parameter in our fits. Instead we fix P to the value reported in Pál et al. (2010), 5.6334729 d. We further minimize correlations in our solutions for e and ω by solving for the Lagrangian orbital elements k ≡ ecos ω and h ≡ esin ω.

In addition to the parameters that define the orbit of HAT-P-2b, we solve for the fractional planetary radius, R_p/R_⋆. Because R_p/R_⋆ < 0.1 for HAT-P-2b (Pál et al. 2010; Bakos et al. 2007a), there exists a strong correlation between i and a/R_⋆ in the transit solution (Winn et al. 2007a; Pál 2008). We instead solve for the parameters b² and ζ/R_⋆ as suggested by Bakos et al. (2007b) and Pál et al. (2010) where

$$\begin{equation} b=\frac{a}{R_{\star }}\frac{1-e^2}{1+h}\cos i \end{equation} \tag{ 5 }$$

and

$$\begin{equation} \frac{\zeta }{R_{\star }}=\frac{a}{R_{\star }}\frac{2\pi }{P}\frac{1}{\sqrt{1-b^2}}\frac{1+h}{\sqrt{1-e^2}}. \end{equation} \tag{ 6 }$$

For the transit portion of the light curves we use four-parameter nonlinear limb-darkening coefficients for each bandpass calculated by Sing (2010), where we assume a stellar atmosphere with T_eff = 6290 K, log (g) = 4.138, and [Fe/H] = +0.14 (Pál et al. 2010). For the secondary eclipse portion of the light curves we treat the planet as a uniform disk and scale the ingress and egress to match the amplitude of the phase curve (Section 2.6), which varies significantly over the duration of the eclipse and is therefore poorly approximated by a constant value. The secondary eclipse depth is defined by the average of the ingress and egress amplitudes.

For the 3.6 and 4.5 μm data sets, we use nine free parameters to constrain the properties of the planetary orbit, transit, and secondary eclipse. The 8.0 μm observations only include one secondary eclipse and therefore only require eight free parameters to constrain the same system properties. The 5.6 μm data set includes only the secondary eclipse portion of the light curve. We therefore elect to only allow the secondary eclipse depth and timing to vary for the 5.6 μm data and fix a/R_⋆, i, e, ω, and R_p/R_⋆ to the average values from the 3.6, 4.5, and 8.0 μm data sets.

The functional form of the phase curve for a planet on an eccentric orbit is not well defined. Unlike close-in, tidally locked planets on circular orbits, eccentric planets experience time-variable heating and non-synchronous rotation rates. Previous studies by Langton & Laughlin (2008) and Cowan & Agol (2011) have investigated theoretical light curves for planets on eccentric orbits using two-dimensional hydrodynamic simulations and semi-analytic model atmospheres, respectively. We also developed a three-dimensional atmospheric model for HAT-P-2b that couples radiative and dynamical processes to further investigate possible phase curve functional forms that will be presented in a future paper. The functional forms for the phase curves described here all provide a reasonable fit to the theoretical light curves presented in Langton & Laughlin (2008), Cowan & Agol (2011), and from our three-dimensional atmospheric model.

To first order, the flux from the planet is proportional to the inverse square of the distance between the planet and host star, r(t). This assumes that the planet has a constant albedo and responds instantaneously and uniformly to changes in the incident stellar flux. We know that there must exist a lag in the peak of the incident stellar flux and the peak of the planet's temperature since atmospheric radiative timescales are finite (Langton & Laughlin 2008; Iro & Deming 2010; Lewis et al. 2010; Cowan & Agol 2011). Our first functional form for the phase variation of HAT-P-2b is a simple 1/r(t)² with a phase lag:

$$\begin{equation} F(f)=F_0+c_1(1+\cos (f-c_2))^2 \end{equation} \tag{ 7 }$$

where f is the true anomaly and c₁–c₂ are free parameters in the fit. The 1 + cos (f − c₂) is a simplified form of the denominator of Equation (4) for r(t). We also test a simpler form of Equation (7) given by

$$\begin{equation} F(f)=F_0+c_1\cos (f-c_2) \end{equation} \tag{ 8 }$$

where c₁–c₂ are free parameters in the fit. This functional form of the phase curve is similar to a simple sine or cosine of the orbital phase angle, λ or ξ, used in the case of a circular orbit.

We also find that a Lorentzian function of time provides a reasonable representation of the expected shape of the orbital phase curve for HAT-P-2b. This is not surprising given that we expect the flux from the planet to vary as 1/r(t)² with a time lag between the minimum of the planetary distance and the peak of the planetary flux. We test both symmetric and asymmetric Lorentzian functions of the time from periapse passage, t, given by

$$\begin{equation} F(t)=F_0+\frac{c_1}{u(t)^2+1} \end{equation} \tag{ 9 }$$

where

$$\begin{equation} u(t)=(t-c_2)/c_3 \end{equation} \tag{ 10 }$$

in the symmetric case and

$$\begin{equation} u(t)= \left\lbrace \begin{array}{@{} ll} (t-c_2)/c_3 & {\rm if\; {t < c_2;}} \\ (t-c_2)/c_4 & {\rm if \; { t > c_2 }} \end{array} \right. \end{equation} \tag{ 11 }$$

in the asymmetric case. In Equations (9), (10), and (11), c₁–c₄ are free parameters in the fit. The Lorentzian functional form for the phase curve is especially useful since the c₂ parameter gives the offset between the time of periapse passage and the peak of the planet's flux and the c₃–c₄ parameters gives an estimate of relevant atmospheric timescales.

We also find that the preferred functional form for the phase curve of planets on circular orbits described in Cowan & Agol (2008) provides a reasonable fit to our theoretical light curves if the orbital phase angle, ξ, is replaced by the true anomaly, f. In this case

$$\begin{eqnarray} F(\theta) &=& F_0+c_1\cos (\theta)+c_2\sin (\theta)\nonumber\\ && +c_3\cos (2\theta)+c_4\sin (2\theta) \end{eqnarray} \tag{ 12 }$$

where c₁–c₄ are free parameters in the fit and θ = f + ω + π such that transit occurs at θ = −π/2 and secondary eclipse occurs at θ = π/2.

In addition to the functional forms for the phase curve presented above, we also test a flat phase curve, F(t) = F₀, to make sure that we are not over-fitting the data with the phase curve parameters. In all cases we tie the F₀ parameter to the secondary eclipse depth to give the appropriate average value of the phase curve during secondary eclipse. We also set any portion of the phase curve that falls below the secondary eclipse depth to zero. During secondary eclipse we no longer see flux from the planet, only the star, therefore the combined star and planetary flux should always be greater than or equal to the flux at the bottom of the secondary eclipse.

We choose the optimal phase curve solution to be the one that gives the lowest BIC value (Equation (2)) using a Levenberg–Marquardt nonlinear least-squares fit to the data (Markwardt 2009). For the 3.6 and 4.5 μm data sets, we also compare the BIC values from each of the three data segments separated by the downlink periods to make sure that the solution is robust at all orbital phases. For the 3.6, 4.5, and 8.0 μm data sets, we find that all of the functional forms presented Section 2.6 provide a significant improvement in the BIC values over the "flat" phase curve BIC values (BIC_3.6 μm = 1,282,575; BIC_4.5 μm = 1,459,141; BIC_8.0 μm = 11,183). For the 3.6 μm data we obtain the lowest BIC value with a functional form for the phase curve defined by either Equation (8) or Equation (12) with c₃–c₄ fixed at zero (BIC = 1,278,521). This is not surprising since basic trigonometric identities make Equation (12) equivalent to Equation (8) if c₃–c₄ can be assumed to be zero. For the 4.5 μm data we obtain the lowest BIC value with a functional form for the phase curve given by Equation (12) with the c₁ and c₃ terms fixed at zero (BIC = 1,458,842).

The phase curve model given by Equation (12) also gives us a reasonable improvement in the BIC for our 8.0 μm data set. However, we find that the second harmonic in Equation (12) is highly degenerate with our ramp correction at 8.0 μm such that we cannot reliably constrain the significance of this harmonic in the fit. We instead use the phase curve functional form given by Equation (7), which gives us the lowest BIC for the 8.0 μm data set (BIC = 10,945). For the 5.8 μm data, which only span 10 hours, we find no statistical difference between solutions with and without parameterizations for planetary phase variations.

HAT-P-2 is a rapidly rotating F star (vsin i = 20.7 km s⁻¹). Measurements of the Ca ii H and K lines of HAT-P-2 by Knutson et al. (2010) do not detect any significant emission in the line cores, which would suggest a chromospherically quiet stellar host for HAT-P-2b. However, this indicator provides relatively weak constraints on chromospheric activity for F stars, which have strong continuum emission in the wavelengths of the Ca ii H and K lines. Ground-based monitoring of HAT-P-2 in the Strömgren b and y bands over a period of more than a year indicates that it varies by less than 0.13% at visible wavelengths. We would expect the star to vary by substantially less than this amount in the infrared, where the flux contrast of spots and other effects is correspondingly reduced. These observations rule out both periodic variability that could be associated with the rotation rate of HAT-P-2 or longer-term trends (G. Henry 2010, private communication). Previous observations find the rotation rate of the star to be on the order of 3.8 d (Winn et al. 2007b) based on the line-of-sight stellar rotation velocity (vsin i_⋆) and assuming sin i_⋆ = 1, which is roughly 0.6 times the orbital period of HAT-P-2b. We do not expect stellar variability to be significant on the timescales of our observations, but for the sake of completeness we also check this assumption directly using our Spitzer data.

We employ two models for stellar variability. The first model is a simple linear function of time given by

$$\begin{equation} F_0(t) = d_1 t \end{equation} \tag{ 13 }$$

where t is measured from the predicted center of transit in each observation and d₁ is a free parameter in the fit. The second model we test has the form

$$\begin{equation} F_0(t)=d_1\sin \,((2\pi /d_2)t-d_3) \end{equation} \tag{ 14 }$$

where t is measured from the predicted center of transit in each observation and d₁ − d₃ are free parameters in the fit. Equation (14) attempts to capture stellar variability that is associated with star spots that rotate in and out of view with the d₂ parameter representing the rotation rate of the star.

We find that the linear model for stellar variability does not improve the χ² of the fits significantly and in fact increases the BIC (Equation (2)). Inclusion of the stellar model given by Equation (14) does improve both the χ² and BIC in our fits. However, we find that the solutions for the "sine-curve" model for stellar variability are often degenerate with our models for the phase curve and the residual ramp at the start of the 3.6 μm observations. We also find that the stellar rotation rate predicted by the d₂ term in Equation (14) differs significantly between the 3.6 and 4.5 μm observations with a rotation rate of 4.3 d preferred for the 3.6 μm data and 3.9 d preferred for the 4.5 μm data. Although these predicted rotation periods are near to the expected 3.8 d rotation period of HAT-P-2, the amplitudes of the predicted stellar variations in the mid-infrared seem spurious. The amplitudes of the predicted stellar flux variations at 3.6 and 4.5 μm are on the order of ∼0.1%, which is comparable to our upper limit on stellar variability at visible wavelengths (0.13%). We would expect star spots to display much larger variations in the visible than at mid-infrared wavelengths. The amplitudes of these stellar variations are also similar to the amplitudes of our predicted phase variations. We therefore conclude that our best-fit solutions for the stellar variability are physically implausible, and likely the result from a degeneracy with the other terms in our fits. As we have no convincing evidence for stellar variability, we assume that the star's flux is constant in our final fits.

Since the discovery of HAT-P-2b by Bakos et al. (2007a), the California Planet Search (CPS) team has continued to obtain regular RV measurements of this system using the HIRES instrument (Vogt et al. 1994) on Keck. We used the CPS pipeline (see, e.g., Howard et al. 2009) to measure precise RVs from the high-resolution spectra of HAT-P-2 using a superposed molecular iodine spectrum as a Doppler reference and PSF calibrant. Previous studies of this system (Bakos et al. 2007a; Winn et al. 2007b; Pál et al. 2010) have included HIRES RV measurements through 2008 May (BJD = 2454603.932112). Here we present 16 additional RV measurements from the HIRES instrument on Keck for a total of 71 data points spanning six years (Table 1). We exclude measurements obtained during transit (Winn et al. 2007b) from our fits in order to avoid measurements affected by the Rossiter–McLaughlin effect.

We simultaneously fit for the RV semi-amplitude (K), zero-point (γ), and long term velocity drift ($\dot{\gamma }$) along with the orbital, transit, eclipse, and phase parameters in the 3.6, 4.5, and 8.0 μm data sets separately. We also test for possible curvature in the RV signal ($\ddot{\gamma }$), but find that our derived values for $\ddot{\gamma }$ were consistent with zero. As discussed in Winn (2011), the transit to secondary eclipse timing strongly constrains the ecos ω term in our fits, but the esin ω term is better constrained by the inclusion of these RV data. Rapidly rotating stars are known to have an increased scatter in their RV velocity distribution beyond the reported internal errors (as discussed in Bakos et al. 2007a; Winn et al. 2007b). We therefore estimate a stellar jitter term (σ_jitter) as described in Section 3.

We find that our estimates for the orbital and RV parameters listed in Table 2 fall within the range of the values from previous studies presented in Table 3. We note that there is often a more than 3σ discrepancy between orbital parameters for the HAT-P-2 system presented in previous studies (Table 3). We conclude that either previous studies underestimate their error bars for the discrepant parameters, or the planet's orbital properties are varying in time. If we compare the orbital parameters we derived from the 3.6, 4.5, and 8.0 μm data sets, we find that our estimates are within 3σ of each other. The previous studies presented in Table 3 included only ground-based transit and RV data. By the inclusion of secondary eclipse in our data we were able to improve the estimate of the orbital eccentricity of HAT-P-2b by an order of magnitude over the value presented in Pál et al. (2010).

Table 3. Results from Previous HAT-P-2 Studies

Parameter	Bakos et al. (2007a)	Winn et al. (2007b)	Loeillet et al. (2008)	Pál et al. (2010)
Transit parameters
Rp/R⋆	0.0684 ± 0.0009	0.0681 ± 0.0036a	0.06891$^{+0.00090}_{-0.00086}$	0.07227 ± 0.00061
i (°)	>84.6	>86.8	90.0$^{+0.85}_{-0.93}$	86.72$^{+1.12}_{-0.87}$
a/R⋆	9.77$^{+1.10}_{-0.02}$	9.90 ± 0.39a	10.28$^{+0.12}_{-0.19}$	8.99$^{+0.39}_{-0.41}$
Tc (BJD −2,400,000)	54212.8563 ± 0.0007	54212.8565 ± 0.0006	...	54387.49375 ± 0.00074
T14 (d)	0.177 ± 0.002	...	...	0.1787 ± 0.0013
T12 (d)	0.012 ± 0.002	...	...	0.0141$^{+0.0015}_{-0.0012}$
Secondary eclipse parameters
ϕsec	0.1847 ± 0.0055a	0.1969 ± 0.0040a	0.1896 ± 0.0016a	0.1868 ± 0.0019
T14 (d)	...	...	...	0.1650 ± 0.0034
Tc (BJD −2,400,000)	...	...	...	54388.546 ± 0.011
Orbital and RV parameters
P (d)	5.63341 ± 0.00013	5.63341 (fixed)	5.63341 (fixed)	5.6334729 ± 0.0000061
ecos ω	...	...	...	−0.5152 ± 0.0036
esin ω	...	...	...	−0.0441 ± 0.0084
e	0.520 ± 0.010	0.501 ± 0.007	0.5163$^{+0.0025}_{-0.0023}$	0.5171 ± 0.0033
ω (°)	179.3 ± 3.6	187.4 ± 1.6	189.92$^{+1.06}_{-1.20}$	185.22 ± 0.95
Tp (BJD −2,400,000)	54213.369 ± 0.041	...	54213.4798$^{+0.0053}_{-0.0030}$	...
K (m s−1)	1011 ± 38	883 ± 57a	966.9 ± 8.3	983.9 ± 17.2
Planetary parameters
Mp (MJ)	9.04 ± 0.50	8.04 ± 0.40	8.62$^{+0.39}_{-0.55}$	9.09 ± 0.24
Rp (RJ)	0.982$^{+0.038}_{-0.105}$	0.98 ± 0.04	0.951$^{+0.039}_{-0.053}$	1.157$^{+0.073}_{-0.062}$
ρp (g cm−3)	11.9$^{+4.8}_{-1.6}$	10.60 ± 0.55a	12.5$^{+2.6}_{-3.6}$	7.29 ± 1.12
gp (m s−2)	227$^{+44}_{-16}$	207 ± 20a	237$^{+30}_{-41}$	168 ± 17
a (AU)	0.0677 ± 0.0014	0.0681 ± 0.0014a	0.0677$^{+0.0011}_{-0.0017}$	0.06878 ± 0.00068
Noise parameters
σjitter (m s−1)	60	31	17	...

Note. ^aParameter value not directly quoted in reference, but calculated from quoted parameter values and errors.

Download table as:  ASCIITypeset image

From our orbital and RV parameters we can estimate the mass of HAT-P-2b using

$$\begin{equation} M_p=\frac{2\pi }{P}\frac{K\sqrt{1-e^2}}{G\sin i}\left(\frac{a}{R_{\star }}\right)^2 R_{\star }^2, \end{equation} \tag{ 16 }$$

where the orbital period (P) and stellar radius (R_⋆) are assumed to be the values presented in Pál et al. (2010), 5.6334729 ± 0.0000061 d and $1.64^{+0.09}_{-0.08}R_{\odot }$, respectively. Values for the RV semi-amplitude (K), eccentricity (e), inclination (i), and normalized semi-major axis (a/R_⋆) are taken as the error-weighted average of the values from the 3.6, 4.5, and 8.0 μm observations, which are presented in Table 4. We estimate the mass of HAT-P-2b to be 8.00 ± 0.97 M_J, which is within 1σ of the previous estimates of M_p for HAT-P-2b (Table 3).

Table 4. HAT-P-2b Parameters from a Weighted Average of the Values from the 3.6, 4.5, and 8.0 μm Fits

Parameter	Value
Rp/R⋆	0.06933 ± 0.00045
Etransit (BJD)	2455288.84923 ± 0.00037
a/R⋆	8.70 ± 0.19
i (°)	85.97$^{+0.28}_{-0.25}$
e	0.50910 ± 0.00048
ω	188.09 ± 0.39
ϕsec	0.19253 ± 0.00018
Esec (BJD)	2455289.93211 ± 0.00066
Eperiapse (BJD)	2455289.4721 ± 0.0038
Mp (MJ)	8.00 ± 0.97
Rp (RJ)	1.106 ± 0.061
ρp (g cm−3)	7.3 ± 1.6
gp (m s−2)	162 ± 27
a (AU)	0.0663 ± 0.0039

Download table as:  ASCIITypeset image

Our RV data span a period of nearly six years, which allows us to test for long-term trends in the RV signal. We find a non-zero value for the linear term in our RV fit ($\dot{\gamma }$), which indicates the presence of a second body orbiting in the system. If we assume that this second companion to HAT-P-2 ("c") is on a circular orbit and is much less massive than its host star, we can relate $\dot{\gamma }$ to the mass and orbital semi-major axis of "c" by the expression presented in Winn et al. (2009)

$$\begin{equation} \dot{\gamma }\approx \frac{G M_c \sin i_c}{a_c^2}. \end{equation} \tag{ 17 }$$

Figure 11 shows the range of values for M_csin i_c and a_c that are allowed for the HAT-P-2 system given that $\dot{\gamma }=-0.0886\pm 0.0030$. We know the orbital period of "c" must be significantly longer than six years since we do not detect any significant curvature in our long-term RV trend, so we set a lower limit of M_csin i_c ∼ 15 M_J and a_c ∼ 10 AU based on an orbital period for "c" of 24 years. We further employ adaptive optics (AO) imaging to search for "c" and set an upper limit on the values of M_csin i_c and a_c.

Zoom In Zoom Out Reset image size

4.2.1. Adaptive Optics Imaging

In an attempt to directly image the body responsible for causing the linear RV trend, we observed HAT-P-2 on 2012 May 29 using NIRC2 (PI: Keith Matthews) and the Keck II AO system at Mauna Kea (Wizinowich et al. 2000). Our observations consist of dithered images taken with the K' (λ_c = 2.12 μm) filter. Using the narrow camera setting, which provides fine spatial sampling of the NIRC2 PSF (10 mas pixel⁻¹), we acquired nine frames each having 13.2 s of on-source integration time. The seeing was estimated to be 04 at visible wavelengths at the time of observations using AO wavefront sensor telemetry. We note that significant wind-shake degraded the quality of correction for some images but only by a marginal amount.

The data were processed using standard techniques to correct for hot pixels, remove background radiation from the sky and instrument optics, flat-field the array, and align and co-add individual frames. Figure 12 shows the final reduced AO image along with the corresponding (10σ) contrast levels achieved as a function of angular separation. No companions were detected in raw or processed frames.

Zoom In Zoom Out Reset image size

We can use the limits from a non-detection to rule out the presence of companions as a function of M_csin (i_c) and a_c. Using HAT-P-2's parallax distance of 119 ± 8 pc and estimated age of 2.6 ± 0.5 Gyr as determined by Pál et al. (2010) through a combined isochrone, light curve, and spectroscopic analysis, we find that our diffraction-limited observations are sensitive to companions on the hydrogen-fusing boundary for separations beyond ≈1''. Interior to this region, the combination of imaging and RV data eliminates most low-mass stars, though late-type M-dwarf tertiaries located at ∼40 AU could cause the long-term Doppler drift yet simultaneously evade direct detection (Figure 11).

4.2.2. Orbital Evolution

The likely presence of an M/L/T/Y dwarf at an orbital distance, a_c, of 10 to 40 AU from HAT-P-2b lends credence to the possibility that HAT-P-2b owes is current orbit to a history of Kozai cycling (Kozai 1962) combined with long-running orbital decay generated by tidal friction (Eggleton et al. 1998; Wu & Murray 2003; Fabrycky & Tremaine 2007). At first glance, this combination of Kozai Cycling and Tidal Friction (KCTF) seems more plausible than disk migration for producing the current orbital configuration of HAT-P-2b and indeed, long-distance disk migration for HAT-P-2b seems quite problematic. The protostellar disk itself would have needed to be exceptionally massive and long lived, and an ad hoc mechanism must be invoked to explain the large observed orbital eccentricity of HAT-P-2b. Furthermore, HAT-P-2b has planetary and orbital parameters that place it significantly outside the well-delineated population of "conventional" hot Jupiters with P ∼ 3 days, M ∼ M_Jup, and e ∼ 0.

In the context of the KCTF process, HAT-P-2b is envisioned to have formed well outside its current orbit, perhaps via gravitational instability in the original protostellar disk. KCTF can operate if the mutual orbital inclination between companion "c" and planet b was larger (or better, substantially larger) than the Kozai critical angle, $i_c = \arccos [(3/5)^{1/2}]\sim 40^{\circ }$ (assuming M_c ≫ M_b, and an initially circular orbit for b). In the event that KCTF did operate, planet b experienced periodic cycling between successive states of high orbital eccentricity and high mutual inclination. The timescale for these cycles was

$$\begin{equation} \tau _{\rm Kozai}=\frac{2 P_{c}^{2}}{3\pi P_{\rm b}}\frac{M_{\star }+M_{\rm b}+M_{\rm c}}{M_{\rm c}} \left(1-e_{\rm c}^2 \right)^{3/2}. \end{equation} \tag{ 18 }$$

With plausible values of P_c = 25 years, e_c = 0.5, and M_c = 20 M_Jup, τ_Kozai = 300 Kyr. By contrast, the general relativistic precession rate for HAT-P-2b is currently

$$\begin{equation} \dot{\omega} _{\rm GR}=\frac{3 G^{3/2}(M_{\star }+M_{b})^{3/2}}{a_{\rm b}^{5/2}c^2 \left(1-e_{b}^2\right)}, \end{equation} \tag{ 19 }$$

which generates a full 2π circulation of the apsidal line in 20 Kyr (τ_GR). Because τ_Kozai ≫ τ_GR, the magnitude of the Kozai cycling is strongly suppressed at present.

To a good approximation, tidal friction is currently the dominant contributor to the orbital evolution of HAT-P-2b, and it is plausible that the current state of the system has undergone dissipative evolution from an earlier epoch where τ_Kozai = τ_GR. If we assume that the tidal evolution has roughly conserved HAT-P-2b's periastron distance, d_peri = a_b(1 − e_b) = 0.033 AU, then in order for τ_Kozai = τ_GR implies that HAT-P-2b has evolved to its current orbit from an orbit with an initial period, $P_{o_{{\rm min}}}$, given by

$$\begin{equation} P_{o_{{\rm min}}}= \left(1-e_{\rm c}^{2} \right)^{3/4}\frac{P_{\rm c}M_{\star }}{c} \left(\frac{G}{M_{\rm c}d_{\rm peri}} \right)^{1/2}, \end{equation} \tag{ 20 }$$

or

$$\begin{equation} P_{o_{{\rm min}}}\sim 33\,{\rm d}\,\, \left(\frac{P_{\rm c}}{25\, {\rm yr}} \right) \left(\frac{20\,M_{\rm Jup}}{M_{\rm c}} \right)^{1/2}. \end{equation} \tag{ 21 }$$

Additionally, the RV-derived constraint on the unseen companion "c" indicates that $M_{\rm c}\sim P_{\rm c}^{4/3}$, which allows us to simply write

$$\begin{equation} P_{o_{{\rm min}}} \sim 33\,{\rm d}\, \left(\frac{P_{\rm c}}{25 \,{\rm yr}} \right)^{1/3}. \end{equation} \tag{ 22 }$$

Given that direct imaging constrains the maximum period for companion "c" to be of order 250 years, the KCTF scenario gives a lower limit on the possible initial periods for HAT-P-2b to be $30\,{\rm days}<P_{o_{{\rm min}}}<60\,{\rm days}$.

The KCTF process delivers planets into orbits in which the planetary orbital angular momentum and the stellar spin angular momentum are initially largely uncorrelated. As mentioned earlier in the text, initial measurements of the projected spin-orbit misalignment angle, λ, suggested that HAT-P-2b's orbit lies in the stellar equatorial plane. Recently, however, an reassessment by Albrecht et al. (2012) reports λ = 9° ± 5°, indicating a modest projected misalignment. In addition, HAT-P-2b's parent star is almost precisely on the T_eff = 6250 K boundary at which stars empirically appear to be able to maintain misalignment over multi-Gyr time scales (Albrecht et al. 2012).

We calculate ephemeris for the HAT-P-2 system given the center of transit and eclipse times presented in Table 2 as

$$\begin{equation} T_c(n)=T_c(0)+n\times P \end{equation} \tag{ 23 }$$

where T_c is the predicted transit, eclipse, or periapse time and n is the number of orbits that have elapses since T_c(0). From our transit data we calculate T_c(0) = 2455288.84923 ± 0.00037 BJD and P = 5.6334754 ± 0.0000026 d, which is consistent with the orbital period for HAT-P-2b presented in Pál et al. (2010). We calculate T_c(0) = 2455289.93211 ± 0.00066 BJD and P = 5.6334830 ± 0.0000086 d from our secondary eclipse data and T_c(0) = 2455289.4721 ± 0.0038 BJD and P = 5.633479 ± 0.000026 d from our estimated times of periapse passage. The orbital periods estimated from the secondary eclipse and periapse passage timings are within 1σ of the value derived from our transit timings. Figure 13 presents our transit, secondary eclipse, and periapse times compared to our derived constant ephemeris. We define the center of transit/eclipse to occur when the projected planet–star distance given by Equation (3) is minimized. The definition of transit/eclipse center is not always clearly stated in studies for eccentric systems. Pál et al. (2010) estimate a difference between RV and photometric transit centers for the HAT-P-2 system of nearly two minutes. Some of the spread in the measured versus predicted transit times could be accounted for by inconsistent definitions of transit center.

Zoom In Zoom Out Reset image size

The separate visits to HAT-P-2b for the 3.6, 4.5, 5.8, and 8.0 μm observations allow constraints to be placed on possible transit timing variations (TTVs) for HAT-P-2b. Previous transit timing measurements for this system (Bakos et al. 2007a; Winn et al. 2007b; Pál et al. 2010) had relatively low timing resolution (Δt ∼ 50 → 500 s) and appeared to be consistent with a constant ephemeris (Pál et al. 2010). Somewhat surprisingly, the Spitzer data from orbits 0 and 83 appear to be inconsistent with a constant ephemeris at the ∼3.5σ level, with a typical deviation of 150 s. The deviations are anti-correlated between primary transit and secondary eclipse, and they switch sign on the two successive visits. This could be due to either an astrophysical cause, or an as-yet unmodeled systematic error. In the sections below we consider two potential astrophysical explanations for the TTVs: a planetary satellite or an external perturber.

4.3.1. TTVs Generated by a Planetary Satellite

HAT-P-2b is expected to be in pseudo-synchronous rotation, in which its spin period is close to the orbital frequency in the vicinity of periapse. Given P_orbit = 5.63347 d, Hut's (1981) treatment gives

$$\begin{eqnarray} &&P_{\rm spin}=\left[\frac{\frac{5 e^6}{16}+\frac{45 e^4}{8}+\frac{15 e^2}{2}+1}{(1-e^2)^{3/2} \left(\frac{3 e^4}{8}+3 e^2+1\right)}\right]^{-1}P_{\rm orbit} = 1.89\,{\rm d}. \nonumber\\ \end{eqnarray} \tag{ 24 }$$

In order to maintain orbital dynamical stability, a prospective moon for HAT-P-2b must have an orbital semi-major axis, a_sat, such that a_sat ≲ F_critR_Hill, where R_Hill is the Hill sphere radius at periapse,

$$\begin{equation} R_{\rm Hill}\sim a_{\rm pl}(1-e) \left(\frac{M_{\rm pl}}{3\,M_{\star }} \right)^{1/3} = 0.0042\, {\rm AU}, \end{equation} \tag{ 25 }$$

and F_crit ≲ 0.5 (see, e.g., Barnes & O'Brien 2002 for a discussion of satellite stability for planets with low orbital eccentricity). At distance a_sat < F_critR_Hill from HAT-P-2b, the orbital period of the satellite is P_sat = 0.385 d, which is substantially shorter than the P_spin = 1.89 d spin period of the planet. Tidal evolution will therefore cause the satellite's orbit to gradually spiral in to meet the planet (Sasaki et al. 2012).

Based on the fairly uniform satellites-to-planet mass ratios exhibited by the regular satellites of the Jovian planets in the solar system, and in the absence of any concrete information regarding exomoons, it is not unreasonable to expect M_sat = M_pl/10⁴ ∼ 0.28 M_⊕. The characteristic time for a satellite's orbital decay, assuming an equilibrium theory of tides with a frequency-independent tidal quality factor, Q_p, is (Murray & Dermott 1999)

$$\begin{equation} \tau = \frac{2}{13}(F_{{\rm crit}} R_{\rm Hill})^{13/2}\frac{Q_{\rm pl}}{3k_{2p}M_{\rm s}R_{\rm pl}^5} \left(\frac{M_{\rm pl}}{G} \right)^{1/2}. \end{equation} \tag{ 26 }$$

Adopting Q_pl = 10⁶, and k_2p = 0.1 we find τ ∼ 2 Gyr, which is comparable to the estimated stellar age, τ_⋆ = 2.7 Gyr (Pál et al. 2010). The small value for the apsidal Love number, k_2p, stemming from the fact that HAT-P-2b at 8 M_Jup is quite centrally condensed in comparison to Jupiter (see Bodenheimer et al. 2001; see Batygin et al. 2009 for a discussion of the relation between k_2p and interior models). Given the substantial range of possible values for Q_pl, it would therefore not be unreasonable to find that HAT-P-2b harbors a fairly massive moon.

A moon with the above qualities would produce TTVs of magnitude (Kipping 2009)

$$\begin{equation} \delta _{{\rm TTV}_{\rm sat}}=\frac{a_{\rm pl}^{1/2} F_{{\rm crit}} R_{\rm Hill} M_{\rm sat}}{M_{\rm sat}+M_{\rm pl}{(2{ GM}_{\star })}^{1/2}\Lambda }\,, \end{equation} \tag{ 27 }$$

where the geometric factor Λ appropriately amplifies or damps the TTVs in accordance with the eccentricity and orientation of the elliptical orbit of the planet about the parent star

$$\begin{eqnarray} &&\Lambda =\cos \left[ \tan ^{-1} \left(\frac{-e \cos (\omega)}{1+e\sin \omega } \right)\right]\left(\frac{2(1+e\sin \omega)}{1-e^2}-1\right)^{1/2}. \nonumber\\ \end{eqnarray} \tag{ 28 }$$

Substituting the various values discussed above, we find

$$\begin{equation} \delta _{{\rm TTV}_{\rm sat}}=0.15\,{\rm s}, \end{equation} \tag{ 29 }$$

which is several orders of magnitude too small to be of current interest. So while (somewhat surprisingly) it is distinctly possible that HAT-P-2b could currently harbor a massive satellite, any such satellite cannot be responsible for TTVs of the size that appear to be present.

4.3.2. TTVs Generated by an External Perturber with a Long Period

Perturbations from an as-yet undetected perturbing body present another potential explanation for the apparent TTVs. After the signature of HAT-P-2b's orbit has been removed from the existing RV observations, a secular acceleration that can be attributed to a distant companion remains. The constancy of the acceleration implies an orbital period for the companion of at least several years; for example, a body with mass, M_c ∼ 18 M_Jup, and semi-major axis, a_c = 10 AU, would suffice. For a configuration in which a_c ≫ a_pl (Holman & Murray 2005),

$$\begin{equation} \delta _{{\rm TTV}_{\rm perturber}} \sim \frac{45\pi }{16}\frac{M_c}{M_{\star }}P_{\rm pl} \alpha _e^3 \left(1-\sqrt{2}\alpha _{e}^{3/2} \right)^{-2}. \end{equation} \tag{ 30 }$$

where α_e = a_pl/(a_c(1 − e_c)). For a perturber with M_c ∼ 18 M_Jup, e_c = 0.3, and a_c = 10 AU, we find $\delta _{{\rm TTV}_{\rm perturber}}=0.1\,{\rm s}$. In addition, direct numerical integrations indicate that for external perturbing planets that are consistent with the observed secular acceleration, the expected TTVs are invariably very small, and furthermore, do not exhibit the rapid variation shown by the timing reversal observed between orbit 0 and orbit 83.

4.3.3. A Low-mass Resonant Perturber?

From our three transit observations we obtain planet–star radius ratios of 0.06821 ± 0.00075, 0.07041 ± 0.00060, and 0.0668 ± 0.0016 at 3.6, 4.5, and 8.0 μm, respectively. Our estimates of the planet/star radius ratio, R_p/R_⋆, are significantly smaller than the value presented in Pál et al. (2010), but are fairly well aligned with the values presented in Bakos et al. (2007a), Winn et al. (2007b), and Loeillet et al. (2008). Although Pál et al. (2010) incorporate the I band observations from Bakos et al. (2007a) they also include follow-up observations that utilize the z and Strömgren b+y bandpasses. While the I and z bands probe a similar wavelength range (∼0.8–0.9 μm), the Strömgren b and y bands probe a slightly shorter wavelength range (∼0.5 μm) where atmospheric scattering may become important.

The difference between our values of R_p/R_⋆ at 3.6, 4.5, and 8.0 μm could point to enhanced atmospheric opacity in the 4.5 μm bandpass due to CO, but our values of R_p/R_⋆ differ by less than 2σ. Given the large value of the average gravitational acceleration of HAT-P-2b (162 ± 27 m s⁻², Table 4), we would expect atmospheric scale heights to be small, which supports our wavelength independent values for the planetary radius. If we assume a value of $R_{\star }=1.64^{+0.09}_{-0.08} R_{\odot }$ (Pál et al. 2010), then we find an average radius value for HAT-P-2b of R_p = 1.106 ± 0.061R_J. The radius of the planet, given its mass, is in line with models of the thermal evolution of massive strongly irradiated planets. Planetary radii of 1.13 to 1.22 R_J for 8 M_J planets are expected for 1–10 Gyr ages (Fortney et al. 2007).

The estimates for the secondary eclipse depths at 3.6, 4.5, 5.8, and 8.0 μm presented here are the first secondary eclipse measurements for HAT-P-2b. We compare these results to the predictions of atmospheric models for this planet to better understand the atmospheric properties of HAT-P-2b. Figure 14 presents predicted day side emission as a function of wavelength from a one-dimensional radiative transfer model similar to those described in Burrows et al. (2008). These models assume a solar-metallicity atmosphere in local thermochemical equilibrium and include a parameterization for the redistribution of heat from the day side to the night side of the planet (P_n) and the possible presence of a high-altitude optical absorber (κ) It is important to note that these one-dimensional models assume instantaneous radiative equilibrium and therefore do not account for the expected phase lag between the incident flux and planetary temperatures for eccentric planets.

Zoom In Zoom Out Reset image size

Models that also include an additional high-altitude optical absorber to produce a dayside inversion best match our 3.6, 4.5, and 8.0 μm data points. We find some difficulty in matching the 5.8 μm data point. We note that even with our best attempts to remove the ramp from this data that the secondary eclipse depth appears to be slightly underestimated (Figure 8). Given the large uncertainty in the 5.8 μm secondary eclipse depth it is within 2.5σ of the flux ratio predicted by our models that include an inversion. Both the 4.5 and 8.0 μm secondary eclipse depths are more than 14σ larger than those predicted by atmospheric models without an inversion. Even if the planet were assumed to be ∼100 K warmer to account for the phase lag in planet-wide temperatures, models without an inversion would still underestimate the planetary flux in the 4.5 and 8.0 μm bands.

We calculate brightness temperatures in each band using a PHOENIX stellar atmosphere model for the star (Hauschildt et al. 1999). The 3.6 μm eclipse depths have an average value of 0.0996% ± 0.0072%, which corresponds to a brightness temperature of 2392 ± 84 K. We find an average 4.5 μm eclipse depth of 0.1031% ± 0.0061% and a corresponding brightness temperature of 2117 ± 65 K. Our eclipse depth at 5.8 μm corresponds to a brightness temperature of $1613^{+335}_{-150}$ K, while our eclipse depth at 8.0 μm corresponds to a brightness temperature of 2258 ± 106 K. Our secondary eclipse measurements in the 3.6 and 4.5 μm bandpasses agree at the 1σ level, and therefore do not provide any convincing evidence for variability in the planet's flux. Further observations of HAT-P-2b at 3.6 and 4.5 μm could place better constraints on possible orbit-to-orbit variability in the planet's thermal structure.

The overall shape and timing of the maximum and minimum of the planetary flux from our phase curves reveals a great deal about the atmospheric properties of HAT-P-2b. For planets on circular orbits, phase curve observations are generally related to day/night brightness contrasts on the planet. In the case of HAT-P-2b, the phase variations in the planetary flux are more indicative of thermal variations that result from the time-variable heating of the planet. In both the circular and eccentric cases the thermal phase amplitude and phase lag are determined by the radiative and advective timescales of the planet. By comparing the 3.6, 4.5, and 8.0 μm phase variations we can gain further insights into the thermal, wind, and possible chemical structure of HAT-P-2b's atmosphere.

We find that the peak of the observable planetary flux at 3.6 μm occurs 4.39 ± 0.28 hr after periapse passage with a peak value of 0.1139% ± 0.0089%, which corresponds to a brightness temperature of 2556 ± 100 K. The exact timing and level of the minimum planetary flux at 3.6 μm is much more uncertain. We find the planetary flux at 3.6 μm drops below observable levels for a period of ∼1 day 18.36$^{+0.18}_{-1.0}$ hr before transit. As such we can only put an upper limit of 1040 K on the minimum brightness temperature of HAT-P-2b at 3.6 μm.

The observed 4.5 μm HAT-P-2b phase curve exhibits a peak in the observable planetary flux 5.84 ± 0.39 hr after periapse passage with a value of $0.1162\%^{+0.0089\%}_{-0.0071\%}$ and corresponding brightness temperature of $2255^{+92}_{-74}$ K. We detect emission from HAT-P-2b over the entirety of its orbit at 4.5 μm with a minimum in the planet/star flux ratio of $0.0372\%^{+0.0086\%}_{-0.0096\%}$, which corresponds to a brightness temperature of $1345^{+119}_{-133}$ K. This minimum in the planetary flux occurs 6.71 ± 0.43 hr after transit. Roughly speaking, the shift of the minimum of the observed phase curve at 4.5 μm away from the region between apoapse and transit points to a minimum in the planetary temperature that is shifted west from the antistellar longitude and/or that the night side of the planet is still cooling even after the transit event. We would expect the planet to have a temperature minimum east of the anti-stellar point near the sunrise terminator assuming a super-rotating flow as shown for HD 189733b and HD 209458b in Showman et al. (2009). This dip in the planetary flux after transit is indeed puzzling and will be investigated further in a future study.

Our 8.0 μm observations cover only the portion of HAT-P-2b's orbit between transit and secondary eclipse, so we can only constrain the behavior of the 8.0 μm phase curve near the peak of the planetary flux. We find that the peak in the planetary flux at 8.0 μm occurs 4.64 ± 0.33 hr after the predicted time for periapse passage. The maximum of the planetary flux at 8.0 μm is 0.1888% ± 0.0072%, which corresponds to a brightness temperature of 2806 ± 79 K.

It is significant that we obtain a good fit to the data using Equation (12), which is similar to the functional form used to fit the light curves of HD 189733b and WASP-12b in Knutson et al. (2012) and Cowan et al. (2012a), respectively, and produce a longitudinal map of the planet's thermal variations. The "longitudinal" direction of thermal phase maps, ϕ, must be understood to mean "local stellar zenith angle," Φ. Thermal phase variations probe the diurnal cycle, T(Φ). A tidally locked planet has a one-to-one correspondence between local stellar zenith angle and longitude (e.g., ϕ = Φ), but phase maps can be made regardless of rotation rate (e.g., for Earth; Cowan et al. 2012b).

Nonetheless, there are two reasons why phase mapping should not work for eccentric planets: (1) the time-variable incident flux makes the brightness map time-variable (T(Φ, t)) and (2) the time-variable orbital angular velocity of the planet dictates that longitudinally sinusoidal variations in brightness on the planet would not correspond to sinusoidal variations in the light curve. Equation (12) is sinusoidal in the true anomaly (f) rather than in time, which implicitly accounts for (2).

The fact that Equation (12) fits the phase variations of HAT-P-2b implies that the diurnal brightness profile of the planet is constant throughout the orbit: dT(Φ)/dt = 0. This is entirely different from claiming that the longitudinal brightness profile of the planet is constant: since the planet is on an eccentric orbit, there is no fixed correspondence between longitude and stellar zenith angle. Rather, our data seem to indicate that the planet maintains a constant brightness as a function of stellar zenith angle, in marked disagreement with expectations for such an eccentric planet. This is almost certainly due to a geometric coincidence; however HAT-P-2b shows us its day side shortly after periapse. The day-side brightness of the planet likely changes throughout its orbit, but we are not privy to it.

4.5.1. Interpreting the Flux Maximum

4.5.2. Model of HAT-P-2b and Circular Analog

The top panels of Figure 15 shows how τ_rad and v_wind affect the 4.5 μm peak flux ratio and phase lag for HAT-P-2b. The dependencies are qualitatively similar for the other two wavebands. For comparison, we also show in bottom panels of Figure 15 the same dependencies for a hypothetical circular twin to HAT-P-2b with semi-major axis fixed at the actual planets periapse separation. There are a number of conclusions we can draw from these models.

• 1.  
  For circular planets, the radiative time scale and wind velocity are entirely degenerate: both the peak flux ratio and the phase lag depend solely on the product τ_radv_wind∝, the "advective efficiency" of Cowan & Agol (2011). Furthermore, the peak flux does not depend on the direction of zonal winds. The phase lag, on the other hand, would have the same amplitude but opposite sign for eastward versus westward zonal winds.
• 2.  
  The peak flux ratio for the eccentric planet HAT-P-2b depends approximately on the advective efficiency, , as for circular planets, but the dependence is no longer monotonic. The peak flux ratio exhibits a maximum for radiative times of ∼6 hr and wind velocities of 2 km s⁻¹. This is because eastward advection of heat brings the hot spot into view shortly after periapse, as discussed in Cowan & Agol (2011). There is a local minimum in peak flux for τ_rad = 0, the radiative equilibrium case. The global minimum, however, still occurs in the limit of long radiative times and rapid winds, which results in a planet with zonally uniform, time-constant temperature, as in the circular case.
• 3.  
  Comparing the top-right and bottom-right panels of Figure 15 shows that zonal winds have qualitatively the same effect, regardless of orbital eccentricity: eastward winds make the peak flux occur early, while westward winds cause a delay in the peak flux. The major difference between the two geometries is the phase-lag expected in the absence of winds, the null hypothesis. For a circular, tidally locked planet there is no phase lag in this limit, whereas the eccentric planet exhibits a large positive phase lag in the absence of winds. As a result of this different zero-point, the amplitude of phase lag from periapse decreases nearly monotonically for increasing in the eccentric model, exactly the opposite behavior from a circular planet. Such counterintuitive behavior is precisely why it is more useful to compare phase lags to the windless scenario rather than quote absolute phase lags from periapse (Cowan & Agol 2011).

Zoom In Zoom Out Reset image size

4.5.3. Constraints on Model Parameters

In Figure 16 we show exclusion diagrams in the v_wind versus τ_rad plane for each waveband. Since we use a one-layer model, this can roughly be thought of as constraining the properties of the photospheres at each waveband, with the understanding that the mid-IR vertical contribution functions span approximately a factor of 10 in pressure (e.g., Knutson et al. 2009a; Showman et al. 2009). Blue lines in Figure 16 show the 1σ and 2σ confidence intervals from the peak flux value. Red lines show the same for the phase offset. The gray scale shows the combined confidence intervals.

Zoom In Zoom Out Reset image size

At 3.6 and 4.5 μm, the peak fluxes are consistent with a broad range of eastward zonal wind scenarios. Only very high values of (top-right of plot) and most westward winds (bottom of plot) are excluded. At 8 μm, only the highest peak flux values are favored, but even these are a very poor fit. If we adopt a Bond albedo of zero, the 8 μm peak flux still disagrees with any models by >5σ. We therefore conclude that the high flux at 8 μm is not due solely to atmospheric dynamics, but also to chemistry and the planet's vertical temperature profile. This waveband falls within a water vapor absorption feature and is therefore expected to originate higher up in the atmosphere than either the 3.6 and 4.5 μm flux. The high flux in the 8.0 μm channel therefore suggests a temperature inversion, which is also supported by our measured secondary eclipse depths in Section 4.4.

The constraints from the phase lag of peak flux (red lines in Figure 16) are stronger, albeit still degenerate. The data favor a narrow range of advective efficiencies (τ_radv_wind). The bottom panel of Figure 16 gives the impression that the peak flux and phase lag at 8 μm combine to give a nice constraint on the model parameters, but the peak flux constraints should be taken with a grain of salt, as described above. It is important to note the decaying exponential-like trend of the preferred values of the zonal wind speeds with τ_rad in Figure 16 that pairs short values for τ_rad with higher values for v_wind and longer values for τ_rad with lower values for v_wind. Three-dimensional models that consistently couple radiative and advective processes will help to further constrain the relevant timescales in HAT-P-2b's atmosphere over the full range of its orbit.