3.6 AND 4.5 μm PHASE CURVES OF THE HIGHLY IRRADIATED ECCENTRIC HOT JUPITER WASP-14b

Each full-orbit observation contains one transit and two secondary eclipses. We model these events using the formalism of Mandel & Agol (2002). The transit light curve includes four free parameters: the scaled orbital semimajor axis a/R_*, the inclination i, the center of transit time t_T, and the planet–star radius ratio R_p/R_*, which is the square root of the relative transit depth. Each secondary eclipse event is defined by a center of eclipse time t_E and a relative eclipse depth d, measured with respect to the value of the phase curve at mid-transit, in addition to a/R_* and i, thus yielding four additional free parameters: t_E1, t_E2, d₁, and d₂. We ensure continuity between the phase curve and secondary eclipse light curves by scaling the amplitudes of eclipse ingress and egress (when the planet is partially occulted by the star) appropriately to match the out-of-eclipse phase curve values at the start and end of the eclipse. The host star WASP-14 has an effective temperature T_* = 6462 ± 75 K, a specific gravity and a metallicity of [Fe/H] = −0.13 ± 0.08 (Torres et al. 2012). We model the limb-darkening in each bandpass using a four-parameter nonlinear limb-darkening law with parameter values calculated as described in Sing (2010) for a 6500 K star with and [Fe/H] = −0.10: at 3.6 μm and at 4.5 μm.

WASP-14b has a relatively low eccentricity of 0.08, and therefore the variation in the planet's apparent brightness throughout an orbit can be modeled to first order as a simple sinusoidal function of the true anomaly f (Lewis et al. 2013), analogous in form to a simple sine or cosine of the orbital phase angle that is used in the case of a circular orbit (Cowan & Agol 2008):

where F₀ is the star's flux (assumed to be constant and normalized to one), c₁ is the amplitude of the phase variations, and c₂ represents the lag between the peak of the planet's temperature and the time of maximum incident stellar flux due to the finite atmospheric radiative timescale of the planet. Here, c₁ and c₂ are free parameters that are computed in our fits. We also experimented with other functional forms of the phase curve that included higher harmonics, but all of them resulted in higher values of the Bayesian Information Criterion (BIC) (see Section 3.4 for more information).

Photometric data obtained using Spitzer/IRAC in the 3.6 and 4.5 μm bandpasses exhibit a well-studied instrumental effect due to intrapixel sensitivity variations (Charbonneau et al. 2005). Small changes in the telescope pointing during observation cause variations in the measured flux from the target, resulting in a characteristic sawtooth pattern in the raw extracted photometric series. In our analysis, we decorrelate this instrumental systematic in two ways.

Our first approach to removing the intrapixel sensitivity effect is called pixel mapping (Ballard et al. 2010; Lewis et al. 2013). In an image j, the location of the target on the array is given by the measured centroid position (x_j, y_j), and the sensitivity of the pixel at that location is determined by comparing other images with measured centroid positions near (x_j, y_j). The effective pixel sensitivity at a given position is calculated as follows:

Here, F_meas,j is the flux measured in the jth image and F_j is the intrinsic flux; x_j, y_j, and are the measured x position, y position, and noise pixel parameter values. The quantities σ_x,j, σ_y,j, and are the standard deviations of x, y, and over the full range in i. For each image j, the summation in Equation (3) runs over the nearest m = 50 neighbors of the stellar target, where we define distance as

This method in effect adaptively smoothes the raw pixel map, allowing for a finer spatial scale in regions where the density of points is high while using a coarser spatial scale in regions with sparser sampling. We chose this number of neighbors to be large enough to adequately map the pixel response while maintaining a reasonably low computational overhead (Lewis et al. 2013). Several previous studies of Spitzer phase curves (e.g., Knutson et al. 2012; Zellem et al. 2014) do not include the noise pixel parameter term in Equation (3); we find that for our WASP-14b data, including the noise pixel parameter term produces ∼5%–10% smaller residual scatter from our best-fit light curve solution in both bandpasses.

The second approach is a recently proposed technique known as pixel-level decorrelation (PLD; see Deming et al. 2015, for a complete description). Unlike most other treatments of the intrapixel sensitivity effect, PLD does not attempt to relate the variations in the calculated position of the target on the pixel to the apparent intensity fluctuations. Rather, it utilizes the actual measured intensities of the individual pixels spanning the stellar PSF to provide an expression of the total measured flux. We consider pixels lying in a 3 × 3 box centered on the star, which have pixel intensities P_k(t), k = 1, ...9. We divide each 3 × 3 pixel box by the summed flux over all nine pixels in order to remove (at least to zeroth order) any astrophysical flux variations, giving the following relation:

The intrapixel sensitivity effect is modeled as a linear combination of the arrays and thus the total measured intensity S is given by

where F(t) is the astrophysical model, comprising the phase curve, transit, and eclipses. The parameters b_k are the linear coefficients that are determined through least-squares fitting, and h is a free normalization parameter that corrects for the overall numerical offset introduced by the linear sum of pixel intensity arrays. Following Deming et al. (2015), we also include a linear ramp in time with a slope parameter b_t.

In Deming et al. (2015), PLD was applied to fitting Spitzer secondary eclipses of several exoplanets, including WASP-14b. PLD was found to be generally more effective in removing time-correlated (i.e., red) noise, resulting in lower residual scatter from the best-fit solution when compared with other decorrelation techniques. The best results were obtained when the photometric time series were binned, pixel by pixel, prior to PLD fitting, since binning improves the precision of the measured intensities in pixels at the edge of the stellar PSF and can be adjusted to reduce the noise on the timescale of interest. When fitting our full-orbit WASP-14b photometric series using the PLD technique, we experimented with fitting either unbinned or binned data, with bin sizes equal to powers of two up to 256 (∼8.5 minutes). We found that larger bin sizes, comparable to or exceeding the occultation ingress/egress timescale of roughly 20 minutes, cause excessive loss of temporal resolution and yielded unsatisfactory secondary eclipse and transit light curve fits. After fitting the model light curve to the data using PLD, we subtract the best-fit solution from the raw unbinned data to produce the residual series, which we use to evaluate the relative amount of time-correlated noise remaining in the data.

In Figure 3, we compare the noise properties of each version of the 3.6 μm residual time series to the ideal scaling we would expect for the case of independent (i.e., "white" noise) Gaussian measurement errors, where n is the bin size. The white noise trend is normalized at bin size n = 1 to the photon noise level corresponding to the median photon count over the observation data set (with sky background included). Comparing the various PLD fits, we find that at small bin sizes, unbinned PLD gives the lowest residual scatter, while at larger bin sizes approaching the duration of eclipse ingress or egress, PLD with larger bins results in lower residual scatter. The same trend is seen when fitting the 4.5 μm data. Deming et al. (2015) found that the optimal PLD performance is achieved when the range of star positions is lower than 0.2 pixels. The range of pixel motion in our full-orbit phase curve observations modestly exceeds this limit, and we indeed find that the residual scatter from the pixel mapping fits is smaller than the scatter from any of the PLD fits. We conclude that the pixel mapping technique produces the lowest residual scatter for our full-orbit observation data sets, and we therefore use this technique in the final version of our analysis.

Zoom In Zoom Out Reset image size

In addition to having higher residual scatter than the pixel mapping solutions, the PLD fits yield eclipse depth and phase curve parameter estimates that often differ strongly from the corresponding values derived from fits with pixel mapping; these discrepancies sometimes exceed the 3σ level. The best-fit values derived from the PLD fits also display a higher level of variation across different choices of binning, photometric aperture, and exponential ramp type (see Section 3.4) than in the case of pixel mapping fits. This points toward an inherent instability in the PLD method when fitting full-orbit phase curves that may be related to the larger range in star motions characteristic of such data sets.

Previous studies using Spitzer/IRAC have noted a short-duration ramp at the beginning of each observation, and again after downlinks (e.g., Knutson et al. 2012; Lewis et al. 2013). The ramp has a characteristic asymptotic shape that typically decays to a constant value on timescales of an hour or less in the 3.6 and 4.5 μm bandpasses. We first experimented with removing the first 30 or 60 minutes of data from each phase curve observation, selecting the removal interval that minimizes the residual rms from the best-fit solution binned in five-minute intervals. We find that in both bandpasses, we obtain the best results when we do not trim any data from the start of the observations. When examining the residual time series, we noticed a small ramp visible at the start of the 3.6 μm observations. We therefore considered whether or not our fits might be further improved by the addition of an exponential function.

We experimented with including an exponential ramp in our phase curve model, using the formulation given in Agol et al. (2010):

where t is the time since the beginning of the observation, and a₁–a₄ are correction coefficients. To determine whether this function is necessary and if so, how many exponential terms to include in the ramp model, we use the BIC, defined as

where k is the number of free parameters in the fit, and N is the number of data points. By minimizing the BIC, we select the type of ramp model that yields the smallest residuals without "over-fitting" the data. For the 3.6 μm data, we find that using a single exponential ramp gives a marginally lower BIC compared to the no-ramp case, while for the 4.5 μm data, no ramp is needed at all. The residuals from the best-fit full phase curve solution, shown in Figures 4 and 5, do not appear to display any uncorrected ramp-like behaviors.

Zoom In Zoom Out Reset image size

We use a Levenberg–Marquardt least-squares algorithm to fit each full-orbit photometric series to our total model light curve, with the intrapixel sensitivity correction calculated via pixel mapping. In the final version of these global fits, we use the updated values for e and ω obtained from our radial velocity analysis (see Section 4.1) as well as the updated orbital period calculated from fitting all published transits. The best-fit transit, secondary eclipse, and phase parameters are listed in Table 1 along with their uncertainties. Figures 4 and 5 show the full-orbit data in the 3.6 and 4.5 μm bandpasses, respectively, with instrumental variations removed. The individual eclipse and transit light curves are shown in Figures 6 and 7.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1.  Best-fit Parameters

Parameter	3.6 μm	4.5 μm
Transit Parameters
Rp/R*
tT (BJD)a
Eclipse Parameters
1st eclipse depth, d1 (%)
tE1 (BJD)a,b
2nd eclipse depth, d2 (%)c
tE2 (BJD)a,b,c
Orbital Parameters
Inclination, i (°)
Scaled semimajor axis, a/R*		5.98 ± 0.13
Phase Curve Parameters
Amplitude, c1 (×10−4)
Phase shift, c2 (°)
Maximum flux offset (h)d	−1.43 ± 0.21	−1.01 ± 0.21
Minimum flux offset (h)d
Ramp Parameters
a1 (×10−4)		...
a2 (d)		...

Notes.

^aAll times are listed in BJD_UTC for consistency with other studies; to convert to BJD_TT, add 66.184 s (Eastman et al. 2010). ^bThe center of secondary eclipse times are not corrected for the light travel time across the system, Δt = 35.9 s. ^cThese values are computed from fitting the second 3.6 μm secondary eclipse separately using the pixel-level decorrelation (PLD) method. This was done in order to remove an anomalous signal that occurs during the eclipse in the global 3.6 μm phase curve fit using pixel mapping. ^dThe maximum and minimum flux offsets are measured relative to the center of secondary eclipse time and center of transit time, respectively, and are derived from the phase curve fit parameters c₁ and c₂. The maximum flux offset reported is the error-weighted mean of the flux offsets relative to the first and second secondary eclipses.

Download table as:  ASCIITypeset image

In the 3.6 μm light curve data, there is an anomalous signal in the residuals from the global phase curve fit that occurs during the last secondary eclipse. This anomaly is characterized by a short ∼20 minute dip in the middle of the eclipse and was not removed by the pixel mapping method, resulting in a 3.1σ discrepancy between the eclipse depths from the first and second events. Similar short-duration anomalies (both positive and negative) have been reported in Spitzer 3.6 μm data and are usually attributed to variations in the width of the stellar target (e.g., Lanotte et al. 2014). Our pixel mapping technique accounts for these variations by incorporating the noise-pixel parameter pixel sensitivity calculation in Equation (3), so we conclude that the anomaly in our data is likely attributable to some other instrumental effect.

In order to recover the second 3.6 μm eclipse depth, we experimented with fitting the eclipse event separately. Selecting a short ∼0.2-day segment of the phase curve observation surrounding the eclipse, we fit the data to a simplified secondary eclipse light curve model using both pixel mapping and PLD. This model has three free parameters—the center of eclipse time t_E2, the eclipse depth d₂, and a linear slope c_t, which accounts for the out-of-eclipse variation in the planet's brightness. The inclination, scaled semimajor axis, and planet–star radius ratio are fixed at the values derived from the global 3.6 μm phase curve fit and listed in Table 1. In these fits, we use the same choice of aperture as in the full phase curve fit; in the case of PLD, we optimize for the bin size based on the residual scatter and find that a bin size of 128 yields the lowest residual rms. Comparing the results of our fits using pixel mapping and PLD, we see that while the residuals from the best-fit solution with pixel mapping show a significant anomalous signal similar to the one present in the global 3.6 μm phase curve fit, no residual anomaly is evident in the best-fit solution with PLD. The range of star positions during this segment of data is much less than 0.2 pixels, so we expect PLD to perform optimally. The best-fit parameter values from the PLD fit are (BJD_UTC), and The data with instrumental effects removed and the best-fit secondary eclipse solution are shown in Figure 8. We note that the best-fit eclipse depth from our individual fit of the second 3.6 μm eclipse is consistent with the depth of the first 3.6 μm eclipse derived from the global phase curve fit at the 0.15σ level. Henceforth, we use the values from our individual secondary eclipse fit with PLD for the center of eclipse and eclipse depth of the second 3.6 μm eclipse.

Zoom In Zoom Out Reset image size

Taking the error-weighted average of the eclipse depths listed in Table 1 at each wavelength, we arrive at 0.1882% ± 0.0048% and 0.2247% ± 0.0086% for the 3.6 and 4.5 μm bandpasses, respectively. These values are consistent with the ones reported in Blecic et al. (2013) in their analysis of previous Spitzer secondary eclipse observations to better than 1σ (0.19% ± 0.01% at 3.6 μm and 0.224% ± 0.018% at 4.5 μm).

We estimate the uncertainties in our best-fit parameters in two ways. The first approach is the "prayer-bead" (PB) method (Gillon et al. 2009), which gives an estimate of the contribution of time-correlated noise to the uncertainty. This method entails extracting the residuals from the best-fit solution, dividing the residuals into segments, and cyclically permuting the residual series segment by segment, each time adding the new residual series back to the best-fit solution and recomputing the parameters using the least-squares algorithm. For each free parameter, we create a histogram of the best-fit values from every permutation and calculate the uncertainties based on the 1σ upper and lower bounds from the median.

The second approach is a Markov chain Monte Carlo (MCMC) routine with 10⁵ steps, where we initiate each chain at the best-fit solution from the least-squares analysis. The uncertainty on individual data points is set to be the standard deviation of the residuals from the best-fit solution. We discard an initial burn-in on each chain of length equal to 20% of the chain length, which we found ensured the removal of any initial transient behavior in a chain, regardless of the choice of initial state. The distribution of values for each parameter is close to Gaussian, and there are no significant correlations between pairs of parameters. As in the PB method, we set the uncertainties in the fitted parameters to be the 1σ upper and lower bounds from the median. For each parameter, we choose the larger of the two errors and report it in Table 1. We find that the PB errors are consistently larger and range between 1.0 and 2.7 times that of the corresponding MCMC errors.

The rms scatter in the best-fit residual series exceeds the predicted photon noise limit by a factor of 1.16 at 3.6 μm and 1.14 at 4.5 μm. We estimate the level of red noise by calculating the standard deviation of the best-fit residuals for various bin sizes, shown in Figure 9 along with the inverse square-root dependence of white noise on bin size for comparison. On timescales relevant for the eclipses and transits (e.g., the ingress/egress timescale—∼15 minutes), the red noise increases the rms by a factor of approximately 1.8 at 3.6 μm and approximately 1.4 at 4.5 μm.

Zoom In Zoom Out Reset image size

We combine the two transit times calculated from our phase curve observations with all other published values (Joshi et al. 2009; Blecic et al. 2013; Raetz et al. 2015) to arrive at an updated ephemeris for the WASP-14 system. Here, we define the zeroth epoch as that of the transit nearest in time to the error-weighted mean of all measured transit times. The transit observations span more than 5.3 years, and by fitting a line through the transit times, we derive new, more precise estimates of the orbital period P and mid-transit time T_c,0:

Figure 10 shows the observed minus calculated transit times derived from these updated ephemeris values.

Zoom In Zoom Out Reset image size

We use the secondary eclipse times to obtain a second, independent estimate of the orbital period. Carrying out a linear fit through the four secondary eclipse times calculated from our phase curve data and the two secondary eclipse times published in Blecic et al. (2013), we arrive at a best-fit value of P = 2.2437660 ± 0.0000017 days. This period is consistent with the best-fit transit period at the 0.5σ level. In Figure 11, we plot the orbital phase of secondary eclipse for all published secondary eclipse times using the updated transit ephemeris values in Equation (9). The error-weighted mean orbital phase of secondary eclipse is 0.48412 ± 0.00013.

Zoom In Zoom Out Reset image size

By combining the updated transit ephemeris and secondary eclipse times derived from the global phase-curve fits with the radial velocity measurements analyzed in Knutson et al. (2014), we can obtain new estimates of the orbital eccentricity and pericenter longitude: and . These values are consistent with those reported in Knutson et al. (2014). We use the updated values of orbital eccentricity and pericenter longitude in our final phase curve fits (Table 1). Results from our radial velocity fits are shown in Table 2 and Figure 12. These fits provide new estimates of the orbital period (P_b) and center of transit time (T_c,b), as well as the semi-amplitude of the planet's radial velocity (K_b), the radial velocity zero-point offsets for data collected by each of the different spectrographs from which radial velocity measurements of the system were obtained (γ_1–4), and the slope () of the best-fit radial velocity acceleration. The radial velocity slope is consistent with zero, indicating no evidence for additional planets in the WASP-14 system.

Zoom In Zoom Out Reset image size

Table 2.  Results from Radial Velocity Fit with Priors on Transit Ephemeris and Eclipse Times

Parameter	Value	Units
RV Model Parameters
Pb	2.24376524 ± 4.4E − 07	days
Tc,b	2456034.21261 ± 0.00015
eb	0.0830	⋯
ωb	252.67	degrees
Kb	986.4	m s−1
γ1	155.8 ± 2.6	m s−1
γ2	−70.0	m s−1
γ3	−73.7 ± 6.8	m s−1
γ4	187.4	m s−1
	0.0025 ± 0.0032	m s−1day−1
RV-derived Parameters
	−0.02474	⋯
	−0.0792	⋯

Note. Radial velocity zero-point offsets (γ_1–4) derived from four separate RV data sets: 1—Keck/HIRES (Knutson et al. 2014), 2—FIES (Joshi et al. 2009), 3—SOPHIE (Joshi et al. 2009), 4—SOPHIE (Husnoo et al. 2011). See the text for a description of other variables.

Download table as:  ASCIITypeset image

From the RV fits we also arrive at updated values of the orbital semimajor axis and planet mass: a = 0.0371 ± 0.0011 AU and M_p = 7.76 ± 0.47 M_Jup. Using the error-weighted best-fit values of a/R_* and R_p/R_* from both bandpasses, we obtain a new estimate of the planet's radius: R_p = 1.221 ± 0.041 R_Jup. A full list of updated planetary parameters is given in Table 3.

Table 3.  Updated Planetary Parameters

Parameter	Value
Rp/R*	0.09419 ± 0.00043
a/R*	5.99 ± 0.09
i (°)	84.63 ± 0.24
e
ω (°)
Mp (MJup)	7.76 ± 0.47
Rp (RJup)	1.221 ± 0.041
ρp (g cm−3)	5.29 ± 0.62
gp (m s−2)	129 ± 12
a (AU)	0.0371 ± 0.0011

Download table as:  ASCIITypeset image

In this section, we combine the results of our global fits with model-generated spectra and light curves to provide constraints on the atmospheric properties of the planet. To compute the relative planet–star flux ratio phase curve F_p/F_*, we subtract the secondary eclipse depth from the best-fit phase curve and divide by the remaining flux measured at the center of secondary eclipse, which represents the brightness of the star alone. For the secondary eclipse depths in each bandpass, we take the error-weighted average of the depths listed in Table 1. The resulting phase curves are shown in Figure 13 along with the corresponding ±1σ brightness bounds.

Zoom In Zoom Out Reset image size

Some of the main quantitative characteristics of the phase curves are summarized in Table 4. Notably, the fluxes at both Spitzer wavebands peak significantly before secondary eclipse, implying that the hottest regions are shifted eastward from the substellar point. This behavior appears to be common on hot Jupiters and, in addition to WASP-14b, has also been clearly observed on HD 189733b (Knutson et al. 2007, 2009a, 2012), HD 209458b (Zellem et al. 2014), HAT-P-2b (Lewis et al. 2013), WASP-43b (Stevenson et al. 2014), WASP-12b (Cowan et al. 2012), and Ups And b (Crossfield et al. 2010).

Table 4.  Phase Curve Comparison

Source	3.6 μm	4.5 μm
Maximum flux ratio [%]
Measured
Eq. chem. model	0.2124	0.2600
No TiO model	0.1461	0.1667
Minimum flux ratio [%]
Measured	<0.0175a
Eq. chem. model	0.0178	0.0318
No TiO model	0.0403	0.0512
Phase curve amplitude [%]b
Measured	>0.1702a
Eq. chem. model	0.1946	0.2282
No TiO model	0.1058	0.1155
Maximum flux offset [h]c
Measured	−1.43 ± 0.21	−1.01 ± 0.21
Eq. chem. model	+0.3	+0.2
No TiO model	−0.9	−1.0
Minimum flux offset [h]c
Measured
Eq. chem. model	−3.7	−2.9
No TiO model	−6.2	−6.2

Notes.

^aValues based off of the 2σ upper limit of the flux ratio minimum. ^bDifference between maximum and minimum flux ratios. ^cThe maximum and minimum flux offsets are measured relative to the center of secondary eclipse time and center of transit time, respectively, and are derived from the phase curve fit parameters. Negative time offsets of the maximum or minimum flux indicate an eastward shift in the location of the hot or cold region in the planet's atmosphere.

Download table as:  ASCIITypeset image

This behavior was predicted before the Spitzer era (Showman & Guillot 2002; Cooper & Showman 2005) and has now been reproduced in a wide variety of general circulation models (GCMs, e.g., Showman et al. 2009, 2015; Lewis et al. 2010; Rauscher & Menou 2010, 2012, 2013; Heng et al. 2011a, 2011b; Perna et al. 2012; Dobbs-Dixon & Agol 2013). In these models, the eastward hotspot displacement results from eastward advection due to a fast, eastward-flowing jet stream centered at the equator. Given such a jet, a significant hot spot offset occurs under conditions when the radiative timescale at the photosphere is comparable to timescales for air to advect horizontally over a planetary radius. The eastward offsets in our observations thus provide evidence that WASP-14b exhibits an eastward equatorial jet stream.

It is interesting to compare our observed day–night flux differences with those of other planets observed to date. As mentioned previously, Cowan & Agol (2011) and Perez-Becker & Showman (2013) inferred that, in general, planets that receive higher stellar fluxes exhibit larger fractional day–night temperature differences and less efficient day–night heat redistribution than planets that receive lower stellar fluxes. Averaging the phase curve amplitudes at 3.6 and 4.5 μm (listed in Table 4), we find that WASP-14b is roughly consistent with this trend, with a day–night flux difference that is smaller than those of highly irradiated planets like WASP-12b and WASP-18b, but larger than those of more weakly irradiated planets like HD 189733b and HD 209458b.

We consider four types of atmosphere models. First, we use an interpolated PHOENIX spectrum for the host star WASP-14 (Husser et al. 2013) to calculate the brightness temperatures of WASP-14b in each band from the retrieved secondary eclipse depths. The best-fit brightness temperatures are at 3.6 μm and at 4.5 μm. We find that the 3.6 and 4.5 μm eclipse depths are consistent with a single blackbody temperature and derive an effective temperature of T_eff = 2402 ± 35 K from a simultaneous fit to both bandpasses. The predicted equilibrium temperature for WASP-14b assuming zero albedo is 2220 K if incident energy is re-radiated from the dayside only and 1870 K if the planet re-radiates the absorbed flux uniformly over its entire surface. We can also compare the computed brightness temperatures to the effective temperature of the dayside in the no-albedo, no-circulation limit assuming each region is a blackbody locally in equilibrium with the incident stellar flux: 2390 K (Cowan & Agol 2011). The high brightness temperatures from the blackbody fits therefore suggest that WASP-14b has a very hot dayside atmosphere and inefficient day–night recirculation.

Next, we compare our phase curves and emission spectra to theoretical models generated from the three-dimensional Substellar and Planetary Atmospheric Radiation and Circulation (SPARC) model to investigate the global circulation of WASP-14b. The SPARC model was specifically developed with the study of extrasolar planetary atmospheres in mind (Showman et al. 2009). The SPARC model couples the MITgcm (Adcroft et al. 2004) with the non-gray radiative transfer model of Marley & McKay (1999) in order to self-consistently calculate the amount of heating/cooling at each grid point. In this way, the SPARC model does not require advective tuning parameters often employed in one-dimensional radiative transfer models or prescribed pressure-temperature profiles utilized in Newtonian cooling schemes employed in other circulation models. Our use of an atmospheric model that fully couples radiative and dynamical processes is especially important for planets on eccentric orbits, such as WASP-14b, which experience time variable heating (Lewis et al. 2010, 2014; Kataria et al. 2013).

The SPARC models of WASP-14b presented here utilize the cubed-sphere grid (Adcroft et al. 2004) with a horizontal grid resolution of roughly 5625 in latitude and longitude (a so-called low-resolution C16 grid). In the vertical direction, the models span pressures ranging from 200 bar to 0.2 mbar with 39 layers evenly spaced in and a top layer that extends to zero pressure. Here we consider atmospheric compositions in thermochemical equilibrium both with and without the incorporation of the strong optical absorbers TiO and VO (hereafter, "equilibrium chemistry" and "no TiO," respectively). The presence of TiO/VO in our models allows for the development of a dayside inversion layer (Fortney et al. 2008a). All models adopt solar elemental ratios of heavy elements.

We assume that WASP-14b is in a pseudo-synchronous rotation state (P_rot ∼ 2.14 days; Hut 1981). As in the case of GJ 436b, WASP-14b's relatively low eccentricity means that assuming pseudo-synchronous rotation instead of synchronous rotation is not likely to strongly affect the global circulation patterns that develop (Lewis et al. 2010).

Figure 13 compares the best-fit phase curve derived from the Spitzer data in each bandpass with the band-averaged light curves generated from the SPARC model using the methods of Fortney et al. (2006). Model dayside and nightside spectra at the center of secondary eclipse and center of transit times are shown in Figure 14 along with the corresponding measured flux ratios. The corresponding dayside and nightside temperature–pressure profiles generated by the models are shown in Figure 15. For the dayside planetary emission, we combine the 3.6 and 4.5 μm flux ratios reported in the present work and previously in Blecic et al. (2013) to arrive at the error-weighted average values: 0.1886% ± 0.0043% at 3.6 μm and 0.2245% ± 0.0078% at 4.5 μm. We also include the measured 8.0 μm flux ratio (0.181% ± 0.021%) from Blecic et al. (2013). Table 4 compares the model-predicted maximum and minimum flux ratios and time offsets with the values derived from the data.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In both the 3.6 and 4.5 μm bandpasses, the measured dayside planetary emission lies between the equilibrium chemistry and no TiO models. A possible explanation is that the dayside atmosphere of WASP-14b may contain a sub-solar abundance of TiO/VO (possibly due to cold trapping), which would result in a weaker temperature inversion than the one predicted by the solar abundance equilibrium chemistry model (e.g., Fortney et al. 2008a; Spiegel et al. 2009; Madhusudhan et al. 2011; Parmentier et al. 2013). Meanwhile, the measured 8.0 μm dayside brightness is not reproduced by either of the models.

We report offsets in the time of maximum and minimum flux relative to the center of transit and center of secondary eclipse times, respectively, as derived from our phase curve analysis and the models (Table 4). Negative time offsets of the maximum or minimum flux indicate an eastward shift in the location of the hot or cold region in the planet's atmosphere. The no TiO model greatly overestimates the magnitude of the minimum flux offsets while giving good agreement with the maximum flux offsets at both wavelengths. The equilibrium chemistry model yields a less severe overestimation of the minimum flux offset, while predicting a maximum flux offset that is opposite the one measured from the data. We find that overall the equilibrium chemistry model yields the closer match. The no TiO model predicts a significant asymmetry in the phase curve around the time of minimum flux, which is not seen in the Spitzer data.

The most notable discrepancy between the best-fit phase curves and the SPARC model results is in the nightside planetary flux. Both the equilibrium chemistry and the no TiO model overestimate the planet's nightside brightness at 3.6 μm, while underestimating the brightness at 4.5 μm. The planet's very low 3.6 μm nightside emission indicates a higher atmospheric opacity at that wavelength than is predicted by the models. Meanwhile, the higher-than-predicted 4.5 μm nightside emission points toward a slight reduction in the atmospheric opacity at that wavelength relative to the models. One possible explanation for both of these trends is an increased C/O ratio. Moses et al. (2013) demonstrate that increasing the C/O ratio above equilibrium values leads to an excess of CH₄ and a depletion of CO. This enhances the opacity in the 3.6 μm bandpass where CH₄ has strong vibrational bands, while reducing the opacity in the 4.5 μm bandpass where CO has strong vibrational bands.

An enhanced C/O ratio can also be invoked to explain the higher amplitude of the 3.6 μm phase curve as compared with the 4.5 μm phase curve. In this scenario the 3.6 μm bandpass would probe lower pressure regions that have strong day/night contrasts due to a combination of short radiative timescales and the possible presence of a dayside temperature inversion. A relative depletion of CO would shift the 4.5 μm photosphere to higher pressures, where longer radiative timescales facilitate more efficient day–night heat transport, resulting in a reduced phase curve amplitude. We note, however, that recent spectroscopic surveys of M dwarfs suggest that the occurrence of high C/O ratios in stellar atmospheres appears to be low (see Gaidos 2015 and references therein). At the same time, models of disk chemistry predict that the C/O ratio of the solids and gas in the disk may vary as a function of disk radius (e.g., Öberg et al. 2011; Madhusudhan et al. 2014). In particular, different snowlines of oxygen- and carbon-rich ices are expected to result in systematic variations in the C/O ratio of gaseous material across the protoplanetary disk. Gas giants that derive most of their atmosphere from the gas disk outside of the water snowline but inside of the carbon dioxide snowline may therefore end up with high atmospheric C/O ratios. We also note that Madhusudhan et al. (2011) showed that planets with high C/O ratios should naturally have less TiO and correspondingly weaker temperature inversion, consistent with our previous discussion of WASP-14b's dayside emission. In addition, Blecic et al. (2013) show that atmospheric models with high C/O ratios often predict reduced emission at longer wavelengths (6–8 μm) and may explain the low 8.0 μm planetary flux ratio (see Figure 14).

An alternative explanation for the planet's low 3.6 μm nightside emission is the presence of high-altitude silicate clouds on the cold nightside in or above the photosphere. The formation of an equilibrium silicate cloud is possible at the photospheric pressures (∼100 mbar) probed on the nightside of WASP-14b by these observations (e.g., Visscher et al. 2010). However, it is expected that the presence of a thick cloud would suppress both the 3.6 and 4.5 μm nightside fluxes from the planet, and the resulting planetary emission spectrum would resemble a blackbody with an effective temperature corresponding to that of the cloud tops. While the low 3.6 μm nightside brightness temperature of 1079 K (2σ upper limit) is consistent with a cold emitting layer, the two nightside band fluxes together are not consistent with a single blackbody spectrum.

In this scenario, the high 4.5 μm nightside emission (corresponding to a brightness temperature of ) could be explained by introducing an emitting layer of CO in a warmer thermosphere that is situated above the cold cloud tops (Koskinen et al. 2013). However, it is unclear how a high-altitude temperature inversion might arise on the nightside, as there is no incident stellar irradiation and the efficiency of recirculation from the dayside should be relatively weak at low pressures. Therefore, an enhanced C/O ratio in WASP-14b's atmosphere provides a more straightforward explanation for the observed differences between WASP-14b's 3.6 and 4.5 μm phase curves. Further observations of WASP-14b, possibly with the Hubble Space Telescope's Wide Field Camera 3, would provide important additional information to constrain WASP-14b's composition and either support or refute the high atmospheric C/O ratio scenario posited here for WASP-14b. We note that Lewis et al. (2014) also suggested an enhanced C/O ratio to explain the differences seen between the 3.6 and 4.5 μm phase-curve observations of HAT-P-2b. Additional phase-curve observations of hot Jupiters may reveal similar trends and point to a fundamental piece of physics currently missing from exoplanet atmospheric models and/or planet formation theories. The SPARC model utilized in this work does not readily accommodate non-solar C/O ratios. The exploration of the effects of different C/O ratios on the atmospheric dynamics of hot Jupiters will be the topic of future work.

We also compare the Spitzer dayside and nightside emission with 1D models generated using the methods of Burrows et al. (2008) with varying recirculation. These models assume local thermodynamic equilibrium, solar composition, and a plane-parallel atmosphere. A heat sink is included at depth (between 0.003 and 0.6 bars) to redistribute heat from the dayside to the nightside. These models also incorporate a generalized gray absorber at low pressures (0 to 0.03 bars), which is parameterized with an absorption coefficient, κ_e, with units of cm² g⁻¹. This absorber enhances the opacity of the planet's atmosphere at optical wavelengths, raising the local atmospheric temperature and producing a high-altitude temperature inversion. A second dimensionless parameter P_n is used to specify the efficiency of energy redistribution, with P_n = 0.5 indicating complete redistribution and P_n = 0 signifying redistribution on the dayside only.

In Figure 16, we show the dayside and nightside spectra for various values of κ_e and P_n and compare them with the measured relative planetary brightnesses in the two Spitzer bands. Figure 17 shows the corresponding dayside and nightside temperature–pressure profiles generated by the model. In these models the presence of the gray absorber does not affect day–night recirculation, so the nightside spectra depend only on the value of P_n. We define the measured dayside brightness in each bandpass to be the measured secondary eclipse depth, while for the the nightside data points, we used the value of the best-fit relative planetary phase curve F_p/F_* calculated at the mid-transit time in each bandpass. Looking at the dayside spectra, we see that the κ_e = 0.2 and P_n = 0.1 model spectrum comes closest to matching the measured data points at 3.6 and 4.5 μm. This suggests that WASP-14b has poor day–night recirculation and a moderate thermal inversion in the dayside atmosphere.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

On the nightside, as in the SPARC models, none of the models reproduce the low measured 3.6 μm brightness, while the measured 4.5 μm flux ratio is most consistent with the P_n = 0.3 model. As discussed earlier in the context of the SPARC models, this discrepancy is consistent with the hypothesis that WASP-14b might have an enhanced C/O ratio. We also note that large deviations from solar-like equilibrium chemical composition could yield a larger disparity in the pressures probed on the dayside and nightside than is accounted for in our models. Specifically, one could be probing a deeper than expected pressure level on the nightside, where circulation is efficient (high P_n), while probing a higher-up (lower pressure) level on the dayside, where the circulation is less efficient (low P_n). This scenario is especially probable when there is a dayside temperature inversion, as our data suggest.

Finally, we compare the albedo and recirculation derived from our best-fit phase curves to other hot Jupiters for which full-orbit thermal measurements have been obtained. Following the methods described in Schwartz & Cowan (2015), we use the measured eclipse depths and phase curve amplitudes to calculate the error-weighted dayside and nightside brightness temperatures, with corrections for the contamination due to reflected light: T_d = 2312 ± 35 K and T_n = 1299 ± 77 K. Although these methods were developed for planets on circular orbits, the low eccentricity of WASP-14b places it in a regime for which the model is still a reasonable approximation. From these temperatures, we derive a Bond albedo of A_B = 0 (<0.08 at 1σ) and a day–night heat transport efficiency of =0.23 ± 0.06, where the latter is defined such that = 0 means no heat recirculation to the nightside, and = 1 indicates complete redistribution. Figure 18 shows the location of WASP-14b in albedo-recirculation space along with six other exoplanets with measured thermal phase curves. WASP-14b's low day–night heat transport efficiency and high irradiation temperature are consistent with the general observed trend that highly irradiated hot Jupiters have poor heat recirculation (Cowan & Agol 2011; Perez-Becker & Showman 2013).

Zoom In Zoom Out Reset image size

Both WASP-14b and WASP-18b have much lower estimated albedos than those of other Jovian mass planets with thermal phase curve measurements; they are also both significantly more massive (7–10 versus ∼1 M_Jup). The low estimated albedos of these massive planets, implied by their high thermal emission, may instead be indicative of the detection of some amount of thermal radiation from the planet's interior in addition to the re-radiated stellar flux. A different cooling and/or migration history for these massive planets could result in an added contribution to the planetary emission from the deep atmosphere; the age of the WASP-14 system is relatively young (<1 Gyr; Joshi et al. 2009), so this additional internal flux may be due to significant residual heat of formation. In addition, the higher internal fluxes of these higher-mass planets may support stronger magnetic fields than on smaller planets, which may lead to stronger Ohmic dissipation in their atmospheres (e.g., Christensen et al. 2009; Batygin et al. 2013), thereby providing a source of additional emission on the dayside.