Orbital Phase Variations of the Eccentric Giant Planet HATP2b
Abstract
We present the first secondary eclipse and phase curve observations for the highly eccentric hot Jupiter HATP2b in the 3.6, 4.5, 5.8, and 8.0 m bands of the Spitzer Space Telescope. The 3.6 and 4.5 m data sets span an entire orbital period of HATP2b ( d), making them the longest continuous phase curve observations obtained to date and the first fullorbit observations of a planet with an eccentricity exceeding 0.2. We present an improved nonparametric method for removing the intrapixel sensitivity variations in Spitzer data at 3.6 and 4.5 m that robustly maps positiondependent flux variations. We find that the peak in planetary flux occurs at 4.390.28, 5.840.39, and 4.680.37 hours after periapse passage with corresponding maxima in the planet/star flux ratio of , , and in the 3.6, 4.5, and 8.0 m bands respectively. Our measured secondary eclipse depths of , , , and in the 3.6, 4.5, 5.8, and 8.0 m bands respectively indicate that the planet cools significantly from its peak temperature before we measure the dayside flux during secondary eclipse. We compare our measured secondary eclipse depths to the predictions from a onedimensional radiative transfer model, which suggests the possible presence of a transient day side inversion in HATP2b’s atmosphere near periapse. We also derive improved estimates for the system parameters, including its mass, radius, and orbital ephemeris. Our simultaneous fit to the transit, secondary eclipse, and radial velocity data allows us to determine the eccentricity () and argument of periapse () of HATP2b’s orbit with a greater precision than has been achieved for any other eccentric extrasolar planet. We also find evidence for a longterm linear trend in the radial velocity data. This trend suggests the presence of another substellar companion in the HATP2 system, which could have caused HATP2b to migrate inward to its presentday orbit via the Kozai mechanism.
Subject headings:
planets and satellites: general, planets and satellites: individual: HATP2b, techniques: photometric, methods: numerical, atmospheric effects1. Introduction
The supermassive () Jupiter sized () planet HATP2b (aka HD 147506b) was first discovered from transit observations by Bakos et al. (2007b) using the HATNet (Bakos et al., 2002, 2004) network of groundbased telescopes. Followup radial velocity measurements of the HATP2 system revealed that the orbit of HATP2b is highly eccentric ( Bakos et al., 2007b; Loeillet et al., 2008). Only a handful of transiting exoplanets have been shown to posses eccentricities in excess of that of Mercury (), which has made HATP2b an interesting target for many theoretical studies concerning the evolution of the HATP2 system (Jackson et al., 2008; Fabrycky, 2008; Matsumura et al., 2008; Baraffe et al., 2008). Because of its mass, HATP2b represents a class of exoplanets that provides an important link between extrasolar giant planets and lowmass brown dwarfs (Baraffe et al., 2008). Observational constraints on the basic atmospheric properties of HATP2b will provide an important probe into the structure and evolution of not only HATP2b, but an entire class of massive extrasolar planets.
Atmospheric circulation models for planets on eccentric orbits show significant variations in atmospheric temperature and wind speeds that provide an important probe into atmospheric radiative and dynamical timescales (Langton & Laughlin, 2008; Lewis et al., 2010; Cowan & Agol, 2011; Kataria et al., 2012). The incident flux on HATP2b from its stellar host at periapse is ten times that at apoapse, which should cause large variations in atmospheric temperature, wind speeds, and chemistry. Heating and cooling rates for HATP2b can be constrained by measuring planetary brightness as a function of time. Such observations are analogous to previous phase curve observations of HD 189733b using the Spitzer Space Telescope, which provided the first clear observational evidence for atmospheric circulation in an exoplanet atmosphere (Knutson et al., 2007, 2009b, 2012).
HATP2b’s large orbital eccentricity makes it an ideal target for investigating hot Jupiter migration mechanisms. Gas giant planets such as HATP2b are expected to form beyond the ice line in their protoplanetary disk far from their stellar hosts. HATP2b currently resides at a semimajor axis of 0.07 AU from its host star, indicating that it must have migrated inward via a process such as gas disk migration (Lin et al., 1996), planetplanet scattering (Rasio & Ford, 1996), secular interactions (Wu & Lithwick, 2011), or Kozai migration (Weidenschilling & Marzari, 1996; Naoz et al., 2011). HATP2b’s closein and highly eccentric orbit favors one of the latter three mechanisms since disk migration tends to damp out orbital eccentricities. For Kozai migration, the presence of a third body with at least as much mass as HATP2b is needed. In this study we present six years of radial velocity measurements for this system, which allow us to search for the presence of a massive third body in the HATP2 system.
Here we present our analysis of the 3.6 and 4.5 m fullorbit phase curves of the HATP2 system, which include two secondary eclipses and one transit at each wavelength. These fullorbit phase curves represent the longest continuous phase observations ever obtained by the Spitzer Space Telescope. The orbital period of HATP2b (5.6334729 d) is more than 2.5 times that of other exoplanets with published fullorbit phase curve observations: WASP12b (Cowan et al., 2012a) and HD 189733b (Knutson et al., 2012). Additionally, we present an analysis of previous partial orbit phase curve and secondary eclipse Spitzer observations at 8.0 and 5.8 m respectively. We use these observations to characterize the changes in the planet’s emission spectrum as a function of orbital phase and to probe the atmospheric chemistry and circulation regime of HATP2b. The following sections describe our observations and data reduction methods (§2) and the results of our analysis (§3). Section 4 discusses the results from our analysis of the Spitzer data and compares them to predictions from atmospheric models for HATP2b. Additionally, we discuss trends in our radial velocity data that indicate the presence of an additional body in HATP2 system in Section 4. Section 5 overviews the main conclusions from our analysis and presents ideas for future work.
2. Observations
We analyze nearly 300 hours of observation of HATP2 at 3.6 m and 4.5 m obtained during the postcryogenic mission of the Spitzer Space Telescope (Werner et al., 2004) using the IRAC instrument (Fazio et al., 2004) in subarray mode. The observation periods were UT 2010 March 28 to UT 2010 April 03 and UT 2011 July 09 to UT 2011 July 15 for the 3.6 and 4.5 m bandpasses respectively. Both observations cover a period just over 149 hours with two approximately 2hour breaks for data downlink, corresponding to 1.2 million images in each bandpass. Each observation begins a few hours before the secondary eclipse of the planet, continues through planetary transit, and ends a few hours after the subsequent planetary secondary eclipse.
We also analyze observations of the HATP2 system at 5.8 m and 8.0 m obtained during the cryogenic phase of the Spitzer Space Telescope mission. The 5.8 m observations cover a single secondary eclipse of the planet that occurs on UT 2009 March 16. The 8.0 m observations cover a portion of the planet’s orbit that includes transit, periapse passage, and secondary eclipse on UT 2007 September 1011. The 5.8 m observations were obtained using subarray mode, while the 8.0 m observations were obtained using the full IRAC array. We obtain our 3.6, 4.5, 5.8, and 8.0 m photometry from the Basic Calibrated (BCD) files produced by version 18.18.0 of the Spitzer analysis pipeline.
In subarray mode, 3232 pixel images are stored in sets of 64 as a single FITS file with a single header. We calculate the BJD_UTC at midexposure for each image from time stamp stored in the MBJD_OBS keyword of the FITS header, which corresponds to the start of the first image in each set of 64. We assume uniform spacing of the images over the time period defined by the AINTBEG and ATIMMEEND header keywords. The image spacing is roughly equal to the 0.4 s exposure time selected for the 3.6, 4.5, and 5.8 m observations of HATP2 (). For the 8.0 m full array observations, the MBJD_OBS, AINTBEG, and ATIMMEEND header keywords are used to calculate the BJD_UTC at midexposure for each image, which are spaced by roughly the 12.0 s exposure time. In order to convert from UTC to TT timing standards as suggested by Eastman et al. (2010), 66.184 s would be added to our BJD_UTC values. We report all of our timing measurements using BJD_UTC for consistency with other studies.
2.1. 3.6 and 4.5 m Photometry
We determine the background level in each image from the region outside of a ten pixel radius from the central pixel. This minimizes contributions from the wings of the star’s point spread function 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 fluxweighted 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 fluxweighted centroiding estimates and find that for the HATP2 data fluxweighted 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 hours. We find that fluxweighted 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 fluxweighted 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 timevarying aperture sizes based on the noise pixel parameter described in the Appendix We find that we obtain the lowest standard deviation of the residuals from our bestfit solution at 3.6 m using a variable aperture size given by the squareroot 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 ‘sawtooth’ 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 the Appendix.
2.2. 5.8 m Photometry
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 (§2.1). The background flux is 0.6% of the total flux in the science aperture. We tested both fluxweighted 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 fluxweighted centroids. This preference for the Gaussian centroiding method is likely the result of the shorttime 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 the Appendix, but find that including intrapixel sensitivity corrections actually degraded the precision of the final solution.
2.3. 8.0 m Photometry
Unlike the 3.6, 4.5, and 5.8 m data, the 8.0 m data were observed in the fullarray (256256 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 2D Gaussian to a 3x3pixel median filtered version of the frame; we find that this method produces slightly better photometric precision than using a centeroflight 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 the Appendix, but find that including intrapixel sensitivity corrections actually degraded the precision of the final solution.
2.4. Flux Ramp Correction
Our data exhibit a ramplike 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, 2009a); 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 wellunderstood, 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 chargetrapping 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):
(1) 
where is the measured flux, is the stellar flux incident on the array, is the time from the start of the observation in days, and 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 HATP2b 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. (2011). As seen in Figure 4, the 8.0 m data exhibit the well know asymptotic increase in flux (e.g. Knutson et al., 2007, 2009a; 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
(2) 
where is the degrees of freedom in the fit and is the total number of points in the fit (Liddle, 2007). This allows us to determine if we are ‘overfitting’ 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 () 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 than this amount ether gives a higher BIC or reduces the amount of out of eclipse data to less than the eclipse duration.
2.5. Transit and Eclipse Fits
We model our transit and eclipse events using the equations of Mandel & Agol (2002) modified to account for the orbital eccentricity, , and argument of periapse, , of the HATP2 system. For an eccentric system the normalized separation of the planet and star centers, , is given by
(3) 
where is the radial planetstar distance as a function of time, is the stellar radius, is the orbital inclination of the planet, and is the true anomaly as a function of time. The term in Equation 3 is calculated as
(4) 
The true anomaly angle () that appears in both Equations 3 and 4 is determined using Kepler’s equation (Murray & Dermott, 1999) and is a function of , , and the orbital period of the planet, . Because of the degeneracies that exist between , and in determining we elect to not use as a free parameter in our fits. Instead we fix to the value reported in Pál et al. (2010), 5.6334729 d. We further minimize correlations in our solutions for and by solving for the Lagrangian orbital elements and .
In addition to the parameters that define the orbit of HATP2b, we solve for the fractional planetary radius, . Because for HATP2b (Pál et al., 2010; Bakos et al., 2007b), there exists a strong correlation between and in the transit solution (Winn et al., 2007b; Pál, 2008). We instead solve for the parameters and as suggested by Bakos et al. (2007a) and Pál et al. (2010) where
(5) 
and
(6) 
For the transit portion of the light curves we use four parameter nonlinear limbdarkening coefficients for each bandpass calculated by Sing (2010), where we assume a stellar atmosphere with K, , 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 (§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 datasets, 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 , , , , and to the average values from the 3.6, 4.5, and 8.0 m data sets.
2.6. Phase Curve Fits
The functional form of the phase curve for a planet on an eccentric orbit is not well defined. Unlike closein, tidally locked planets on circular orbits, eccentric planets experience timevariable heating and nonsynchronous rotation rates. Previous studies by Langton & Laughlin (2008) and Cowan & Agol (2011) have investigated theoretical light curves for planets on eccentric orbits using twodimensional hydrodynamic simulations and semianalytic model atmospheres respectively. We also developed a threedimensional atmospheric model for HATP2b 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 threedimensional 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, . 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 HATP2b is a simple with a phase lag:
(7) 
where is the true anomaly and are free parameters in the fit. The is a simplified form of the denominator of Equation 4 for . We also test a simpler form of Equation 7 given by
(8) 
where 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 HATP2b. This is not surprising given that we expect the flux from the planet to vary as 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, , given by
(9) 
where
(10) 
in the symmetric case and
(11) 
in the asymmetric case. In Equations 9, 10, and 11 are free parameters in the fit. The Lorentzian functional form for the phase curve is especially useful since the parameter gives the offset between the time of periapse passage and the peak of the planet’s flux and the 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, . In this case
(12)  
where are free parameters in the fit and such that transit occurs at and secondary eclipse occurs at .
In addition to the functional forms for the phase curve presented above, we also test a flat phase curve, , to make sure that we are not overfitting the data with the phase curve parameters. In all cases we tie the 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 LevenbergMarguardt 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=1,282,575; BIC=1,459,141; BIC=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 fixed at zero (BIC=1,278,521). This is not surprising since basic trigonometric identities make Equation 12 equivalent to Equation 8 if 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 and 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 ten hours, we find no statistical difference between solutions with and without parameterizations for planetary phase variations.
2.7. Stellar Variability
HATP2 is a rapidly rotating F star (=20.7 km s). Measurements of the Ca II H & K lines of HATP2 by Knutson et al. (2010) do not detect any significant emission in the line cores, which would suggest a chromospherically quiet stellar host for HATP2b. 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 & K lines. Groundbased monitoring of HATP2 in the Strömgren and 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 HATP2 or longerterm 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., 2007a) based on the lineofsight stellar rotation velocity () and assuming , which is roughly 0.6 times the orbital period of HATP2b. 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
(13) 
where is measured from the predicted center of transit in each observation and is a free parameter in the fit. The second model we test has the form
(14) 
where is measured from the predicted center of transit in each observation and 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 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 ‘sinecurve’ 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 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 HATP2, the amplitude of the predicted stellar variations in the midinfrared seem spurious. The amplitude of the predicted stellar flux variations at 3.6 and 4.5 m are on the order of , which is comparable to our upper limit on stellar variability at visible wavelengths (). We would expect star spots to display much larger variations in the visible than at midinfrared wavelengths. The amplitude of these stellar variations are also similar to the amplitude of our predicted phase variations. We therefore conclude that our bestfit 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.
2.8. Radial Velocity Measurements
Since the discovery of HATP2b by Bakos et al. (2007b), the California Planet Search (CPS) team has continued to obtain regular radial velocity (RV) measurements of this system using the HIRES instrument (Vogt et al., 1994) on Keck. We used the CPS pipeline (see, e.g. ?) to measure precise RVs from the highresolution spectra of HATP2 using a superposed molecular iodine spectrum as a Doppler reference and point spread function calibrant. Previous studies of this system (Bakos et al., 2007b; Winn et al., 2007a; Pál et al., 2010) have included HIRES radial velocity measurements through May 2008 (BJD = 2454603.932112). Here we present 16 additional RV measurements from the HIRES instrument on Keck for a total of 71 data points spanning 6 years (Table 1). We exclude measurements obtained during transit (Winn et al., 2007a) from our fits in order to avoid measurements affected by the RossiterMcLaughlin effect.
BJD  2400000  RV  Uncertainty 

(days)  (m s)  (m s) 
53981.777500  19.17  7.46 
53982.871700  310.87  7.64 
53983.814865  538.73  7.81 
53984.894980  855.62  7.94 
54023.691519  698.64  7.94 
54186.998252  696.78  7.72 
54187.104158  684.36  6.84 
54187.159878  717.17  6.59 
54188.016885  757.70  7.00 
54188.159622  774.61  6.80 
54189.010378  651.86  6.59 
54189.088911  635.16  6.48 
54189.157721  619.70  7.22 
54216.959395  722.59  8.04 
54257.756431  27.91  5.79 
54257.758677  35.49  6.17 
54257.760702  24.33  6.18 
54257.794116  11.70  5.48 
54257.796779  19.01  5.02 
54257.799452  20.94  5.10 
54257.802148  23.48  4.74 
54257.804926  21.89  5.15 
54257.807634  35.65  4.91 
54257.810342  24.40  5.08 
54257.813143  40.82  5.22 
54257.815817  37.08  5.45 
54257.818490  23.98  5.21 
54258.024146  313.00  4.87 
54258.027167  310.40  4.89 
54258.030292  311.56  4.37 
54258.033278  326.04  4.58 
54258.042630  314.19  5.10 
54258.045488  336.27  4.92 
54258.048393  342.10  5.07 
54258.051483  351.92  4.88 
54258.054828  356.80  4.94 
54258.058161  356.72  4.65 
54258.061472  360.49  4.85 
54258.064666  374.47  4.73 
54258.099110  432.29  5.27 
54258.102836  433.53  4.97 
54258.106679  446.91  4.74 
54279.876893  371.51  8.26 
54285.823854  137.15  5.91 
54294.878702  728.72  6.61 
54304.864982  588.92  6.07 
54305.870122  737.55  6.23 
54306.865216  731.39  7.95 
54307.912379  456.62  6.39 
54335.812619  549.60  6.64 
54546.098175  684.07  7.79 
54547.115700  536.15  7.19 
54549.050468  754.51  7.26 
54602.916550  271.30  6.47 
54603.932112  662.48  5.74 
54839.166895  369.59  8.09 
55015.871081  679.68  8.89 
55350.945695  513.11  8.45 
55465.742817  606.08  6.99 
55703.872995  623.50  7.70 
55704.836273  377.02  8.02 
55705.853346  535.77  8.12 
55706.831645  257.19  7.57 
55707.844353  504.50  7.68 
55808.759715  323.97  8.49 
55850.695496  557.65  7.37 
55932.161110  241.35  9.48 
55945.128011  595.28  8.57 
55992.018498  401.30  10.08 
56147.753065  553.56  7.95 
56149.738596  381.71  8.37 
We simultaneously fit for the RV semiamplitude (), zeropoint (), and long term velocity drift () 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 (), but find that our derived values for were consistent with zero. As discussed in Winn (2010), the transit to secondary eclipse timing strongly constrains the term in our fits, but the 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., 2007b; Winn et al., 2007a). We therefore estimate a stellar jitter term () as described in Section 3.
3. Results
Parameter  3.6 m  4.5 m  5.8 m  8.0 m 

Transit Parameters  
0.068210.00075  0.070410.00060  0.06933 
0.06680.0016  
0.122  0.3450.042  ….  0.238  
()  87.37  84.910.47  85.97 
86.0 
9.530.35  8.280.24  8.70 
8.83  
(BJD2400000) 
55288.849880.00060  55756.425200.00067  …  54353.69110.0012 
(d)  12.2210.058  12.2860.057  …  12.210.11 
(d) 
0.17700.0011  0.18130.0013  …  0.17890.0023 
(d) 
0.01280.0010  0.01770.0012  …  0.01440.0021 
Secondary Eclipse Parameters  
0.191770.00029  0.193100.00025  …  0.192530.00048  
(d) 
0.15500.0027  0.16510.0023  0.16100.0043  
(d) 
0.010900.00075  0.014440.00083  0.01210.0016  
1 Eclipse Depth  0.080%0.012%  0.1009%0.0084%  …  … 
(BJD2400000) 
55284.29670.0014  55751.87950.0011  …  … 
2 Eclipse Depth  0.0993%0.0090%  0.1057%0.0090%  0.071%  0.1392%0.0095% 
(BJD2400000) 
55289.93020.0014  55757.51300.0011  54906.8561  54354.77570.0022 
Orbital & RV Parameters  
0.505390.00057  0.503010.00051  …  0.504190.00088  
0.07410.0072  0.07290.0054  …  0.0685 0.0059  
0.510810.00092  0.508290.00068  0.50910 
0.508850.00097  
()  188.340.80  188.250.62  188.09 
187.740.66 
(BJD2400000)  55289.47340.0079  55757.05194 0.0061  …  54354.31090.0065 
(m s)  927.05.9  923.05.8  …  923.16.0 
(m s)  247.43.7  248.03.6  …  248.03.6 
(m s d)  0.08900.0050  0.08810.0049  …  0.08860.0058 
Phase Curve Parameters  
Functional Form  Eq. (12)  Eq. (12)  Flat  Eq. (7) 
0.0379%0.0015%  0 (fixed)  …  0.0555%0.0032  
0.0422%0.0025%  0.0293%0.0042%  …  41.82.7  
0 (fixed)  0 (fixed)  …  
0 (fixed)  0.0163% 0.0035%  …  …  
Amplitude  0.114%0.010%  0.079%0.013%  …  … 
Minimum Flux  0.00014%  0.0372%  …  … 
Minimum Flux Offset (h) 
18.36  6.710.43  …  … 
Maximum Flux  0.1139%0.0089%  0.1162%  …  0.1888%0.0072% 
Maximum Flux Offset (h) 
4.390.28  5.840.39  …  4.680.37 
Ramp Parameters  
Functional Form  Eq. (1)  None  Eq. (1)  Eq. (1) 
0.00134  …  +0.00683  0.00537  
0.078  …  0.0950.015  0.0195  
0 (fixed)  …  …  0.01765  
0 (fixed)  …  …  0.2812  
Noise Parameters  
0.00422090.0000024  0.0057064 0.0000032  0.0153800.000034  0.0021640.000015  
(m s)  26.02.1  25.72.2  …  25.52.1 
We perform a simultaneous fit and calculate uncertainties for the relevant transit, secondary eclipse, phase curve, radial velocity, flux ramp, and intrapixel sensitivity correction parameters in our data sets using a Markov Chain Monte Carlo (MCMC) method (Ford, 2005). For the 3.6, 4.5, and 8.0 m data sets fit parameters include , , , , , transit time, the secondary eclipse depth(s), the phase function coefficients , a photometric noise term (), and RV parameters , , , . For the 3.6 and 8.0 m data sets we additionally fit for the ramp correction coefficients . The free parameters in the fit to the 5.8 m data set are the eclipse time, eclipse depth, ramp correction coefficients , and . The value of the stellar flux, , is inherently accounted for with our intrapixel sensitivity correction method described in the Appendix. Because we do not apply intrapixel sensitivity corrections to our 5.8 and 8.0 m data sets we additionally fit for in those cases.
The only parameter held fixed in our analysis is the orbital period (), the value for which we take from Pál et al. (2010). All other orbital, planetary, phase, and ramp parameters are allowed to vary freely. Initial attempts at fitting just the 3.6 and 4.5 m transits simultaneously produce a value for within 1 of the Pál et al. (2010) value, and we are therefore confident that adopting the Pál et al. (2010) value for has not introduced any significant errors into our analysis.
We initially attempted to fit for the wavelength independent parameters , , , , and simultaneously for the RV, 3.6, 4.5, and 8.0 m data sets. However, given the size of our data sets (over 2.5 million data points) and the time required to create our ‘pixelmaps’ for the 3.6 and 4.5 m data (see Appendix),we found it computationally infeasible. It is possible that further improvements to our analysis code, including parallelization, could make the problem more computationally tractable. Such improvements are left for future iterations of our analysis methods.
We plot the normalized time series for the 3.6, 4.5, and 8.0 m data sets after the bestfit intrapixel sensitivity variations and ramp corrections have been removed in Figure 5. The regions near secondary eclipse and transit for each best fit solution are presented in Figures 6, 7, 8, and 9 for the 3.6, 4.5, 5.8, and 8.0 m data sets respectively. We also present our best fit solution to the RV data in Figure 10.
We use a total of five independent chains with steps per chain in our MCMC analysis. Each chain is initialized at a position in parameter space determined by randomly perturbing the bestfit parameters from a LevenbergMarquardt nonlinear least squares fit to the data (Markwardt, 2009). Instead of using a standard minimization scheme, we instead opt to maximize the of the likelihood () given by
(15) 
where is the relevant error term. This allows us to simultaneously solve for the noise terms and with our other fit parameters (see Carter & Winn, 2009, for further discussion of the maximum likelihood method as applied to exoplanet transit observations). After each chain has reached steps, we find the the point where first surpasses the median value and discard all steps up to that point. We then combine the results from our independent chains and find the range about a median value that contains 68% of the values for a given parameter. We set our bestfit parameters equal to this median value and use this distribution to initially determine the 1 uncertainties in each of our parameters. For most parameters our error distribution was close to being symmetric about the median value. In the cases where the distribution was significantly asymmetric, we have noted both the positive and negative uncertainties in the parameter value.
We find that there is (‘red’) correlated noise in our data even after the best fit intrapixel sensitivity variations have been removed. To account for correlated noise in our data we first employed the waveletbased MCMC method described in Carter & Winn (2009). However, we found that our correlated noise could not be treated as a stationary noise parameter. Often we found that fits to the transit and secondary eclipse portions of our phase curves were degraded to introduce ‘red’ noise consistent with other portions of the light curve. As a result, we instead opt to use the ‘residualpermutation’ or ‘prayerbead’ method to estimate the errors in our parameters in the presence of correlated noise (see, for example Jenkins et al., 2002; Southworth, 2008; Bean et al., 2008; Winn et al., 2008). The ‘prayerbead’ errors are typically 1.53 larger than the errors determined from our MCMC analysis. The largest increases in the uncertainty using the ‘prayerbead’ method were for the planetstar ratio and secondary eclipse depths. In some cases the ‘prayerbead’ errors are slightly smaller () than the errors from the MCMC analysis. In those cases we report the larger errors from the MCMC analysis. Table 2 presents the bestfit parameters for our data sets and their 1 error bars. We find that our photometric errors, , are 1.05, 1.11, 1.11, and 1.15 times higher than the predicted photon noise limit at 3.6, 4.5, 5.8, and 8.0 m respectively.
4. Discussion
In the following sections we discuss the implications of these results for our understanding of the HATP2 system and the atmospheric properties of HATP2b. We compare our results with the previous studies of the HATP2 system from Bakos et al. (2007b), Winn et al. (2007a), Loeillet et al. (2008), and Pál et al. (2010), which were limited to ground based transit and radial velocity data. We also compare our results to predictions from onedimensional radiative transfer and semianalytic models of HATP2b’s atmosphere.
4.1. Orbital and RV Parameters
Parameter  Bakos et al. (2007b)  Winn et al. (2007a)  Loeillet et al. (2008)  Pál et al. (2010) 

Transit Parameters  
0.06840.0009  0.06810.0036 
0.06891  0.072270.00061  
()  84.6  86.8  90.0  86.72 
9.77  9.900.39 
10.28  8.99  
(BJD2400000)  54212.85630.0007  54212.85650.0006  …  54387.493750.00074 
(d)  0.1770.002  …  …  0.17870.0013 
(d)  0.0120.002  …  …  0.0141 
Secondary Eclipse Parameters  
0.18470.0055 
0.19690.0040 
0.18960.0016 
0.18680.0019  
(d)  …  …  …  0.16500.0034 
(BJD2400000)  …  …  …  54388.5460.011 
Orbital & RV Parameters  
(d)  5.633410.00013  5.63341 (fixed)  5.63341 (fixed)  5.63347290.0000061 
…  …  …  0.51520.0036  
…  …  …  0.04410.0084  
0.5200.010  0.5010.007  0.5163  0.51710.0033  
()  179.33.6  187.41.6  189.92  185.220.95 
(BJD2400000)  54213.3690.041  …  54213.4798  … 
(m s)  101138  88357 
966.98.3  983.917.2 
Planetary Parameters  
()  9.040.50  8.040.40  8.62  9.090.24 
()  0.982  0.980.04  0.951  1.157 
(g cm)  11.9  10.600.55 
12.5  7.291.12 
(m s)  227  20720 
237  16817 
a (AU)  0.06770.0014  0.06810.0014 
0.0677  0.068780.00068 
Noise Parameters  
(m s)  60  31  17  … 
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 HATP2 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 groundbased transit and radial velocity data. By the inclusion of secondary eclipse in our data we were able to improve the estimate of the orbital eccentricity of HATP2b by an order of magnitude over the value presented in Pál et al. (2010).
From our orbital and RV parameters we can estimate the mass of HATP2b using
(16) 
where the orbital period () and stellar radius () are assumed to be the values presented in Pál et al. (2010), d and respectively. Values for the RV semiamplitude (), eccentricity (), inclination (), and normalized semimajor axis () are taken as the errorweighted 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 HATP2b to be , which is within 1 of the previous estimates of for HATP2b (Table 3).
Parameter  Value 

0.069330.00045  
(BJD)  2455288.849230.00037 
8.700.19  
()  85.97 
0.509100.00048  
188.090.39  
0.192530.00018  
(BJD)  2455289.932110.00066 
(BJD)  2455289.47210.0038 
()  8.000.97 
()  1.1060.061 
(g cm)  7.31.6 
(m s)  16227 
(AU)  0.06630.0039 
4.2. Linear RV Trend
Our RV data spans a period of nearly six years, which allows us to test for long term trends in the RV signal. We find a nonzero value for the linear term in our RV fit (), which indicates the presence of a second body orbiting in the system. If we assume that this second companion to HATP2 (‘c’) is on a circular orbit and is much less massive than its host star, we can relate to the mass and orbital semimajor axis of ‘c’ by the expression presented in Winn et al. (2009)
(17) 
Figure 11 shows the range of values for and that are allowed for the HATP2 system given that . 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 and AU based on an orbital period for ‘c’ of 24 years. We further employ adaptive optics imaging to search for ‘c’ and set an upper limit on the values of and .
Adaptive Optics Imaging
In an attempt to directly image the body responsible for causing the linear RV trend, we observed HATP2 on May 29, 2012 using NIRC2 (PI: Keith Matthews) and the Keck II adaptive optics (AO) system at Mauna Kea (Wizinowich et al., 2000). Our observations consist of dithered images taken with the K’ (m) filter. Using the narrow camera setting, which provides fine spatial sampling of the NIRC2 pointspreadfunction (10 mas ), we acquired 9 frames each having 13.2 seconds of onsource integration time. The seeing was estimated to be 0.4 at visible wavelengths at the time of observations using AO wavefront sensor telemetry. We note that significant windshake degraded the quality of correction for some images but by only a marginal amount.
The data were processed using standard techniques to correct for hot pixels, remove background radiation from the sky and instrument optics, flatfield the array, and align and coadd individual frames. Figure 12 shows the final reduced AO image along with the corresponding () contrast levels achieved as a function of angular separation. No companions were detected in raw or processed frames.
We can use the limits from a nondetection to rule out the presence of companions as a function of and . Using HATP2’s parallax distance of and estimated age of Gyr as determined by Pál et al. (2010) through a combined isochrone, light curve, and spectroscopic analysis, we find that our diffractionlimited observations are sensitive to companions on the hydrogenfusing boundary for separations beyond . Interior to this region, the combination of imaging and RV data eliminates most lowmass stars, though latetype Mdwarf tertiaries located at AU could cause the longterm Doppler drift yet simultaneously evade direct detection (Figure 11)
Orbital Evolution
The likely presence of an M/L/T/Y dwarf at an orbital distance, , of 10 to 40 AU from HATP2b lends credence to the possibility that HATP2b owes is current orbit to a history of Kozai cycling (Kozai, 1962) combined with longrunning 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 HATP2b and indeed, longdistance disk migration for HATP2b seems quite problematic. The protostellar disk itself would have needed to be exceptionally massive and long lived, and an adhoc mechanism must be invoked to explain the large observed orbital eccentricity of HATP2b. Furthermore, HATP2b has planetary and orbital parameters that place it significantly outside the welldelineated population of ‘conventional’ hot Jupiters with days, , and .
In the context of the KCTF process, HATP2b 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, (assuming , and an initially circular orbit for ). 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
(18) 
With plausible values of years, , and , Kyr. By contrast, the general relativistic precession rate for HATP2b is currently
(19) 
which generates a full circulation of the apsidal line in 20 Kyr (). Because , the magnitude of the Kozai cycling is strongly suppressed at present.
To good approximation, tidal friction is currently the dominant contributor to the orbital evolution of HATP2b, and it is plausible that the current state of the system has undergone dissipative evolution from an earlier epoch where . If we assume that the tidal evolution has roughly conserved HATP2b’s periastron distance, AU, then in order for implies that HATP2b has evolved to its current orbit from an orbit with an initial period, , given by
(20) 
or
(21) 
Additionally, the radial velocityderived constraint on the unseen companion ‘c’ indicates that , which allows us to simply write
(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 HATP2b to be .
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 spinorbit misalignment angle, , suggested that HATP2b’s orbit lies in the stellar equatorial plane. Recently, however, an reassessment by Albrecht et al. (2012) reports , indicating a modest projected misalignment. In addition, HATP2b’s parent star is almost precisely on the K boundary at which stars empirically appear to be able to maintain misalignment over multiGyr time scales (Albrecht et al., 2012).
4.3. Transit Timing Variations
We calculate ephemeris for the HATP2 system given the center of transit and eclipse times presented in Table 2 as:
(23) 
where is the predicted transit, eclipse, or periapse time and is the number of orbits that have elapses since . From our transit data we calculate BJD and d, which is consistent with the orbital period for HATP2b presented in Pál et al. (2010). We calculate BJD and d from our secondary eclipse data and BJD and 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 planetstar 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 HATP2 system of nearly two minutes. Some of the spread in the measured vs. predicted transit times could be accounted for by inconsistent definitions of transit center.
The separate visits to HATP2b for the 3.6, 4.5, 5.8, and 8.0 m observations allow constraints to be placed on possible transit timing variations (TTV) for HATP2b. Previous transit timing measurements for this system (Bakos et al., 2007b; Winn et al., 2007a; Pál et al., 2010) had relatively low timing resolution (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 anticorrelated 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 asyet unmodeled systematic error. In the sections below we consider two potential astrophysical explanations for the TTVs: a planetary satellite or an external perturber.
TTVs generated by a planetary satellite
HATP2b is expected to be in pseudosynchronous rotation, in which its spin period is close to the orbital frequency in the vicinity of periapse. Given , Hut (1981)’s treatment gives
(24) 
In order to maintain orbital dynamical stability, a prospective moon for HATP2b must have have an orbital semimajor axis, , such that , where is the Hill Sphere radius at periapse,
(25) 
and (see, e.g. Barnes & O’Brien (2002) for a discussion of satellite stability for planets with low orbital eccentricity). At distance from HATP2b, the orbital period of the satellite is d, which is substantially shorter than the 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 satellitestoplanet 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 . The characteristic time for a satellite’s orbital decay, assuming an equilibrium theory of tides with a frequencyindependent tidal quality factor, , is (Murray & Dermott, 1999)
(26) 
Adopting , and we find Gyr, which is comparable to the estimated stellar age, Gyr (Pál et al., 2010). The small value for the apsidal Love number, , stemming from the fact that HATP2b at 8 is quite centrally condensed in comparison to Jupiter Bodenheimer et al. (see 2001, for a discussion of the relation between and interior models); Batygin et al. (see 2009, for a discussion of the relation between and interior models). Given the substantial range of possible values for , it would therefore not be unreasonable to find that HATP2b harbors a fairly massive moon.
A moon with the above qualities would produce transit timing variations of magnitude (Kipping, 2009)
(27) 
where the geometric factor appropriately amplifies or damps the transit timing variations in accordance with the eccentricity and orientation of the elliptical orbit of the planet about the parent star
(28) 
Substituting the various values discussed above, we find
(29) 
which is several orders of magnitude too small to be of current interest. So while (somewhat surprisingly) it is distinctly possible that HATP2b could currently harbor a massive satellite, any such satellite cannot be responsible for transit timing variations of the size that appear to be present.
TTVs generated by an external perturber with a long period
Perturbations from an asyet undetected perturbing body present another potential explanation for the apparent transit timing variations. After the signature of HATP2b’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, , and semimajor axis, , would suffice. For a configuration in which (Holman & Murray, 2005),
(30) 
where . For a perturber with , and , we find . 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.
A lowmass resonant perturber?
There are an effectively infinite number of stable orbits for perturbing bodies in meanmotion resonance with HATP2b, and the diversity of such orbits is extended by HATP2b’s substantial eccentricity. A body in loworder resonance with HATP2b can readily induce transit timing variations of the magnitude that are apparently observed, and may be able to produce the curious structure exhibited by the timing variations of the secondary eclipses and the primary transits. Exploratory calculations are currently underway to evaluate this possibility. If we assume that the observed TTV are generated by a gravitational perturber, this appears to be the most promising approach. However, we note that the significance of the reported transit timing variations is less than 4. Further highprecision transit and/or secondary eclipse observations along with additional RV measurements would be needed to confirm the presence of these TTVs and constrain the properties of the perturbing body.
4.4. Transit and Secondary Eclipse Depths
From our three transit observations we obtain planetstar radius ratios of 0.068210.00075, 0.070410.00060, and 0.06680.0016 at 3.6, 4.5, and 8.0 m respectively. Our estimates of the planet/star radius ratio, , 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. (2007b), Winn et al. (2007a), and Loeillet et al. (2008). Although Pál et al. (2010) incorporate the band observations from Bakos et al. (2007b) they also include follow up observations that utilize the and Strömgren bandpasses. While the and bands probe a similar wavelength range (m), the Strömgren and bands probe a slightly shorter wavelength range (m) where atmospheric scattering may become important.
The difference between our values of 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 differ by less than 2. Given the large value of the average gravitational acceleration of HATP2b (16227 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 (Pál et al., 2010), then we find an average radius value for HATP2b of . The radius of 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 for 8 planets are expected for 110 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 HATP2b. We compare these results to the predictions of atmospheric models for this planet to better understand the atmospheric properties of HATP2b. Figure 14 presents predicted day side emission as a function of wavelength from a onedimensional radiative transfer model similar to those described in Burrows et al. (2008). These models assume a solarmetallicity 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 () and the possible presence of a highaltitude optical absorber () It is important to note that these onedimensional 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.
Models that also include an additional highaltitude 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 planetwide 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 , which corresponds to a brightness temperature of K We find an average 4.5 m eclipse depth of and a corresponding brightness temperature of K. Our eclipse depth at 5.8 m corresponds to a brightness temperature of K , while our eclipse depth at 8.0 m corresponds to a brightness temperature of 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 HATP2b at 3.6 and 4.5 m could place better constraints on possible orbittoorbit variability in the planet’s thermal structure.
4.5. Phase Curve Fits
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 HATP2b. For planets on circular orbits, phase curve observations are generally related to day/night brightness contrasts on the planet. In the case of HATP2b, the phase variations in the planetary flux are more indicative of thermal variations that result from the timevariable 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 HATP2b’s atmosphere.
We find that the peak of the observable planetary flux at 3.6 m occurs hours after periapse passage with a peak value of , which corresponds to a brightness temperature of K. The exact timing and level of the minimum planetary flux at 3.6 m is much more uncertain. We find the the planetary flux at 3.6 m drops below observable levels for a period of 1 day 18.36 hours before transit. As such we can only put an upper limit of 1040 K on the minimum brightness temperature of HATP2b at 3.6 m.
The observed 4.5 m HATP2b phase curve exhibits a peak in the observable planetary flux hours after periapse passage with a value of and corresponding brightness temperature of K. We detect emission from HATP2b over the entirety of its orbit at 4.5 m with a minimum in the planet/star flux ratio of , which corresponds to a brightness temperature of K. This minimum in the planetary flux occurs hours after transit. Roughly speaking, the shift of the minimum of the observed phase curve at 4.5 m away from 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 nightside of the planet is still cooling even after the transit event. We would expect the planet to have a temperature minimum east of the antistellar point near the sunrise terminator assuming a superrotating 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 HATP2b’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.640.33 hours after the predicted time for periapse passage. The maximum of the planetary flux at 8.0 m is , which corresponds to a brightness temperature of 280679 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 WASP12b 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, . A tidallylocked planet has 1to1 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 that phase mapping should not work for eccentric planets: 1) the timevariable incident flux makes the brightness map timevariable (), and 2) the timevariable 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 () rather than in time, which implicitly accounts for 2).
The fact that Equation 12 fits the phase variations of HATP2b implies that the diurnal brightness profile of the planet is constant throughout the orbit: . 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: HATP2b shows us its dayside shortly after periapse. The dayside brightness of the planet likely changes throughout its orbit, but we are not privy to it.
Interpreting the Flux Maximum
We use the semianalytic model developed in Cowan & Agol (2011) to interpret the peak amplitudes and phase lags of the
planet’s thermal brightness variations. This is essentially a twodimensional, onelayer energy balance model where the
user specifies Bond albedo, radiative timescale at the substellar point at periapse, and the characteristic
zonal
The Bond albedo is assumed to be 0.1 for all of the simulations shown here. Measured albedos of hot Jupiters range from 0.025 for TrES2b
(Kipping & Spiegel, 2011) to 0.32 for Kepler7b (Demory et al., 2011). Since HATP2b is viewed equatoron, its albedo is degenerate with poleward heat
transport, which we do not explicitly model. Increasing either albedo or meridional
The substellar radiative timescale is specified at periapse, and is assumed to scale as (Showman & Guillot, 2002; Showman et al., 2010), as one would expect for blackbody parcels of gas. The radiative relaxation time therefore varies throughout the orbit and is also a function of the location of a parcel of gas on the planet.
Observationally, the zonal wind speeds on a gas giant like HATP2b are degenerate with its rotation rate. In our model we therefore specify the angular velocity of gas parcels in an inertial frame, rather than in the usual rotating planet frame. By adopting the Hut (1981) prescription for pseudosynchronous rotation of binary stars, we convert the inertial angular frequency into an equatorial zonal wind speed in the rotating planet’s frame. If another prescription is more appropriate for the planet’s rotation rate (e.g. Ivanov & Papaloizou, 2007), then the equatorial wind velocities presented here are off by a uniform offset.
Both zonal wind speeds and albedo are assumed to be constant throughout the planet’s orbit. The constant zonal wind velocity is likely the most problematic assumption given that threedimensional simulations of Kataria et al. (2012) predict that zonal wind speeds at the midIR photosphere (pressures of 0.11 bar) change by tens of percent throughout the orbit of a hot Jupiter with an eccentricity of 0.5. It is also likely that the amount of equatortopole heat transport, which is degenerate with albedo in our model, will vary throughout HATP2b’s orbit. Our assumption of a constant albedo for the planet could also be limiting given that clouds could form near the apoapse of HATP2b’s orbit, then dissipate near periapse (Kane & Gelino, 2010). By focusing on the region near periapse we limit the possible influence of temporal changes in albedo and wind speeds on our results.
Model of HATP2b and Circular Analog
The top panels of Figure 15 shows how and affect the 4.5 m peak flux ratio and phase lag for HATP2b. 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 HATP2b with semimajor axis fixed at the actual planetÕs periapse separation. There are a number of conclusions we can draw from these models:

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 , 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 vs. westward zonal winds.

The peak flux ratio for the eccentric planet HATP2b 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 hours 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 , 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, timeconstant temperature, as in the circular case.

Comparing the topright and bottomright 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 phaselag expected in the absence of winds, the null hypothesis. For a circular, tidallylocked 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 zeropoint, 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).
Constraints on Model Parameters
In Figure 16 we show exclusion diagrams in the vs. plane for each waveband. Since we use a onelayer model, this can roughly be thought of as constraining the properties of the photospheres at each waveband, with the understanding that the midIR vertical contribution functions span approximately a factor of ten in pressure (e.g. Knutson et al., 2009b; 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. Grayscale shows the combined confidence intervals.
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 (topright 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 . 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 (). 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 exponentiallike trend of the preferred values of the zonal wind speeds with in Figure 16 that pairs short values for with higher values for and longer values for with lower values for . Threedimensional models that consistently couple radiative and advective processes will help to further constrain the relevant timescales in HATP2b’s atmosphere over the full range of its orbit.
5. Conclusions
In this paper we present the first secondary eclipse and phase curve observations for the eccentric hot Jupiter HATP2b in the Spitzer 3.6, 4.5, 5.8, and 8.0 m bandpasses. These data include two fullorbit phase curves at 3.6 and 4.5 m, a partialorbit phase curve at 8.0 m, three transit events, and six secondary eclipse events. The timing between transit and secondary eclipse during a single planetary orbit combined with radial velocity measurements allows us to better constrain the eccentricity () and argument of periapse () of HATP2b’s orbit. Longterm linear trends in the RV data indicate the presence of a third body in the system.
A comparison of our secondary eclipse depths with a onedimensional model for the dayside emission of the planet suggests the presence of a dayside inversion in HATP2b’s atmosphere. The timing and magnitude of the peak in the planetary flux at 3.6, 4.5, and 8.0 m are explained by a range of advective and radiative parameters in our twodimensional energy balance model, but suggest the presence of strong eastward equatorial winds ( km s) and short radiative timescales ( hours) at midIR photospheric levels near periapse.
Further work is needed to fully explain our observations of HATP2b. Threedimensional atmospheric models that couple radiative and advective processes and include a range of compositions would help to further explain the phase variations we observe for HATP2b, especially outside of the orbital region near periapse. Additionally, lowresolution spectroscopy taken during secondary eclipse would help to better constrain the atmospheric chemistry of HATP2b. Exoplanets on highly eccentric orbits like HATP2b present observers and modelers with a unique opportunity to disentangle the contributions of radiative, advective, and chemical processes at work in hot Jupiter atmospheres. By refining our understanding of exoplanets like HATP2b we will be able to use that knowledge to better understand the atmospheric processes at work in other exoplanet atmospheres.
Appendix A Noise Pixel Parameter
The IRAC Point Response Function (PRF) results from the convolution of the stellar Point Spread Function (PSF) with the Detector Response Function (DRF). The shape of the PRF is not constant and varies with the DRF at each stellar centroid position and with changes in the stellar PSF determined by the optics. The shape of the PRFs differs substantially in each channel of the IRAC instrument as shown in Figure 17. Given this change in the shape of the PRF with wavelength, we expect that different methods to determine the stellar centroid position and correct for systematic effects may be needed in each of the four IRAC bands. Both the 3.6 and 4.5 m channels of the IRAC instrument are shortward of the Spitzer 5.5 m diffraction limit (Gehrz et al., 2007) and exhibit undersampled PRFs whose shape changes as the stellar centroid moves from the center to the edge of a pixel causing intrapixel sensitivity variations. The PRFs in the 5.8 and 8.0 m channels of the IRAC instrument exhibit a more Gaussianlike shape that is better sampled by the detector resulting in only small intrapixel sensitivity variations.
Ideally, we would like to account for these changes in the PRF in our photometric measurements. We cannot make a direct determination of the exact shape of the PRF as a function of pixel position, but we can measure the normalized effectivebackground area of the PRF (), which is also called noise pixels by the IRAC Instrument Team. If we assume the measured flux in a given pixel () is given by
(A1) 
where is the point source flux and is PRF in the pixel. The noise pixel parameter, , is given by
(A2) 
since can be assumed to be constant. The numerator in Equation A2 is simply the square of the PRF volume () and the denominator is the effective background area (). These quantities are related to the sharpness parameter, , first introduced by Muller & Buffington (1974) for constraining adaptive optics corrections by
(A3) 
The parameter describes the ‘sharpness’ of the PRF and can range from zero (flat stellar image) to one (all the stellar flux in the central pixel).
As discussed in Mighell (2005), the parameter is related to the standard deviation () of the PRF by
(A4) 
where is a constant that depends on the sampling of the PRF on the detector. From Equations A4 and A3 it can be shown that
(A5) 
For both the 3.6 and 4.5 m data we measure with the same circular aperture sizes used to determine the stellar centroid position.
We find that can serve two purposes in improving the signaltonoise of the 3.6 m observations. First, using a circular aperture with a radius proportional to reduces the overall variations in raw unbinned flux from 5% to 2%. Harris (1990) and Pritchet & Kline (1981) note that the optimal aperture radius for a circularlysymmetric Gaussian with a standard deviation of is , which is similar to our optimal aperture radius . Second, we find that using as an additional spatial constraint in the intrapixel sensitivity correction at 3.6 m can improve the accuracy, as defined by the standard deviation of the residuals, in our final solution by 12%. We find that using as an additional constraint for the 4.5 m observations does not significantly improve the accuracy of our results. This is not surprising since the IRAC 4.5 m channel is closer to the Spitzer diffraction limit of 5.5 m (Gehrz et al., 2007). We also find that does not vary with stellar centroid position in the 5.8 and 8.0 m observations, which are longward of the Spitzer of the diffraction limit.
Appendix B Intrapixel Sensitivity Correction
The 3.6 and 4.5 m channels of the Spitzer IRAC instrument both exhibit variations in the measured flux that are strongly correlated with the position of the star on the detector (Figures 1 and 2). These flux variations are the result of well documented intrapixel sensitivity variations that are exacerbated by an undersampled PRF (e.g., Reach et al. (2005), Charbonneau et al. (2005, 2008), MoralesCalderón et al. (2006), Knutson et al. (2008)). The most common method used to correct intrapixel sensitivity induced flux variations is to fit a polynomial function of the stellar centroid position (Reach et al., 2005; Charbonneau et al., 2005, 2008; MoralesCalderón et al., 2006; Knutson et al., 2008). This method works reasonably well on short timescales ( 10 hours) where the variations in the stellar centroid position are small ( 0.2 pixels). Recently, studies by Ballard et al. (2010) and Stevenson et al. (2011) have implemented nonparametric corrections for intrapixel sensitivity variations by creating pixel ‘maps’, which give a smaller scatter in the re