A NEW 24 μm PHASE CURVE FOR υ ANDROMEDAE b

Ian J. M. Crossfield; Brad M. S. Hansen; Joseph Harrington; James Y.-K. Cho; Drake Deming; Kristen Menou; Sara Seager

doi:10.1088/0004-637X/723/2/1436

1. INTRODUCTION

The first thermal characterizations of highly irradiated extrasolar planetary atmospheres were made by measuring the flux decrement that occurs during secondary eclipse, when an extrasolar planet passes behind its host star (Deming et al. 2005; Charbonneau et al. 2005). This decrement gives an estimate of the (hemisphere-averaged) temperature of a planet's star-facing side at the time of eclipse and provides insight into the energy budgets of these hot worlds. Secondary eclipse (or occultation) observations have been widely interpreted as indicating two types of planetary atmospheres, namely planets with and without high-altitude temperature inversions (Burrows et al. 2008; Fortney et al. 2008). However, the number of free parameters in current models is such that in many cases it is still difficult to place strong constraints on a planet's atmospheric structure with the few data points available for most systems (Madhusudhan & Seager 2009).

Burrows et al. (2008) and Fortney et al. (2008) suggested that sufficiently high levels of irradiation prevent strong optically absorbing species from condensing and "raining out" of the upper atmosphere of hot Jupiters, thus directly linking high levels of incident stellar flux to the presence of a temperature inversion. However, subsequent secondary eclipse measurements with Spitzer/IRAC have complicated the picture and a straightforward connection between irradiation and inversions now seems untenable. For example, TrES-3b receives substantially more flux than does HD 209458b, yet the latter has an atmospheric inversion (Knutson et al. 2008) while the former does not (Fressin et al. 2010). Thus, planetary classification will require more subtlety than a simple critical-flux level model can provide. Knutson et al. (2010) have recently suggested a correlation between stellar activity and the absence of a temperature inversion: in this hypothesis high-altitude absorbing species are photodissociated by the ultraviolet flux from an active star. It remains to be seen how this theory addresses the issue of temporally variable stellar activity (Shkolnik et al. 2008).

Phase curves provide complementary insights into planetary atmospheres. If a system's total (star plus planet) infrared flux varies periodically and in phase with the planet's orbit, the variation can be attributed to spatially nonuniform radiation emitted by the planet. Such measurements have the potential to constrain the planet's circulation and heat redistribution patterns. If incident stellar flux were instantaneously reradiated by the planet, the hottest region on the planet would be at the substellar point; such a phase curve is said to have zero phase offset. Nonzero phase offsets thus imply heat transport around the planet; for example, by advection of absorbed stellar energy by global winds (Showman et al. 2009) or by heating induced by atmospheric gravity waves (Watkins & Cho 2010). Interpreting phase curves can be challenging because the brightest atmospheric region will also depend on the opacity structure of the atmosphere and the wavelength at which one observes the system.

The extrasolar planet υ Andromedae b (υ And b) is the first exoplanet for which a phase curve was reported. Harrington et al. (2006, hereafter H06) used Spitzer/MIPS to measure the 24 μm system flux at five epochs over one orbit and reported a finite amplitude phase curve consistent with zero phase offset, though here we report a new analysis and interpretation of our H06 data based on a better understanding of MIPS systematics. The 8 μm observations of Cowan et al. (2007) also found variations with zero phase offset for HD 179949b. Observations of the less intensely irradiated planet HD 189733b at 8 μm and 24 μm (Knutson et al. 2007, 2009) revealed a relatively small temperature contrast between the planet's day and night sides and a 30°–40° phase offset, indicating a moderate level of eastward heat redistribution from the warm day side to the cool night side. It is important to note a possible observational bias: the first two-phase curves were sparsely sampled and are for non-transiting systems with unknown orbital inclinations. Though simulations suggest that inclination should not substantially affect a planet's observed phase offset, the flux amplitude will be directly affected by inclination (Rauscher et al. 2008); furthermore, phase curve interpretations are more ambiguous without the absolute calibration provided by a secondary eclipse (Burrows et al. 2008).

Several groups have hypothesized a connection between temperature inversions and the magnitude of a system's phase offset (e.g., Burrows et al. 2008; Fortney et al. 2008). The favored (though unproven) cause of inversions is a species residing at high altitude that absorbs optically but is transparent to infrared radiation. In this scenario, an inverted atmosphere absorbs stellar energy at lower pressures where it should quickly reradiate to space; in a non-inverted atmosphere the energy is absorbed much deeper, where it may circulate farther around the planet and cause a measurable phase offset. Showman et al. (2009) do a fair job of reproducing the HD 189733b phase curves, but they predict secondary eclipse depths for the more highly irradiated HD 209458b that do not match the observations of Knutson et al. (2008); Burrows et al. (2010) also model HD 209458b and they, too, do not match the observed eclipse depths especially well. Thus, our current understanding of the atmospheric structure and dynamics of even the best-characterized planets still appears to be incomplete.

This is the context in which we obtained the high-cadence 24 μm phase curve of υ And b described below. We introduce the υ Andromedae system and discuss our observations and data analysis in Section 2. In Section 3, we describe the planetary temperature contrast and heat redistribution implied by our analysis. We discuss possible interpretations of our results in Section 4 and conclude in Section 5.

2. OBSERVATIONS AND ANALYSIS

2.1. The υ Andromedae System

The planet υ And b was one of the earliest reported hot Jupiters (Butler et al. 1997), and the three-planet υ And system has been observed numerous times in the years since (Butler et al. 2006; Naef et al. 2004; Wittenmyer et al. 2007; McArthur et al. 2010). The host star has a spectroscopic effective temperature of 6212 ± 64 K (Santos et al. 2004) and a directly measured diameter of 1.631 ± 0.014 R_☉ (Baines et al. 2008). Spectroscopic and isochronal mass estimates generally agree on a mass of ∼1.3 M_☉ (Fuhrmann et al. 1998; Ford et al. 1999; Valenti & Fischer 2005; Takeda et al. 2007). Because this system's planets do not transit we do not know their physical sizes; however, if we assume υ And b is a typical hot Jupiter we can estimate its radius to be ∼1.3 R_J.⁸

Using a combination of radial velocity and astrometry, McArthur et al. (2010) recently determined the orbits of the second and third planets in the υ And system to be mutually inclined by 30°. They suggest this nonplanar system may result from planet–planet scattering that could also have moved the innermost planet, υ And b, into its current orbit at 0.059 AU. The small stellar reflex motion induced by υ And b precluded a direct measurement of its orbital inclination, but their preliminary numerical simulations extending 10⁵ yr suggest υ And b's inclination may lie in the range ∼20°–45° (implying a planetary mass of 2–3 M_J). Though their simulations did not fully explore the available parameter space, the inclination range McArthur et al. (2010) suggest for υ And b is broadly consistent with the constraints we place on its inclination in Section 3.

2.2. Observations

We observed the υ And system with Spitzer's MIPS 24 μm channel (Rieke et al. 2004) with observations spread across 1.2 orbits of υ And b (∼5 days) during 2009 February. The observations consist of seven brief (∼3000 s on target) observational epochs and one long, near-continuous observation ∼28 hr in length and centered at phase 0.5 (secondary eclipse for transiting systems in circular orbits, and the predicted time of flux maximum based on H06). Our integrations total 18.5 hr. Spitzer breaks up observations into blocks of time called astronomical observation requests (AORs) for instrument scheduling purposes: our short observations consist of three sequential AORs, and the long observing sequence consists of 71 AORs. Altogether our data consist of 25,488 frames, each with an integration time of 1.57 s. We also observe a flux calibrator star, HD 9712, in three two-AOR epochs, for a total of 1.3 hr of integration. The observations of υ And by H06, which we reanalyze, consist of five AORs spaced over ∼5 days.

2.3. Data Reduction

We use the basic calibrated data (BCD) files generated by version 18.14 of the MIPS data reduction pipeline (Masci et al. 2005). During MIPS observations, the instrument and spacecraft dither the target star between 14 positions on the detector (SSC 2007, Section 8.2.1.2). As noted previously (Deming et al. 2005; H06; Knutson et al. 2009), the MIPS detector response is spatially nonuniform and so we treat the observations as consisting of 14 separate time series, modeling their systematics separately in the final fit. In each frame, we measure the system flux and the position of the star on the detector by fitting a 100× supersampled model MIPS point-spread-function (PSF)⁹ generated using a 6200 K blackbody spectrum. We shift and scale the model PSF to determine the best-fitting combination of background, stellar flux, and PSF position. Using position-dependent model PSFs does not significantly change our results, so in all our photometry we use a single model PSF generated at the center of the MIPS field of view. We exclude hot pixels from the PSF fit by setting to zero the weight of any frame's pixel that is more than 5σ discrepant from the median value of that pixel for all frames taken at that particular dither position and also exclude bad pixels flagged by the MIPS reduction pipeline.

In each set of ∼40 frames, the first frame has ∼3% higher stellar flux, and the first several frames have a lower background, than the rest of the frames in the set. These effects may be related to the MIPS "first-frame" chip reset effect, though the effect we see is qualitatively different from the one described in the MIPS Instrument Handbook (SSC 2010). We exclude the first frame in each set from the remainder of our analysis, removing 708 frames. We further exclude the first three contiguous AORs (700 frames, ∼40 minutes), which are markedly discrepant from the final time series. These first data were taken soon after a thermal anneal of the MIPS detector, and during these observations we see anomalous readings in the 24 μm detector anneal current, the scan mirror temperature, and the MIPS B side temperature sensor. We do not see such anomalous readings in the rest of the 2009 data set or in the 2006 data set.

The initial extracted photometry reveals a clear sinusoidal flux variation, but with additional variability correlated with the subpixel motion of the star on the detector, as shown in Figures 1 and 2. This effect is likely from imperfect flat fielding rather than an intrapixel sensitivity variation as seen in the Infrared Array Camera's 3.6 and 4.5 μm channels (Reach et al. 2005; Charbonneau et al. 2005), because several time series with similar intrapixel locations exhibit anticorrelated flux variations. We tried creating a flat field by median-stacking all the BCD frames after masking the region containing the target star, but applying this flat field to the data before computing photometry did not reduce the amplitude of the position-correlated photometric variability. As we describe below, we suspect intermittent, low levels of scattered light interfere with our ability to construct a sufficiently accurate flat field.

**Figure 1.** Removal of pointing-correlated photometric variability. The top panel shows photometry after bringing all 14 individual time series to a common level. A large-scale sinusoidal flux variation is evident, but so is a shorter-scale "ripple" (near phase 0.5); this ripple is correlated with the motion of υ And on the MIPS detector, plotted as X and Y in the middle two panels. The bottom plot shows the final, cleaned photometry after removing this correlation. For display purposes, the full data set has been averaged over each AOR in this figure. Uncertainties on the X and Y points are typically ∼10⁻³ pixel, while the uncertainties of the photometric points (not plotted for clarity) are typically ∼1.5 × 10⁻⁴.
Download figure:
Standard image High-resolution image

**Figure 2.** Raw photometry (dots) and the best-fit model (solid lines) at each of the 14 dither positions. The measured flux varies by several percent from one position to another, as evidenced by the different scales in each panel. The downward trend evident in all panels is due to the decreasing flux from the system, shown more clearly in Figure 3. For display purposes, the data have been averaged over each AOR and we plot only the continuous-observing segment of our observations. The precision of individual points (not plotted for clarity) is ∼5 × 10⁻⁴.
Download figure:
Standard image High-resolution image

Instead, we treat this variability by removing a linear function of position from the computed photometry at each of the 14 dither positions. Fitting functions with a higher-order dependence on position, or including cross terms, does not change our final results and increases the Bayesian Information Criterion (BIC = χ² + kln N, where k is the number of free parameters and N is the number of points). We remove the systematic effects, and assess possible phase functions, by fitting the jth data point at the ith dither position (f_ij) with the following equation:

$\begin{equation} f_{ij} = (a + b \sin [\Omega _{{\rm orb}} t - \Phi _0 ]) (1 + c_i + d_i x_{ij} + e_i y_{ij}). \end{equation} \tag{ 1 }$

This equation contains three astrophysical parameters and 42 instrumental parameters that account for the nonuniform detector response. The astrophysical parameters of interest are the average system flux a and a time (t)-dependent sinusoidal phase function with known planetary orbital frequency Ω_orb but unknown amplitude b and phase offset Φ₀. The remaining parameters represent systematic effects to be removed: residual sensitivity corrections c_i and linear dependence on detector position (d_i, e_i) for each of the 14 time series.

To prevent parameters a and b from floating, we artificially set c₁ to satisfy the relation Π_i(1 + c_i) = 1. Ignoring the (d_i, e_i) factors does not change our final result for the 2009 data set but increases the BIC, indicating a poorer fit to the data. For the purposes of backward comparison, we also apply our analysis to the observations of H06. Due to the limited temporal coverage of this data set including the (d_i, e_i) does not improve the fit; therefore, we use only 17 parameters when reanalyzing this earlier data set. We otherwise apply the same data reduction steps as described above.

In this work, we consider only a sinusoidal phase function. Although we recognize that the phase curve will not be truly sinusoidal in shape, such a model is simple to work with and has a straightforward interpretation as a two-hemisphere "orange slice" model (Cowan et al. 2007; Cowan & Agol 2008, see also Section 3.2). More complicated models would necessitate fitting for the system's orbital inclination as well, which would be much more computationally intensive since there is no simple analytic expression for orange slice models of arbitrary inclination (see Cowan & Agol 2008). Because the motion-correlated flux variation only depends on relative motion (rather than on absolute detector position) we normalized the (x_ij, y_ij) in Equation (1) by subtracting the mean position in each of the 14 time series.

2.4. Calibration and Instrument Stability

Our continuous photometry reveals that the MIPS background flux changes discontinuously from one AOR to the next. Stellar photometry is not similarly affected. The background flux varies at the level of 0.1%, comparable in amplitude to the expected planetary signal. We also see these discontinuities in background flux during MIPS observations of HD 189733b (Knutson et al. 2009) and HD 209458b (unpublished; Spitzer Program ID 40280).

Because the measured background level varies with the Spitzer AORs it is extremely unlikely that this variability is of astrophysical origin. It is also unlikely that the background variability results from the calibration process because we see the same effect in both the raw and calibrated data products. We observe no correlation between the background variability and the various reported instrumental parameters, though we cannot rule out either intermittent scatter from other sources or slight changes in the detector bias. A global sensitivity drift does not seem to be the culprit because the changes in background flux are uncorrelated with stellar photometry.

Although we were unable to determine the cause of the background variations, we suspect that they are due to small changes in the amount of scattered light reaching the MIPS detector. Variable scattered light could also explain our inability to remove the motion-correlated photometric jitter with an empirical flat field. Stellar photometry was substantially more stable on short timescales than was the background in all of the extended MIPS 24 μm observations we examined, so we use this photometry in our subsequent analysis.

Our observations include a flux calibrator star to check MIPS's long-term photometric stability. MIPS 24 μm photometry is known to be stable at the 0.4% level over several years (Engelbracht et al. 2007); if it were this variable on short timescales, we would be unable to discern the expected planetary emission. We observed the K1 III star HD 9712, taken from a catalog of bright interferometric calibrator stars (Mérand et al. 2005), during six AORs. In our reduction of these data, we do not apply a correction for pixel motions and achieve a repeatability of ∼10⁻⁴. These observations, shown as the red triangles in Figure 3, imply that the MIPS 24 μm detector was stable over the course of our observations, and so we rule out detector sensitivity drifts as the source of the observed flux variations.

**Figure 3.** Phase curve of the υ And system, phased to the orbit of the innermost planet (black circles). Position-dependent sensitivity effects have been removed and for display purposes we have averaged the data over each AOR. The best-fit sinusoid (black dotted line) exhibits a phase offset of ∼80°, consistent with a planetary "hot spot" advected almost to the planet's day–night terminator. The light green circles are our reanalysis of the data from Harrington et al. (2006), and the light green dashed line shows the best-fit sinusoid to these data; the phase coherence between the 2006 and 2009 data sets is consistent with flux modulated by the innermost planet's orbit. The red triangles show our flux calibrator observations, which are consistent with a constant detector sensitivity.
Download figure:
Standard image High-resolution image

The level of MIPS photometric stability was not well known early in the Spitzer mission, so H06 used the background, attributed to smoothly varying zodiacal emission, as a calibrator to adjust the photometry (see H06, Figure 1B). As we now suspect the MIPS background variations to result from scattered light in the instrument, there is no longer any need, nor possibility, to use the zodiacal light as a calibrator. The H06 data were part of a preliminary Spitzer program to assess variability of several systems for subsequent study in programs such as ours. Its five υ And AORs were taken at separate epochs and there was no calibrator star, so it was impossible to make the assessment described above. We have reanalyzed the H06 data using our new procedure. The original phase variation still persists, as the photometric values are consistent with those plotted in H06, Figure 1A (see our Figure 3). Table 2 provides accurate phase curve parameters for that data set.

3. RESULTS

3.1. System Flux

Although our primary science result—the planetary phase curve described below—is inherently a relative measurement, our observations also allow us to measure precise, absolute 24 μm photometry for the υ And system. We measure F_ν = 0.488 Jy and 0.490 Jy for the 2009 and 2006 data sets, respectively. These fluxes differ by 0.4%, which is at the limit of the MIPS 24 μm precision; we therefore report the mean system flux as 0.489 ± 0.002 Jy. This value is significantly discrepant from the IRAS 25 μm flux of 0.73 ± 0.05 Jy (Moshir 1989), but it is consistent with a Kurucz (1979) stellar spectrum tied to optical and near-infrared photometry of υ And from the Tycho-2 (Høg et al. 2000) and Two Micron All Sky Survey (Skrutskie et al. 2006) catalogs.

3.2. A Two-hemisphere Model

As shown in Figure 3 and discussed below, we detect a flux variation consistent with the orbital period of υ And b. This measurement allows us to put tight constraints on the temperature contrast and phase offset of the planet. We interpret the observed flux variation in the context of a planet composed of two blackbody hemispheres of constant temperature—i.e., a two-wedge "orange slice" model (see Cowan et al. 2007). Sufficiently precise observations of such a bifurcated planet will reveal a flux variation with peak-to-trough amplitude:

$\begin{equation} \frac{\Delta F_P}{\langle F \rangle } = \frac{B_{\nu }(T_{P1}) - B_{\nu }(T_{P2}) }{B_{\nu }(\gamma T_{{\rm eff}})} \left(\frac{R_P}{R_*} \right)^2 \sin i \hbox{,} \end{equation} \tag{ 2 }$

which normalizes the full amplitude of the planetary flux variation ΔF_P by the mean flux from the system 〈F〉; in Equation (1) ΔF_P/〈F〉 = 2b/a. The quantity γ accounts for the star being fainter in the mid-infrared than a blackbody with temperature T_eff; at 24 μm we set γ = 0.8 based on the models of Kurucz (1979). Thus, measuring ΔF_P/〈F〉 gives the hemispheres' brightness temperature contrast relative to the stellar flux, modulo a sin i ambiguity. By assuming a planetary albedo A_B and assuming all emitted radiation is reprocessed starlight, a second (bolometric) relation obtains

$\begin{equation} \left(1 - A_B \right) \frac{R_*^2}{2 a^2} T_{{\rm eff}}^4 = T_{P1}^4 + T_{P2}^4 \hbox{.} \end{equation} \tag{ 3 }$

Thus, we implicitly assume that each hemisphere of our model planet emits a bolometric flux equal to that of a blackbody with the 24 μm brightness temperature in Equation (2). Subject to this assumption and given known or assumed values for the albedo A_B and planetary radius R_P, one can use Equations (2) and (3) to determine the hemispheric temperatures T_P1 and T_P2 at arbitrary orbital inclinations.

Setting T_P2 = 0 and T_P1 = T_P1,max gives the minimum planetary radius capable of reproducing the observed flux variation, ΔF_P/〈F〉, as a function of the orbital inclination. The computed radius will necessarily depend somewhat on the particular temperature distribution; for example, in an instantaneous reradiation model the hemisphere-integrated dayside brightness temperature is significantly higher than calculated here. Because this relation depends on inclination as (sin i)^−1/2, measuring ΔF_P/〈F〉 gives an upper limit to the planet's surface gravity:

$\begin{equation} g \le \frac{G (m \sin i)}{R_*^2} \frac{B_{\nu }(T_{P1,{\rm max}})}{B_{\nu }(\gamma T_{{\rm eff}})} \left(\frac{\Delta F_P}{\langle F \rangle } \right)^{-1} \hbox{.} \end{equation} \tag{ 4 }$

3.3. Model Fits

We determine the best-fit parameters using Powell's (1964) method for multivariate minimization (the SciPy function optimize.fmin_powell) and assess their uncertainties and correlations with a Metropolis–Hastings, Markov Chain Monte Carlo analysis (MCMC; see Press et al. 2007, Section 15.8). Table 1 (for the 2009 data) and Table 2 (for the 2006 data) report these results. Table 3 lists the astrophysical parameters of interest. We list the χ² and BIC values for the fits in Table 4. Parameter uncertainties are estimated from distributions generated using the kernel density method (KDE, implemented using the SciPy function stats.gaussian_kde) by determining the parameter values with equal KDE frequency that enclose 68% of the distribution.

Table 1. Phase Curve and Calibration Parameters from Equation (1) (2009 Data)^a

Parameter	Value and Uncertainty
a^b	0.488436 ± 0.000011 Jy
b	0.000317 ± 0.000017 Jy
Φ^c₀	845 ± 23
c₁	0.005367 ± 0.000054
c₂	0.012926 ± 0.000055
c₃	0.001180 ± 0.000055
c₄	0.009498 ± 0.000055
c₅	0.001548 ± 0.000054
c₆	0.014334 ± 0.000027
c₇	−0.006333 ± 0.000054
c₈	−0.004193 ± 0.000054
c₉	−0.002129 ± 0.000054
c₁₀	−0.007009 ± 0.000055
c₁₁	−0.000549 ± 0.000055
c₁₂	−0.009458 ± 0.000054
c₁₃	−0.002757 ± 0.000054
c₁₄	−0.012002 ± 0.000054
d₁	−0.00429 ± 0.00072
d₂	0.00146 ± 0.00070
d₃	−0.00214 ± 0.00071
d₄	−0.00428 ± 0.00070
d₅	−0.00520 ± 0.00071
d₆	−0.00177 ± 0.00070
d₇	−0.00100 ± 0.00078
d₈	−0.00049 ± 0.00072
d₉	0.00196 ± 0.00066
d₁₀	−0.00252 ± 0.00067
d₁₁	−0.00029 ± 0.00066
d₁₂	−0.00409 ± 0.00067
d₁₃	0.00259 ± 0.00068
d₁₄	−0.00152 ± 0.00068
e₁	0.00066 ± 0.00080
e₂	0.00031 ± 0.00079
e₃	0.00211 ± 0.00040
e₄	0.00218 ± 0.00078
e₅	0.00409 ± 0.00080
e₆	0.00009 ± 0.00080
e₇	0.00211 ± 0.00081
e₈	−0.00259 ± 0.00076
e₉	0.00628 ± 0.00076
e₁₀	−0.00284 ± 0.00039
e₁₁	−0.00168 ± 0.00075
e₁₂	−0.0029 ± 0.0011
e₁₃	−0.00061 ± 0.00076
e₁₄	0.00025 ± 0.00074

Notes. ^aErrors quoted are the 68.3% confidence limits. The first three parameters are the mean system flux a, the phase curve half-amplitude b, and the phase offset Φ₀. The other coefficients represent fits to the nonuniform detector response at each of the 14 MIPS dither positions. The c_i are dimensionless, and the d_i and e_i are in units of pixel⁻¹. ^bThe absolute accuracy of the mean system flux a is 0.4% as discussed by Engelbracht et al. (2007). ^cΦ₀ is measured relative to our computed ephemeris, as discussed in Section 3.5.

Download table as: ASCII Typeset image

Table 2. Phase Curve and Calibration Parameters from Equation (1) (2006 Data)^a

Parameter	Value and Uncertainty
a^b	0.490250 ± 0.000039 Jy
b	0.000221 ± 0.000053 Jy
Φ₀^c	57° ± 21°
c₁	0.00284 ± 0.00030
c₂	0.01423 ± 0.00028
c₃	0.00112 ± 0.00030
c₄	0.01094 ± 0.00028
c₅	0.00235 ± 0.00029
c₆	0.01574 ± 0.00028
c₇	−0.00757 ± 0.00027
c₈	−0.00596 ± 0.00030
c₉	−0.00086 ± 0.00030
c₁₀	−0.00755 ± 0.00029
c₁₁	−0.00097 ± 0.00029
c₁₂	−0.00945 ± 0.00029
c₁₃	−0.00299 ± 0.00027
c₁₄	−0.01138 ± 0.00027

Notes. ^aErrors quoted are the 68.3% confidence limits. The first three parameters are the mean system flux a, the phase curve half-amplitude b, and the phase offset Φ₀. The other coefficients represent fits to the nonuniform detector response at each of the 14 MIPS dither positions. The c_i are dimensionless, and the d_i and e_i are in units of pixel⁻¹. ^bThe absolute accuracy of the mean system flux a is 0.4% as discussed by Engelbracht et al. (2007). ^cΦ₀ is measured relative to our computed ephemeris, as discussed in Section 3.5.

Download table as: ASCII Typeset image

Table 3. Astrophysical Parameters of Interest

Parameter	2006	2009
ΔF/F	0.00090 ± 0.00022	0.001300 ± 0.000074
Φ₀	57° ± 21°	845 ± 23

Download table as: ASCII Typeset image

Table 4. Goodness-of-fit Statistics

Parameter	2006	2009
χ²	814.1	38422.7
k	17	45
N	838	23884
BIC	928.6	38876.4

Download table as: ASCII Typeset image

MCMC analysis evolves an initial set of parameters in a way that is ultimately representative of their underlying probability distributions. For our Markov chain, we choose a step size to give approximately a 30% step acceptance rate. To adequately sample the full parameter space, we found it necessary to run the Markov chains longer for the 2006 data set than for the 2009 data set. For the 2009 data set, we first ran the chain for 10⁶ burn-in steps and discarded these; we then ran the chain for 2 × 10⁷ steps, saving every 1000th step. For the 2006 data set, our procedure is the same but the burn-in phase lasted for 10⁷ steps and the chain was then run for 5 × 10⁷ steps, saving every 1000th step. We inspected correlation plots for all possible parameter pairs in both analyses to ensure adequate coverage of phase space and to assess parameter correlations.

All the one-dimensional parameter distributions are unimodal and approximately Gaussian in shape. We see some correlations between d_i and e_i, which are unsurprising given the degree of correlation between the X and Y components of motion as shown in the middle panels of Figure 1. More surprising is a correlation between the mean system flux, a, and the phase offset, Φ₀, as shown in Figure 4. Using a simulated data set with white noise, we confirmed that when forcing a fit to a sinusoid of known period, a slightly higher mean value be counteracted by a slightly lower phase offset; however, we observe the opposite correlation. In any case no significant correlation is apparent between the phase curve amplitude ΔF_P/〈F〉 and the phase offset, which are the primary quantities of interest for our analysis. The best χ² from the MCMC is consistent to within a small fraction of the uncertainties with the optimizer values.

**Figure 4.** Two-dimensional confidence intervals on the astrophysical parameters of interest in Equation (1): the relative phase curve amplitude Δ*F_P*/F, the absolute mean system flux a, the phase curve half-amplitude b, and the phase offset Φ₀. The solid lines are the contours that enclose 68.3%, 95.4%, and 99.7% of the parameter space from the 2009 data set. The "×" in each panel marks the best-fit parameter listed in Tables 1 and 3.
Download figure:
Standard image High-resolution image

The model in Equation (1) provides a good fit to both the astrophysical flux modulation and the instrumental flux variations at each of the 14 dither positions, as shown in Figure 2. We plot the photometry after removal of the systematic effects in Figure 3, along with the best-fitting sinusoidal phase curve, for both the 2006 and 2009 data sets. Both data sets appear to vary approximately in phase; this coherence is a strong argument that the flux modulation we see is due to the planet υ And b. The goodness-of-fit statistics (χ² and BIC) for both data sets are listed in Table 4. The high χ² for the 2009 data set probably results from a somewhat non-sinusoidal phase curve and from the residual systematics apparent in Figure 3.

3.4. Phase Curve Amplitude

We measure values of ΔF_P/〈F〉 of ΔF_P/〈F〉 = 0.001300 ± 0.000074 for the 2009 data and 0.00090 ± 0.00022 for our reanalysis of the 2006 data, as shown in Table 3. Using the absolute calibration from Section 3.1, we find absolute peak-to-trough phase curve amplitudes of 0.636 ± 0.036 mJy and 0.44 ± 0.11 mJy for the 2009 and 2006 data sets, respectively. Thus, the detection of the phase curve amplitude at both epochs is statistically significant at the >4σ level and is substantially smaller than reported by H06 (due to the calibration issues discussed above). The two epochs' phase curve amplitudes are consistent at the 1.7σ level; thus there is no evidence that the planetary emission exhibits inter-epoch variability. The lack of variability is consistent with the recent results of Agol et al. (2010), who set an upper limit of 2.7% on HD 189733b's dayside flux variations.

As υ And b's orbit is inclined toward face-on, a greater intrinsic temperature contrast is required to generate the observed flux variation. Using Equation (2) and our measurement of the phase curve amplitude, we determine the expected day/night contrast ratio and plot it in the upper panel of Figure 5. We also plot the upper limits on the day/night contrast assuming planetary radii of 1.3 (1.8) R_J; the implication is that the planet's orbital inclination angle is likely ≳28° (14°). These limits complement the preliminary limits on υ And b's orbital inclination from the stability modeling of McArthur et al. (2010), which suggest i ∼ 20°–45°.

**Figure 5.** Day–night temperature contrast as a function of orbital inclination angle, assuming a planet with zero albedo. (a) Our measurement of the phase curve amplitude and Equation (2) constrain the allowed day/night contrast to lie between the dot-dashed 3σ limits. The maximum allowable contrast for a planet of 1.3 (1.8) *R_J* is shown as the dashed (dotted) lines, suggesting a lower limit on the inclination angle of 25° (15°). (b) Considering Equation (3) and assuming a radius of 1.3 *R_J* allow us to determine the temperature of both hemispheres as a function of inclination (solid lines, with dot-dashed 3σ limits). This radius implies a planetary temperature contrast of ≳900 K. The cooler hemisphere is more sensitive to changes in planetary radius, though both temperatures increase as radius increases.
Download figure:
Standard image High-resolution image

Invoking Equation (3) allows us to determine the brightness temperatures of the planetary hemispheres in our model at each inclination angle. We assume zero albedo (see Rowe et al. 2008) and a planetary radius of 1.3 R_J and plot the hemispheres' temperatures and 3σ limits in the lower panel of Figure 5. This sets a lower bound to the temperature contrast between the two hemispheres to be T_P1 − T_P2 ≳ 900 K. The hotter hemisphere's temperature remains in the range ∼1700–1900 K as we vary the radius from 1.0 R_J to 1.8 R_J; however, larger radii result in higher temperatures for the cooler hemisphere to maintain the measured flux and thus decrease the temperature contrast.

Using Equation (4), we find that υ And b's surface gravity is <2100 cm s⁻² with 3σ confidence: this result is independent of assumptions about the planet's radius or orbital inclination. For a hot hemisphere temperature of ∼ 1800 K (see Figure 5), this limit on the surface gravity implies an atmospheric scale height >300 km. In Figure 6, we plot the allowed regions of mass–radius parameter space against the known population of transiting extrasolar planets. Thus, our measurements suggest that υ And b has a lower surface gravity than Jupiter, HD 189733b, and a number of other transiting extrasolar planets. This result demonstrates that υ And b is indeed a gaseous Jovian planet, but we cannot determine whether it is a highly inflated planet or whether it is dominated by a sizeable rocky core.

**Figure 6.** Mass/radius constraints for υ And b from Equation (4). The shaded area is the portion of mass–radius space excluded at the 3σ level. Points are known transiting extrasolar planets; measurement uncertainties have been omitted for clarity. Extrasolar planets that have been observed with MIPS are indicated by a "+" (HD 189733b) and a "×" (HD 209458b), while Jupiter and Saturn are marked with a "J" and "S," respectively. The dot-dashed lines represent lines of constant surface gravity with g = 10³ cm s⁻² and 3.2 × 10³ cm s⁻². Our model and phase curve measurement constrain g < 2100 cm s⁻² at the 3σ level.
Download figure:
Standard image High-resolution image

3.5. Phase Offset

Because υ And b does not transit its host star we know its orbital ephemeris less precisely than we do for transiting planets. We used the Systemic Console¹⁰ (Meschiari et al. 2009) to reanalyze the published radial velocity data of υ And (Butler et al. 1997, 2006; Naef et al. 2004; Wittenmyer et al. 2007; McArthur et al. 2010) using Systemic's Levenberg–Marquardt algorithm and ignoring system stability constraints. We obtain orbital parameters consistent with those of McArthur et al. (2010). By providing the covariances between the various fit parameters, this reanalysis allows a substantially more precise estimate of the planetary ephemeris than is available from the literature. We compute a time of zero relative radial velocity (phase 0.5 or secondary eclipse in a circular transiting system) of JD = 2,454,868.78 ± 0.07 (1σ), which corresponds to an uncertainty of ∼6° in determining the phase offset.

Assuming a two-hemisphere model, we find a phase offset of 84 fdg 5 ± 2 fdg 3 relative to the computed ephemeris. This uncertainty may be an underestimate since we artificially constrain our phase curve to be sinusoidal. Knutson et al. (2009) discuss the artificially low uncertainties obtained from fitting to an arbitrarily chosen model, though here we have substantially broader phase coverage than was available in that study. The system flux reaches a maximum before phase 0.5, indicating that the brighter hemisphere is offset to the east of the substellar point, as observed for HD 189733b (Knutson et al. 2007, 2009; Agol et al. 2010).

This large phase offset is strikingly different from the near-zero phase offset reported by H06; this difference is due to the choice of system calibration as discussed above. From our reanalysis of the 2006 data, we find a phase offset of 57° ± 21° relative to our ephemeris. The phase offsets at the two epochs are consistent at the 1.3σ level, and thus there is no evidence for inter-epoch variability in the phase offset.

4. DISCUSSION

The most striking result of our analysis is the large phase offset evident in the light curve. The direction of the phase offset is broadly consistent with the prediction by many circulation models of a large-scale, high-velocity, eastward-flowing jet on hot Jupiters (Showman & Guillot 2002; Cooper & Showman 2005; Cho et al. 2008; Showman et al. 2009; Burrows et al. 2010; Rauscher & Menou 2010; Thrastarson & Cho 2010) and as seen on HD 189733b (Knutson et al. 2007, 2009). However, the magnitude of the phase offset is far larger than is predicted at the low pressures characteristic of the 24 μm photosphere (Showman et al. 2009).

A partial explanation for such a large phase offset could be that υ And b's atmosphere is substantially transparent to the incident stellar flux, with the result that the insolation is deposited at sufficient depth for substantial advection to occur. In this case, the planet might lack an atmospheric temperature inversion. Indeed, the large phase offset is more consistent with expectations for a planet lacking a temperature inversion than with expectations for a planet with an inverted atmosphere (cf. HD 189733b versus HD 209458b; Showman et al. 2009). Cooper & Showman (2005, 2006) predicted maximum hemisphere-averaged temperatures to be offset by ≳60° for an irradiated planet with no high-altitude absorbers, but it is unlikely that this would translate into high phase offsets at 24 μm due to the high altitude of the photosphere at this wavelength. The absence of a temperature inversion is at odds with the recent proposal by Knutson et al. (2010) that due to the lack of activity apparent in optical spectra of υ Andone would expect υ And b to have a temperature inversion. Thrastarson & Cho (2010) also show clear shifted hot regions at higher pressure; the deep vortices produced by their simulations are several hundreds of degrees hotter than the surroundings.

Setting aside the phase offset for a moment, the phase curve amplitude we measure could be directly interpreted in the context of the models of Showman et al. (2009) as a typical phase curve for a planet with orbital inclination of 50°–60°. However, simulations with greater phase offsets tend to show lower phase curve amplitudes (e.g., Showman et al. 2009; Burrows et al. 2010); thus it is difficult to simultaneously explain both amplitude and offset and it is difficult to reconcile current theory with our observations.

Alternatively, the large phase offset we see could represent reradiation of thermal energy deposited in shock fronts in the planet's atmosphere. Many simulations predict supersonic equatorial jets on hot Jupiters that carry substantial kinetic energy; shocks could manifest themselves where the jet transitions to subsonic speeds. Rauscher & Menou (2010) note that their simulations, as well as those of Showman et al. (2009), exhibit structures which could be interpreted as shocks—however, these models do not explicitly treat shock physics. Dobbs-Dixon et al. (2010) observe similar features in their models using an additional artificial viscosity factor to simulate shock behavior, and note that in high-altitude, high-velocity regions the energy carried by kinetic energy becomes comparable to the enthalpic energy. Watkins & Cho (2010) have recently suggested that gravity waves in the atmosphere of a hot Jupiter can also heat the planet's upper atmosphere. It is unclear whether either shocks or gravity waves can deposit sufficient energy in the 24 μm photosphere to cause the large phase offset we see. We look forward to further research into these topics to determine whether these or other phenomena can explain our observations.

We must also consider the possibility that the periodic flux modulation we see is intrinsic to the star rather than emanating from the planet. Shkolnik et al. (2005, 2008) report evidence for intermittently periodic stellar activity in this system correlated with υ And b's orbital period. They interpret this periodicity as a possible magnetospheric star–planet interaction, but they detect this periodicity at only some of their observational epochs. The consistency between our analysis of the 2006 and 2009 MIPS data suggests that we are not seeing such a transient phenomenon, although analysis of additional activity measurements (e.g., from the Keck Observatory Archive) spanning the time of our Spitzer observations would help solidify this claim. Shkolnik et al. (2005) also report variations consistent with zero phase offset, which disagree with our observed phase offset. In addition, 24 μm stellar variability at the level we observe would imply much greater variability at optical wavelengths: this is not observed (Henry et al. 2000). Finally, as discussed by H06, energy considerations indicate that such an intense star–planet interaction would cause the planet to spiral into the star in ≲10⁷ yr. Thus, the 24 μm flux variations are likely to be of planetary origin.

Although a general framework exists within which to interpret observations of exoplanetary atmospheres, our understanding is still extremely limited. Madhusudhan & Seager (2009) demonstrate that some planets with claimed inversions can also be fit by non-inverted atmospheric models due to much greater number of free model parameters versus the limited number of observational constraints. In addition, there may be reason to question the reliability of some of the circulation models currently in use. Thrastarson & Cho (2010), in their extensive exploration of initial condition parameter space, have recently shown the extreme susceptibility of at least one circulation model to minute variations in initial conditions, resulting in substantial variability in final "steady-state" temperature contrasts and phase offsets. This chaotic behavior calls into question the ability of at least some models to make accurate, qualitative predictions about any of the quantities we are interested in. Other models (e.g., Dobbs-Dixon et al. 2010) exclude planets' polar regions from their simulations, where large vortices are often seen to form (Cho et al. 2003, 2008; Rauscher et al. 2008). While the polar regions present a small cross-sectional area in transiting systems, for non-transiting systems such as υ And b polar emission will constitute a larger component of the observed system flux. A fully consistent three-dimensional circulation geometry is essential for comparison to our observations. Nevertheless, it is important to remember that general circulation models have trouble predicting global weather patterns even for relatively well-studied solar system planets; thus, a comprehensive, quantitative understanding of extrasolar planetary atmospheric dynamics will likely remain elusive for some time to come.

5. CONCLUSIONS AND FUTURE WORK

We have described a new 24 μm phase curve, which we interpret as being due to emission from the planet υ Andromedae b modulated by the planet's orbit. Using a simple two-hemisphere model, we determine the peak-to-trough phase curve amplitude to be 0.001300 ± 0.000074. This result suggests an average "hot side" temperature of ∼1800 K; for an average-sized hot Jupiter (1.3 R_J), this implies a hemisphere-averaged planetary temperature contrast of ∼900 K and an orbital inclination i ≳ 28°.

We find a phase offset of 84° ± 2° ± 6°, where we break the uncertainty into the error relative to our ephemeris and the error in our ephemeris, respectively. Such a large phase offset is difficult to reconcile with most current models. The phase curve is consistent with the hotter hemisphere offset to the east, as previously observed for HD 189733b (Knutson et al. 2007, 2009).

We reanalyze our earlier (H06) observations of this system and find a phase curve amplitude of 0.00090 ± 0.00022 and a phase offset of 57° ± 21° ± 6°. There is no evidence for inter-epoch variability in the planetary phase curve. This is primarily due to the large uncertainties from the 2006 data set and demonstrates the difficulty in measuring such variability with sparsely sampled phase curves.

There are substantial challenges in interpreting a phase curve observed at only a single wavelength due to the degeneracies between planetary radius, orbital inclination, and atmospheric composition and structure (Burrows et al. 2008). Some of these difficulties could be mitigated with phase curve measurements at additional wavelengths—ideally from space (i.e., warm Spitzer) but also potentially from ground-based near-infrared observations (e.g., Barnes et al. 2010). For example, commensurate phase curve amplitudes at both 3.6 μm and 4.5 μm would suggest that υ And b does indeed lack a temperature inversion. On the other hand, differing phase curve amplitudes at these wavelengths could suggest an inversion and be more difficult to reconcile with the large phase offset at 24 μm. Phase offsets in these or other scenarios would also depend on the particular atmospheric temperature structure of the planet. Phase curves at 24 μm currently exist for only two planets, υ And b and HD 189733b (though unpublished data exist for HD 209458b). Cool, inversionless HD 189733b is probably the best-characterized extrasolar planet, and additional phase curves are already being observed for this object at 3.6 μm and 4.5 μm. Phase curves of additional planets at multiple wavelengths are essential to ensure that our evolving views of the atmospheres of hot Jupiters are not biased by unbalanced data sets.

Whatever the cause of the substantial energy transport implied by the large phase offset we measure, if this phenomenon occurs in other (transiting) systems there are important implications for transmission spectroscopy. Line-of-sight effects cause optical transmission spectra to probe pressures comparable to those probed by mid-infrared emission; temperatures near the planetary terminator of ∼1800 K (as we observe) should easily be detectable with ground- or space-based spectra. Sing et al. (2008) and Lecavelier des Etangs et al. (2008) infer a terminator temperature on HD 209458b of 2200 ± 260 K at 33 ± 5 mbar—roughly at the expected 24 μm photosphere—which is not dissimilar from the hot terminator-centered hemisphere we observe on υ And b. Further phase curve and transmission spectra of additional systems are needed to determine whether this hot, high-altitude terminator measurement results from a mechanism similar to what we observe on υ And b.

We note that υ And b is too bright to be observed photometrically at shorter wavelengths with the planned James Webb Space Telescope (JWST; though spectroscopy may be feasible). Given JWST's magnitude limits and the expected high demand for its observing time, the community should consider a dedicated space-based mid-infrared photometry and spectroscopy mission (Vasisht et al. 2008). Such a mission would allow uninterrupted long-term monitoring of nearby hot Jupiter systems. This would provide high-precision measurements of these systems' thermal emission and energy distributions and possibly provide the first definitive evidence of dynamical meteorological processes—weather—on extrasolar planets.

We thank J. Colbert, C. Engelbracht, and G. Rieke for help in interpreting MIPS systematics, T. Loredo for helpful discussions of statistics, and our anonymous referee for a timely report. This work is based in part on observations made with the Spitzer Space Telescope, which is operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA. Support for this work was provided by NASA through an award issued by JPL/Caltech. We received free software and services from SciPy, Matplotlib, and the Python Programming Language. This research made use of Tiny Tim/Spitzer, developed by John Krist for the Spitzer Science Center; the Center is managed by the California Institute of Technology under a contract with NASA. Part of this work was performed while in residence at the Kavli Institute for Theoretical Physics, funded by the NSF through grant number PHY05-51164.

A NEW 24 μm PHASE CURVE FOR υ ANDROMEDAE b

Article metrics

Permissions

Author e-mails

Author affiliations

Dates

ABSTRACT

1. INTRODUCTION