Mass Loss from the Exoplanet WASP-12b Inferred from $\textit{Spitzer}$ Phase Curves

TL;DR

This study analyzes Spitzer phase curves of exoplanet WASP-12b, revealing infrared signatures of mass loss, possibly due to atmospheric escape or accretion, with implications for its orbital decay and stellar irradiation effects.

Contribution

First detection of infrared gas emission indicating mass loss from WASP-12b using phase curve analysis, highlighting unique atmospheric escape processes in ultra-hot Jupiters.

Findings

Infrared emission from stripped gas detected at 4.5 μm

Evidence of phase offset variability at 3.6 μm

Potential link between mass loss and orbital decay

Abstract

The exoplanet WASP-12b is the prototype for the emerging class of ultra-hot, Jupiter-mass exoplanets. Past models have predicted---and near ultra-violet observations have shown---that this planet is losing mass. We present an analysis of two sets of 3.6 $\mu$m and 4.5 $\mu$m $\textit{Spitzer}$ phase curve observations of the system which show clear evidence of infrared radiation from gas stripped from the planet, and the gas appears to be flowing directly toward or away from the host star. This accretion signature is only seen at 4.5 $\mu$m, not at 3.6 $\mu$m, which is indicative either of CO emission at the longer wavelength or blackbody emission from cool, $\lesssim$ 600 K gas. It is unclear why WASP-12b is the only ultra-hot Jupiter to exhibit this mass loss signature, but perhaps WASP-12b's orbit is decaying as some have claimed, while the orbits of other exoplanets may be more stable; alternatively, the high energy irradiation from WASP-12A may be stronger than the other host stars. We also find evidence for phase offset variability at the level of $6.4\sigma$ ($46.2^{\circ}$) at 3.6 $\mu$m.

Tables

Table 1. Table 1: Key fiducial light curve parameters

		1st Order	2nd Order
		Phase Offset†	Phase Offset†	$F_{\mathrm{day}}/F_{\ast }$ ‡
Data Set	$R_p/R_{\ast }$ ‡	(degrees)	(degrees)	(ppm)
2010, 3.6 $\mu$m	$0.11642\pm 0.00063$	$-32.6\pm 6.2$	—	$3870\pm 130$
2013, 3.6 $\mu$m	$0.11327\pm 0.00068$	$13.6\pm 3.8$	—	$3840\pm 120$
2010, 4.5 $\mu$m	$0.10656\pm 0.00085$	$-9.5\pm 2.3$	$94.7\pm 1.6$	$4360\pm 140$
2013, 4.5 $\mu$m	$0.1049\pm 0.0010$	$-19.1\pm 3.9$	$93.2\pm 1.9$	$3920\pm 150$

Table 2. Table A1: A summary of all the priors used in the three independent analyses. Uniform priors were used where there are inequalities below, Gaussian priors were used where uncertainties are indicated, variables were fixed where only a value is indicated, and parameters were unconstrained where Free is written.

	Zhang PLD	Bell SPCA	Cubillos POET
$t_0$ (BMJD)		56176.16825800 $\pm$ 0.00007765	Free
$R_p/R_{\ast }$	$>0$		$>0$
$a/R_{\ast }$	3.039	3.039 $\pm$ 0.034	3.039 $\pm$ 0.034
$i$ (degrees)	83.37	83.37 $\pm$ 0.68	83.37 $\pm$ 0.68 (fitted $\cos i$)
$P$ (days)	1.09142245	1.0914203	1.0914203
$F_p/F_{\ast }$	Free		$>0$
$C_1$	Free	Positive Phasecurve	Positive Phasecurve
$D_1$	Free	Positive Phasecurve	Positive Phasecurve
$C_2$	Free	Positive Phasecurve	Positive Phasecurve
$D_2$	Free	Positive Phasecurve	Positive Phasecurve
$R_{p,2}/R_{\ast }$	$>0$		$>0$
${\sigma }_F/F_{\ast }$ (white noise)		$>0$	Free
Limb Darkening	Sing 2010 Model	;	;
$e$	0	0	$t_{\mathrm{eclipse}}$: Free; $t_{14,\mathrm{eclipse}}>0$; $t_{12,\mathrm{eclipse}}>0$; $t_{34}=t_{12}$
Instrumental Variables	Free (PLD coefficients, Slope in time)	Free (Polynomial coefficients, Slope in time for 2013 3.6 $\mu$m)	Free (Slope in time for 2013 4.5 $\mu$m)

Equations

$$F_{\mathrm{model}}(t)=A(t)\times \tilde{D}(t)$$

$$A(t)=F_{\ast }(t)+F_p(t),$$

$$F_{p}(t)=F_{\rm day}\Phi\big{(}\psi(t)\big{)},$$

$$\Phi_{\rm 1}(\psi)=1+C_{1}\bigg{(}\cos(\psi)-1\bigg{)}+D_{1}\sin(\psi),$$

$$\Phi_{\rm 2}(\psi)=\Phi_{\rm 1}(\psi)+C_{2}\bigg{(}\cos(2\psi)-1\bigg{)}+D_{2}\sin(2\psi),$$

$${\Phi }_{1,\mathrm{ellipsoid}}(\psi )=S(\psi ){\Phi }_1(\psi ),$$

$$\displaystyle S(\psi)=\bigg{[}$$

$$\displaystyle+\bigg{(}\frac{R_{p,2}}{R_{p}}\bigg{)}^{2}\cos^{2}(i)\bigg{]}^{1/2},$$

$$C_{\mathrm{corr}}(\lambda )=1+g(\beta ,\lambda ){\alpha }_{\mathrm{comp}}(\lambda ),$$

$$R_{p,\mathrm{corr}}(\lambda )=\sqrt{C_{\mathrm{corr}}(\lambda )}{R}_{p,\mathrm{meas}}(\lambda ),$$

$$F_{\mathrm{day},\mathrm{corr}}(\lambda )=C_{\mathrm{corr}}(\lambda )F_{\mathrm{day},\mathrm{meas}}(\lambda ),$$

$$T_{b}(\lambda)=\frac{hc}{\lambda k_{B}}\bigg{[}\ln\bigg{(}1+\frac{\exp(hc/\lambda k_{B}T_{*,\rm b})-1}{\psi(\lambda)}\bigg{)}\bigg{]}^{-1},$$

$${\psi }_1=-(180/\pi )\arctan (D_1/C_1)$$

$${\psi }_2=180-0.5(180/\pi )\arctan (D_2/C_2).$$

$$\Omega (r,\theta ,\varphi )=\frac{1}{r}$$

$$+\frac{q+1}{2}{r}^2{\sin }^2\theta ,$$

$$\mathrm{BIC}=-2\ln (L)+N_{\mathrm{par}}\ln (N_{\mathrm{dat}}),$$

$$\ln (L)=-\frac{{\chi }^2}{2}-N_{\mathrm{dat}}\ln ({\sigma }_F)-\frac{N_{\mathrm{dat}}}{2}\ln (2\pi )$$

$$\chi^{2}=\frac{\sum_{i}\big{(}F_{{\rm obs},i}-F_{{\rm model},i}\big{)}^{2}}{\sigma_{F}^{2}}$$

Full text

††pubyear: 2019††pagerange: Mass Loss from the Exoplanet WASP-12b Inferred from Spitzer Phase Curves–Appendix G: Model Selection

1 Introduction

The exoplanet WASP-12b (Hebb et al., 2009) is one of the hottest planets known to date and, as a result of its exceedingly tight orbit and inflated radius ($a/R_{*}=3.039$, $R_{p}=1.900~{}R_{J}$; Collins et al., 2017), it is one of the best-studied exoplanets. WASP-12b is also the archetype of an emerging class of exoplanets called ultra-hot Jupiters (UHJs). Planets in this regime are so strongly irradiated by their host star that many of the molecules (e.g., H2 & H2O) in their dayside atmospheres thermally dissociate (Bell et al., 2017; Bell & Cowan, 2018; Arcangeli et al., 2018; Kreidberg et al., 2018; Lothringer et al., 2018; Mansfield et al., 2018; Parmentier et al., 2018) and may recombine nearer the nightside (Bell & Cowan, 2018; Komacek & Tan, 2018; Parmentier et al., 2018). UHJs also bear some similarities to cataclysmic variable star (CV) systems and may undergo significant tidal distortion and mass loss, depending on the specifics of the star-planet system (e.g., Bisikalo et al., 2013a; Burton et al., 2014).

While tidal distortion is expected for WASP-12b, a 2010 Spitzer Infrared Array Camera (IRAC) thermal phase curve observation of WASP-12b at 4.5 $\mu$m demonstrated second order sinusoidal variations (with two maxima per planetary orbit) that were far greater than predicted (Cowan et al., 2012). The substellar axis would have to be $1.8$ times as long as the dawn–dusk and polar axes if the observed variations were entirely due to the tidally distorted shape of the planet. Additionally, no evidence of these second order sinusoidal variations was found in a Spitzer/IRAC 3.6 $\mu$m phase curve also taken in 2010 (Cowan et al., 2012).

In this paper, we combine a new set of phase curves taken in 2013 with a reanalysis of the data from 2010 to determine the source of the unusually strong second order sinusoidal variations at 4.5 $\mu$m reported by Cowan et al. (2012). The observations are described in Section 2. Our three astrophysical models are described in Section 3.1, and our three independent reduction and decorrelation methods are described in Section 3.2. Results and their physical implications are presented in Section 4, and Section 5 contains the discussion and conclusion.

2 Observations

We combine two sets of two-channel (3.6 $\mu$m and 4.5 $\mu$m) Spitzer/IRAC observations taken in 2010 (PID 70060, PI Machalek) and 2013 (PID 90186, PI Todorov), all during the Post-Cryogenic Spitzer Mission. In all four phase curves, the system was observed nearly-continuously for $\sim$33 hours (breaking only once or twice to repoint the telescope), beginning shortly before one secondary eclipse and ending shortly after the subsequent secondary eclipse. The reduced and detrended observations are shown in Figure 1.

For both data sets, the sub-array mode was used with 2 s exposures which produced data cubes of 64 images with $32\times 32$ pixel ($39^{\prime\prime}\times 39^{\prime\prime}$) dimensions. The 2010 observations were divided into 2 Astronomical Observation Requests (AORs) with a total of 902 data cubes (57 728 exposures), while the 2013 observations were divided into 3 AORs with a total of 909 data cubes (58 176 exposures). The 2010 full-phase observations were published by Cowan et al. (2012), the eclipse timings from the 2013 observations were published by Patra et al. (2017), and some derived parameters from all four phase curves were published as part of a broad comparison between different planets (Zhang et al., 2018).

Past observations of WASP-12 show a nearby M-dwarf binary system WASP-12B,C 1.*′′*06 away from WASP-12A (Bergfors et al., 2011; Crossfield et al., 2012b; Bechter et al., 2014). As this binary system lies too close to WASP-12A to be resolved by Spitzer, we correct for blended light after analyzing the light curves, following past work (Stevenson et al., 2014a); see Appendix A for more details.

3 Light Curve Analysis

3.1 Astrophysical Models

We model the observations as

[TABLE]

where $\tilde{D}(t)$ is the normalized detector model; see Section 3.2 for details on the specific models used which consist of both parametric (2D polynomials and pixel level decorrelation) and non-parametric models (bilinear interpolated subpixel sensitivity mapping). The astrophysical model is

[TABLE]

where $F_{*}$ is the flux from the host star (assumed to be constant except during transits) and $F_{p}$ is the planetary flux. Transits and eclipses are modelled using batman (Kreidberg, 2015), assuming a quadratic limb-darkening model for the host star and a uniform disk for the planet. The planetary flux is modelled as

[TABLE]

where $F_{\rm day}$ is the instantaneous eclipse depth at phase 0.5 (assumed to be constant over each $\sim$33 hour phase curve), $\Phi$ describes the phase variations, and the orbital phase with respect to eclipse is $\psi(t)=2\pi(t-t_{e})/P$, where $t_{e}$ is the time of eclipse and $P$ is the planet’s orbital period.

We consider three different models for the astrophysical phase variations in the lightcurve. The simplest astrophysical model we consider is a first order sinusoid

[TABLE]

and we also consider a second order sinusoid

[TABLE]

where $C_{1}$, $D_{1}$, $C_{2}$, and $D_{2}$ are all constants. If the previously reported double peaked phase curve (Cowan et al., 2012) is astrophysical in nature, one potential interpretation is that some/all of the power in the second order sinusoidal variations is from tidal distortion of the planet. We model this scenario with

[TABLE]

where $S(\psi)$ describes the projected area of an ellipsoid as it rotates. Rather than model a triaxial ellipsoid, we constrain the polar and dawn–dusk axes to share the same length since rotational deformation is expected to be negligible compared to tidal deformation (Leconte et al., 2011a). To find the deviations in the projected area of this biaxial ellipsoid, we adapt an equation from past work (Leconte et al., 2011b),

[TABLE]

where $i$ is the orbital inclination, $R_{p}$ is the planetary radius along the polar and dawn-dusk axes (the two axes observed during transit and eclipse if $i=90^{\circ}$), and $R_{p,2}$ is the planetary radius along the line connecting the planet and star (the sub-stellar axis).

3.2 Decorrelation Procedures

To ensure our results are robust and independent of the methods used, we perform three independent reductions and analyses following previously employed methods (Zhang et al., 2018; Dang et al., 2018; Cubillos et al., 2014) which are summarized below. Each analysis considers all three phase variation models. The model priors for each analysis are described below and summarized in Table A1 for convenience. Within each analysis pipeline, models are selected based on the Bayesian Information Criterion (BIC). We cannot choose our fiducial models between our three analyses using the BIC as there are significant differences between the number of data used in each analysis because of different $\sigma$-clipping and binning. Instead, we choose to discriminate between the three analyses by selecting the model with the largest log-likelihood per datum, $\ln(L)/N_{\rm data}$; we therefore adopt the preferred models from M. Zhang’s analyses as our fiducial models. The fiducial reductions of the four data sets are presented in Figure 1 and Table 1 (see also the Appendix and Supplementary Information).

3.2.1 Fiducial Reduction and Decorrelation Procedure

For reasons described below, M. Zhang’s analyses were selected as our fiducial analyses and follow their previous work (Zhang et al., 2018). In this analysis, we perform aperture photometry with a radius of 2.7 pixels on the Spitzer BCD files to get the raw flux for all frames. The background is calculated by excluding all pixels within a radius of 12 pixels from the star, rejecting outliers using sigma clipping, and then calculating the biweight location of the remaining pixels. We then bin the background-subtracted raw fluxes with a bin size of 64, discard the first 0.05 days of data, and perform fitting with emcee (Foreman-Mackey et al., 2013). The fitting uses 250 walkers that walk for 20 000 burn-in steps and 20 000 post-burn-in steps. Our instrumental model uses first order PLD for all data except the 2010 3.6 $\mu$m data, in which case we find that second order PLD minimizes BIC. Aside from PLD, the instrumental model also includes a linear slope with respect to time. We fit for the following parameters, all with uniform priors: transit time, eclipse time, $R_{p}/R_{*}$, eclipse depth (assumed to be constant over each ${\sim 33}$ hour phase curve), sinusoidal phase variation amplitudes ($C_{1}$ and $D_{1}$ for the first order sinusoid, and $C_{2}$ and $D_{2}$ if running a second order sinusoidal model), photometric error, slope in flux with time, and PLD coefficients. We fixed $P$, $a/R_{*}$, and $i$ to the highly precise values from the literature (Collins et al., 2017) as they are poorly constrained by our observations. As limb darkening is not that important in the Spitzer bands, we adopt the closest model from a grid of 1D stellar models (Sing, 2010).

Our fiducial analyses find that the photon noise limits are 652 ppm and 637 ppm for the 2010 and 2013 3.6 $\mu$m observations, respectively, and the limits for the 2010 and 2013 4.5 $\mu$m observations are 866 ppm and 860 ppm, respectively. The differences between these two is likely due to the star falling on parts of the detector with slightly different sensitivities, as well as varying aperture sizes. The fitted photometric standard deviation from our fiducial analyses are 950 ppm and 976 ppm for the 2010 and 2013 3.6 $\mu$m observations (1.46 and 1.53 times greater than the photon noise limit). For the 2010 and 2013 4.5 $\mu$m observations, the fitted photometric standard deviations are 1134 ppm and 1158 ppm (1.31 and 1.35 times greater than the photon noise limit). Figures showing the normalized raw, decorrelated, and residual fluxes from all four phase curves analyzed with M. Zhang’s pipeline can be found in the Appendix (Figures A5 and A6).

3.2.2 T. Bell’s Reduction and Decorrelation Procedure

Reduction and decorrelation of these data follow Dang et al. (2018) and are summarized here. We convert the pixel intensity from MJy/str to electron counts and mask bad pixels, i.e., $4\sigma$ outliers with respect to the median of that pixel in the datacube as well as any NaN pixels. We discard all frames with a bad pixel within the aperture used for photometry. We also discard every first frame from each data cube from the 2010 observations and every first and second frame from each data cube for the 2013 observations because these frames consistently show the presence of significant outliers compared to other frames within the same data cube. The effect of this sigma clipping is minimal, given that model fitting is performed on the median binned values from each data cube. There is another star (other than WASP-12A,B,C) that falls on the detector but lies outside the considered photometric apertures (${\sim}10^{\prime\prime}$ away); we place a $3\times 3$ pixel mask around this star to ensure that it does not bias the background subtraction.

We then perform aperture photometry on each individual frame, with an aperture at the fixed pixel-location (15,15), and centroids were found using a flux-weighted mean algorithm and later used for decorrelation. Apertures ranging from 2 to 5 pixels in radius were considered as well as two different aperture edges: hard (the pixel’s flux is included if the centre of the pixel lies within the aperture) and soft (each pixel is weighed by the exact fraction of its area included within the aperture). While some flux will be lost by smaller apertures, a smaller aperture better allows us to remove intra-pixel sensitivity variations, which are the dominant source of noise in our data. We select the aperture radius and edge which resulted in the lowest RMS after a copy of the raw data were smoothed by a boxcar filter of width 5 data cubes (${\sim}11$ minutes which is approximately half the ingress/egress duration) to remove features such as transits, eclipses, and phase variations. Tests run with apertures centred on the flux-weighted mean derived centroids showed that the RMS was $>100$ ppm higher than the fixed position apertures. For the 2010 data, we selected a hard-edged 2.5 pixel radius aperture for the 4.5 $\mu$m data and an soft-edged 4.3 pixel radius aperture for the 3.6 $\mu$m data; the previous analysis of these data (Cowan et al., 2012) used IDL’s approximation on a soft-edged 2.5 pixel radius aperture for both wavelengths. For the 2013 data, we selected a hard-edged 3.2 pixel radius aperture for the 4.5 $\mu$m data and an soft-edged 2.9 pixel radius aperture for the 3.6 $\mu$m data. Before decorrelating and analyzing the data, we first bin the flux and centroid measurements from all 64 frames within a data cube using a median to reduce noise and decrease computation time. On average, each of our models take ${\sim}0.5$ hour to fit to the binned data, and computation time grows linearly with the number of data points, so running each of the different models on unbinned data is not feasible.

T. Bell’s analyses used various systematic models as implemented in the open-source Spitzer Phase Curve Analysis (SPCA; Dang et al., 2018) pipeline111https://github.com/lisadang27/SPCA. In particular, we used two-dimensional polynomials of order 2 through 5 and BiLinear Interpolated Subpixel Sensitivity (BLISS) mapping. The two-dimensional polynomials (Charbonneau et al., 2008) assume the sensitivity of the detector can be described by an $n$th-order 2D polynomial in the measured centroid. BLISS mapping (Stevenson et al., 2012a; Ingalls et al., 2016; Schwartz & Cowan, 2017) is a non-parametric method to account for the intra-pixel sensitivity variations which requires accurate centroid measurements; when fitting BLISS models we adopt an 8$\times$8 grid of knots. For the 2013 observations at 3.6 $\mu$m, we also needed to add a slope in time to remove residual red noise.

Models were fit using the Markov Chain Ensemble Sampler emcee (Foreman-Mackey et al., 2013). The orbital parameters of WASP-12b in the literature (Collins et al., 2017) have smaller errors than we can achieve with our photometry. Additionally, numerous searches for eccentricity have found that WASP-12b’s orbit is best described by a circular orbit (Campo et al., 2011; Croll et al., 2011; Bailey & Goodman, 2019), so we set the orbital eccentricity to zero. Several orbital parameters are poorly constrained by a single phase curve observation compared to the literature values, so we adopted the following Gaussian priors to marginalize over the uncertainties in the literature values: $t_{0}=56176.16825800\pm 0.00007765$ (BMJD), $a/R_{*}=3.039\pm 0.034$, $i=83.37^{\circ}\pm 0.68^{\circ}$ (Collins et al., 2017). The orbital period is known to within 12 ms, so we simply fixed it at $1.09142030$ days (Collins et al., 2017). The parameters that were always fitted were $t_{0}$, $R_{p}/R_{*}$, $a/R_{*}$, $i$, $F_{\rm day}/F_{*}$, two quadratic limb darkening parameters (Kipping, 2013) $q_{1}$ and $q_{2}$, and the first order sinusoidal amplitudes $C_{1}$ and $D_{1}$. In some models, we also fitted $R_{p,2}/R_{*}$ or $C_{2}$ and $D_{2}$. A number of detector parameters were also fitted, with the exact number depending on the detector model used.

For T. Bell’s apertures, the photon noise limits are 578 ppm and 566 ppm for the 2010 and 2013 3.6 $\mu$m observations, respectively, and the limits for the 2010 and 2013 4.5 $\mu$m observations are 795 ppm and 791 ppm, respectively. The fitted photometric standard deviation from T. Bell’s analysis are 1493 ppm and 1246 ppm for the 2010 and 2013 3.6 $\mu$m observations (2.58 and 2.20 times greater than the photon noise limit). For the 2010 and 2013 4.5 $\mu$m observations, the fitted photometric standard deviations are both 1440 ppm (1.81 and 1.82 times greater than the photon noise limit). The fitted data and red noise tests from these analyses can be found in the Supplementary Information (Figures A7 and A8).

3.2.3 P. Cubillos’ Reduction and Decorrelation Procedure

The models run by P. Cubillos use the Photometry for Orbits, Eclipses, and Transits (POET) pipeline (Stevenson et al., 2010, 2012a, 2012b; Campo et al., 2011; Nymeyer et al., 2011; Cubillos et al., 2013, 2014). The POET pipeline starts by flagging bad pixels from the Spitzer BCD files using the permanent bad pixel masks and performing a sigma-rejection routine. Next, it estimates the target center position either fitting a two-dimensional Gaussian function or calculating the least asymmetry (Lust et al., 2014). Then it obtains raw light curves by applying a circular interpolated aperture photometry, testing several aperture radii between 2.0 and 4.0 pixels.

To determine the optimal centroiding method and photometry aperture, POET minimizes the standard deviation of the residuals, and minimizes time-correlated noise at timescales equal and larger than the transit duration (estimated through the time-averaging method). Least asymmetry centroiding outperformed Gaussian centring for all datasets, except the 2013 4.5 $\mu$m observation. The optimal apertures were 2.5 and 3.0 pixels (2010) and 4.0 and 2.0 pixels (2013) for the 3.6 and 4.5 $\mu$m observations, respectively. In any case, all relevant astrophysical parameters vary within their uncertainties as we vary the centroiding and photometry.

POET models the unbinned light curves, simultaneously fitting the astrophysical phase curve and the telescope systematics. The systematics model consists of the non-parametric BLISS intrapixel model, for which we set the map’s bin size equal to the RMS of the frame-to-frame target position (0.01 pixels), and require at least 8 points per bin. For the 2013 4.5 $\mu$m observation, we also apply a linear time-dependent ramp with the slope as a free parameter.

The astrophysical model consists of transit and eclipse models (Mandel & Agol, 2002), combined with the sinusoidal and ellipsoid models described in the Methods section. The transit free fitting parameters are the epoch, ratio between the planetary and stellar radii, cosine of inclination, semi-major axis to stellar radius ratio, stellar flux, and quadratic limb-darkening coefficients. The eclipse free fitting parameters are the midpoint, duration, depth, and ingress duration (setting the egress duration equal to the ingress duration). We adopt uniform priors for all parameters, except for $\cos(i)$ and $a/R_{\rm*}$, which have Gaussian priors, and kept the orbital period fixed (same values as in T. Bell’s reduction and decorrelation procedure; Collins et al., 2017).

POET incorporates the MC3 statistical package (Cubillos et al., 2017) to find the best-fitting parameter values (using Levenberg-Marquardt optimization) and uncertainties (using a differential-evolution Markov Chain Monte Carlo algorithm; ter Braak & Vrugt, 2008), requiring the Gelman-Rubin statistic (Gelman & Rubin, 1992) to be within 1% of unity for each free parameter for convergence. POET uses Bayesian hypothesis testing to select the model best supported by the data, selecting the lowest BIC model. The POET results support the independent results of the other pipelines. Both 4.5 $\mu$m observations strongly favour the second order sinusoidal model, while both 3.6 $\mu$m observations strongly favour the first order sinusoidal model.

For P. Cubillos’ apertures, the photon noise limits are 4886 ppm and 6664 ppm for the 2010 and 2013 3.6 $\mu$m observations, respectively, and the limits for the 2010 and 2013 4.5 $\mu$m observations are 8273 ppm and 7988 ppm, respectively. The fitted photometric standard deviation from P. Cubillos’ analysis are 6915 ppm and 7360 ppm for the 2010 and 2013 3.6 $\mu$m observations (1.41 and 1.10 times greater than the photon noise limit). For the 2010 and 2013 4.5 $\mu$m observations, the fitted photometric standard deviations are 9130 ppm and 8658 ppm (1.10 and 1.08 times greater than the photon noise limit). The fitted data and red noise tests from these analyses can be found in the Supplementary Information (Figures A9 and A10).

4 Results

4.1 Comparison Between Pipelines and Epochs

All three independent analyses confirm the presence of strong and persistent second order sinusoidal variations at 4.5 $\mu$m and the non-detection of these variations at 3.6 $\mu$m. The fitted phase curves parameters for the preferred models from all three independent pipelines are summarized in Figure 1 and Table A2. See the Supplementary Information for tabulated values for all considered models. The astrophysical parameters at 3.6 $\mu$m and 4.5 $\mu$m are mostly consistent between all three analyses, with the preferred models from the three analyses generally differing by $<2\sigma$. In the few cases where one model differs from the others by more than $2\sigma$, the other two models are consistent with each other at a level of $<1\sigma$. Also, there is low-frequency noise in the 2010 4.5 $\mu$m residuals between the first eclipse and the end of the transit that is seen by all three analysis pipelines; the source of these variations is not understood.

From 2010 to 2013, the three pipelines show that all 4.5 $\mu$m phase curve and planetary parameters remain constant within $<2\sigma$. Most of the phase curve and planetary parameters at 3.6 $\mu$m also remain constant between the two observing epochs, with the main exception being the phase offset calculated using the first order sinusoidal terms. M. Zhang’s, P. Cubillos’, and T. Bell’s pipelines find that it changes by $6.4\sigma$ ($46.2^{\circ}$), $7.7\sigma$ ($46.6^{\circ}$), and $3.1\sigma$ ($28.1^{\circ}$), respectively. All three pipelines also agree that the sign of the hotspot offset changes between the two observing epochs, with the offset being “eastward” (before eclipse) in 2010 and “westward” (after eclipse) in 2013. It is interesting to note, however, that over this same time span both the first and second order sinusoidal phase offsets from the 4.5 $\mu$m observations do not change. No other parameter is found by all three analyses to vary by more than $3\sigma$ between the two observing epochs. Finally, if the first order sinusoidal phase variations are entirely attributable to WASP-12b’s temperature map, our 2013 observations at 3.6 $\mu$m exhibit a $13\hbox to0.0pt{.\hss}^{\circ}6\pm 3\hbox to0.0pt{.\hss}^{\circ}8$ westward hotspot offset. This may be a demonstration that eastward hotspot offsets are less ubiquitous than previously believed, with westward hotspot offsets reported for planets with irradiation temperatures spanning 2200–3700 K (Dang et al., 2018; Zhang et al., 2018; Wong et al., 2016).

4.2 Physical Sources

The previously favoured explanation for the double peaked phase curve reported for WASP-12b by Cowan et al. (2012) was detector systematics, but this hypothesis is now strongly disfavoured. To date, 23 papers have been published with new Spitzer phase curves of 18 different exoplanets (Harrington et al., 2006; Knutson et al., 2007; Cowan et al., 2007; Knutson et al., 2009b, a; Laughlin et al., 2009; Crossfield et al., 2010; Cowan et al., 2012; Knutson et al., 2012; Crossfield et al., 2012a; Lewis et al., 2013; Maxted et al., 2013; Zellem et al., 2014; Wong et al., 2015; de Wit et al., 2016; Wong et al., 2016; Krick et al., 2016; Demory et al., 2016; Wong et al., 2016; Stevenson et al., 2017; de Wit et al., 2017; Zhang et al., 2018; Dang et al., 2018; Kreidberg et al., 2018). Of these numerous observations, WASP-12 is the only system which has shown strong a double peaked phase curve not once, but twice. The observing strategy also differed between these two sets of WASP-12b phase curves, with the number and timing of AORs changing and the addition of PCRS Peak-Up before the 2013 observations. The consistency between the two sets of phase curves suggests that the observations probe an astrophysical source which does not vary significantly over a $\sim$3 year timescale. Cowan et al. (2012) suggested that tidal distortion and/or mass loss might be able to explain the Spitzer observations, if this signal was indeed astrophysical in nature. We explore these and other potential sources of emission below.

4.2.1 Tidal Distortion

One potential cause of second order sinusoidal variations is tidal deformation of the host star, as is seen at optical wavelengths for HAT-P-7 (Welsh et al., 2010) and WASP-18 (Shporer et al., 2019). However, stellar distortion is expected to be negligible for WASP-12 (Leconte et al., 2011a). We verified this by numerically solving for the equipotential stellar/planetary surfaces using the dimensionless Roche potential (see Appendix C). However, since the star contributes significantly more flux than the planet, we ran simple simulations of both the star and planet including the effects of gravity darkening to assess their expected amplitudes of ellipsoidal variations. We find that the stellar ellipsoidal variations are approximately the same amplitude at 3.6 $\mu$m and 4.5 $\mu$m and the amplitude of the stellar variations are far smaller than the observed amplitudes; we therefore conclude that tidal bulges on the host star cannot be the source of the strong second order variations observed at 4.5 $\mu$m. Our predicted ellipsoidal variations for WASP-12b are consistent with our limits on second order sinusoidal variations at 3.6 $\mu$m, but significantly under-predict the observed amplitude at 4.5 $\mu$m (also see the implied dimensions of the best-fit ellipsoidal variation model shown in the right panel of Figure 2). If we interpret the second order sinusoidal variations at 4.5 $\mu$m as planetary ellipsoidal variations, this would require the 4.5 $\mu$m photosphere to be significantly higher up than 3.6 $\mu$m as the layers nearer the Roche lobe are more distorted. However, this increased radius at 4.5 $\mu$m is inconsistent with the smaller transit depth at 4.5 $\mu$m compared to 3.6 $\mu$m. We therefore conclude that tidal distortion of the planet is also not the source of the strong second order variations observed at 4.5 $\mu$m.

4.2.2 Stellar Variability and Inhomogeneities

Stellar variability is also unlikely to be the cause of these observations given the comparable phase of the second order variations in the two data sets. For reference, the WASP-12BC dilution correction term is ${\sim}400$ ppm while the observed amplitude of the second order sinusoidal variations is ${\sim}2000$ ppm. Additionally, variability in WASP-12A is only predicted to modulate the planetary signal at a level of $\sim$1 ppm, and variability in WASP-12B,C should also only contribute at the level of ${\sim}1$ ppm (while these M-dwarfs should be ${\sim}10\times$ more variable, they contribute ${\sim}10\times$ less flux; Zellem et al., 2017). We therefore rule out standard stellar variability as the source of the strong second order sinusoidal variations seen at 4.5 $\mu$m. If the second order variations were produced by unusually strong inhomogeneities on the host star, both the sub-planet longitude and the anti-planet longitude would need to be darker than intermediate longitudes — this would imply star–planet interactions. However, these inhomogeneities would also need to be much more pronounced at 4.5 $\mu$m which would not be expected for the ${\sim}6000$ K star.

4.2.3 Mass Loss

There is significant observational evidence from near ultra-violet (NUV) transit observations that WASP-12b is undergoing mass loss and that there is a bow shock in the system (Fossati et al., 2010; Haswell et al., 2012; Nichols et al., 2015) (see the Supplementary Information). A potential explanation for the unusual 4.5 $\mu$m phase curve is that there is gas being stripped from the planet which emits more strongly within the 4.5 $\mu$m bandpass than the 3.6 $\mu$m bandpass. The observations favour a stream of dense gas stripped from the planet flowing directly toward/away from the host star or some other elongated patch of hot gas whose long axis is parallel to the star–planet axis, such as an accretion hot spot. Double-peaked phase curves have been seen for dwarf novae CVs, such as WZ Sge (Skidmore et al., 1997), however this feature was seen through the ultra-violet to infrared; for CVs these variations have been attributed to tidal distortion or an optically thick hot spot in an otherwise optically thin accretion disk (e.g., Skidmore et al., 1997)

The source of the 4.5 $\mu$m variations in the WASP-12 system must lie near the star–planet axis since there is no significant detection of an occultation of the source when the planet is not in transit or eclipse. Additionally, our Spitzer observations demonstrate that the planetary radius appears ${\sim}8$% ($11\sigma$) smaller at 4.5 $\mu$m than at 3.6 $\mu$m which is in disagreement with model predictions (Burrows et al., 2007; Burrows et al., 2008; Cowan et al., 2012); this rules out the transit of a large exosphere that is opaque at 4.5 $\mu$m as this would make the planetary radii at the two wavelengths even more discrepant.

As shown in Table 1, the fitted second order offsets at 4.5 $\mu$m are consistent with being oriented along the star–planet axis ($90^{\circ}$). However, previously published 3D magnetohydrodynamic (MHD) numerical simulations of hypothetical exoplanet systems mostly produced gas flows that significantly lead the star–planet axis (Matsakos et al., 2015). Indeed, gas streaming from the planet’s L1 Lagrange point on a ballistic trajectory should flow $\theta_{\rm gas}=53$∘$.\,4$ ahead of the star–planet axis (Lai et al., 2010) as angular momentum is conserved (see Figure 2 for a schematic depiction). Assuming our observations probe the gas stream, this prediction is $27\sigma$ discrepant with our offset of $4$∘. $0\pm 2$∘$.\,1$ behind the star–planet axis found by averaging the offsets from the two fitted second order sinusoids at 4.5 $\mu$m. This discrepancy could potentially be explained if the 4.5 $\mu$m emitting area is much closer to the planet and is still aligned along the star–planet axis and then becomes more diffuse and flows ahead of the planet as it continues to fall toward the host star.

Alternatively, stellar effects could channel the infalling stream directly toward the star, but this may be inconsistent with past NUV transit observations (Fossati et al., 2010; Haswell et al., 2012; Nichols et al., 2015). One previously published 3D MHD model (Matsakos et al., 2015) did exhibit a stream of gas directly along the star–planet axis (their name for this model was ‘FvrB’). This model has high stellar ultra-violet (UV) flux, a low escape speed from the planet, the planet near to its host star, and a strong planetary magnetic field. In this model, the planet is experiencing Roche lobe overflow with a planetary wind that is weak compared to the stellar wind, producing an approximately linear stream of gas along the star–planet axis as well as a lower density tail trailing behind the planet (Matsakos et al., 2015). The non-detection of the gas trailing behind the planet could be explained if the gas has a lower density and/or has a lower temperature. As the dense gas stream in the ‘FvrB’ model is aligned along the star–planet axis, it may not contribute significantly to the transit depth and may remain consistent with the smaller apparent radius at 4.5 $\mu$m compared to 3.6 $\mu$m. Radiative transfer simulations based on the ‘FvrB’ mass-loss model (Matsakos et al., 2015) would allow for this hypothesis to be tested.

4.3 Radiation Mechanisms

4.3.1 Blackbody Emission

The discrepant second order sinusoidal amplitudes at 3.6 $\mu$m and 4.5 $\mu$m can be explained by one of two emission mechanisms. First, blackbody emission could allow for greater flux at 4.5 $\mu$m than at 3.6 $\mu$m if the gas is sufficiently cool that the 3.6 $\mu$m bandpass lies on the Wien side of the blackbody curve; this scenario would allow us to place an upper limit on the temperature and spatial extent of the emitting gas, which we pursue below.

Using the host star’s effective temperature of $6300\pm 150$ K (Hebb et al., 2009), we assume the host star emits as a blackbody and convert the second order sinusoidal curves from units of $F_{\rm day}/F_{*}$ to $B_{\lambda}$ as shown in the middle panel of Figure 3. We adopt the fiducial 4.5 $\mu$m parameters from 2010, but set the phase offset to 90*∘* since there is no evidence for an appreciable offset from the star–planet axis. We then take the best-fit and the $1\sigma$ and $2\sigma$ upper limits on the amplitude of the 3.6 $\mu$m second order sinusoidal variations from M. Zhang’s analysis using the second order astrophysical model. We assume that none of the flux seen during planetary transit/eclipse is from emission by the gas. By assuming the emitting area is the same at 3.6 $\mu$m and 4.5 $\mu$m, we can use the relative amounts of flux at these two wavelengths to determine the blackbody temperature of the gas.

Given the assumption that our observations are explained by blackbody emission, we can then place a $2\sigma$ upper limit of $619$ K on the gas temperature. For reference, a temperature of $816$ K would provide equal flux in both bandpasses. Attributing any of the “nightside” flux to emission from the gas only lowers this limit further. Also, as WASP-12b’s skin temperature (Goody & Walker, 1972) is $0.5^{0.25}\,T_{\rm b,day}\approx 2500$ K, this gas cannot be the upper layers of the planet’s atmosphere.

By taking the ratio between the flux emitted by the gas and that emitted by the star, we can determine the effective emitting area required to produce the phase curve observations. As shown in the bottom panel of Figure 3, less emission at 3.6 $\mu$m requires lower temperature gas and therefore a larger emitting area. We can therefore place a $2\sigma$ lower-limit on the effective emitting area of the gas of $0.98$ times the planet’s transiting area when seen at planetary quadrature, given the assumption that our observations are explained by blackbody emission. Attributing any of the nightside flux to emission from the gas slightly increases this limit and allows for a non-zero emitting area during planetary transit and eclipse.

4.3.2 CO Emission

An alternative explanation for the increased flux at 4.5 $\mu$m is emission by CO which has its strong $\Delta V=1$ band around 4.5 $\mu$m (see Figure A1 in the Appendix for the CO line intensities); CO emission has previously been predicted for gas lost from WASP-12b (Li et al., 2010; Deming et al., 2011). The CO molecule should be dissociated in the planetary upper atmosphere due to the strong UV and X-ray flux from the host star which also drives most of the observed atmospheric escape seen at NUV wavelengths (Fossati et al., 2010; Haswell et al., 2012; Nichols et al., 2015); the dissociation energy of CO corresponds to a wavelength of roughly 110 nm. However, the atomic carbon and oxygen from the upper layers of the planet’s atmosphere could recombine in a gas stream where the density gets higher because of stellar wind confinement and the “shadow effect” from the material in the stream closer to the star. Given a gas temperature profile (Salz et al., 2016) and our calculations of the thermal dissociation fraction of CO using the Saha equation (Bell & Cowan, 2018), we find that any CO emission must either be produced within ${\sim}0.1~{}R_{p}$ of the planet’s surface or beyond $2.5~{}R_{p}$. In the case of a bow shock supported by mass loss from the planet, gas temperatures are predicted to reach $10^{3}$–$10^{4}$ K (Turner et al., 2016) which should allow for stable CO, provided there is sufficient UV shielding from gas nearer to the star. Simulations of the behaviour of CO in these environments are required to determine the feasibility of this molecule recombining once in a stream and emitting sufficiently strongly to explain our observations.

4.4 A Note on Eclipse Depths

It is important to note that our reported “eclipse depths” ($F_{\rm day}/F_{*}$) are measured with respect to the phase curve value expected at the centre of eclipse and are not measured with respect to pre-ingress and post-egress flux measurements as would be the case for observations of only the eclipse. Given our fitted phase curve parameters for WASP-12b, the difference between our reported value and using the average of pre-ingress and post-egress baselines is ${\sim}9$% of $F_{\rm day}/F_{*}$ at both Spitzer bandpasses (assuming these baseline durations are both the same duration as the eclipse duration). This bias in eclipse observations occurs because the phase variations before ingress and after egress are flattened out by most decorrelation routines when solely observing the eclipse. For most exoplanets whose phase variations will be concave down around eclipse (like WASP-12b when seen at 3.6 $\mu$m), eclipse observations will underestimate $F_{\rm day}/F_{*}$. For the unusual case of WASP-12b’s 4.5 $\mu$m phase variations which are concave up near eclipse, eclipse observations will overestimate $F_{\rm day}/F_{*}$. This effect is particularly important for short period planets which undergo significant rotation throughout the duration of eclipse observations and whose strong day–night temperature contrast cause strong phase variations over this time span. Among other things, this may explain the discrepancies between reported 3.6 $\mu$m and 4.5 $\mu$m eclipse depths from full-orbit phase curves (Cowan et al., 2012) and eclipse-only observations (Madhusudhan et al., 2011; Stevenson et al., 2014b), and the associated inference of C/O ratio. See Figure 4 for a demonstration of this effect.

5 Discussion and Conclusions

By independently analyzing and then combining two sets of 3.6 $\mu$m and 4.5 $\mu$m Spitzer phase curves of the UHJ WASP-12b, we have conclusively detected strong and persistent second order sinusoidal variations at 4.5 $\mu$m and placed stringent upper limits on these variations at 3.6 $\mu$m. These observations of WASP-12b raise several questions which will require further study to resolve.

Our two emission hypotheses could be distinguished with phase curve observations of the $\sim$1.6 $\mu$m and/or 2.29 $\mu$m CO emission bands and/or with phase curve observations at wavelengths longer than 5 $\mu$m which should exhibit strong second order sinusoidal variations if the 3.6 $\mu$m vs. 4.5 $\mu$m amplitude discrepancy is the result of blackbody emission. The high precision and wavelength coverage achievable with the James Webb Space Telescope should allow these two emission hypotheses to be tested. The $\sim$1.6 $\mu$m CO emission band also lies within the Hubble/WFC3 bandpass and may be detectable with phase curve observations.

Critically, future models must also address the fact that the fitted planetary radius is significantly smaller at 4.5 $\mu$m than at 3.6 $\mu$m; this may be the result of unocculted emitting gas. Combined hydrodynamic and radiative transfer simulations are required to fully understand this system. These simulations will allow us to determine the location and spatial extent of the emitting gas, and they may resolve the apparent tension between the constraint from these observations that the gas is well aligned with the star–planet axis, while NUV observations which probe lower density gas show that the gas flows significantly ahead of the planet. Understanding the nature of the increased emission at 4.5 $\mu$m will also require modelling the mass loss and the UV dissociation and potential recombination of CO molecules as they flow from the planet’s upper atmosphere through a gas stream and potentially experience a shock. These models may also assist in understanding the observed hot spot variability seen at 3.6 $\mu$m.

Finally, while WASP-12b is one of the exoplanets closest to overflowing it’s Roche lobe (see Figure A2), there are several other UHJs with similar characteristics with published Spitzer phase curves that do not show strong second order sinusoidal variations at 4.5 $\mu$m: particularly WASP-19b (Wong et al., 2016), WASP-33b (Zhang et al., 2018), and WASP-103b (Kreidberg et al., 2018). One potential explanation is that WASP-12b’s orbit may be decaying (Maciejewski et al., 2016; Patra et al., 2017) while the other exoplanets may be more stable; this could potentially be explained if WASP-12b was locked in a high obliquity state due to a resonance with a perturbing planet which could drive orbital decay and inflate the planet beyond it’s Roche lobe (Millholland & Laughlin, 2018). Alternatively, the high energy irradiation from WASP-12A may be stronger than the other UHJ host stars. Further research is required to understand why WASP-12b is the only exoplanet known to be exhibiting these exceptionally strong second order sinusoidal variations at 4.5 $\mu$m.

Acknowledgements

T.J.B. acknowledges support from the McGill Space Institute Graduate Fellowship, the Natural Sciences and Engineering Research Council of Canada’s Postgraduate Scholarships-Doctoral Fellowship, and from the Fonds de recherche du Québec – Nature et technologies through the Centre de recherche en astrophysique du Québec. The research leading to these results has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement no. 679633; Exo-Atmos). We have also made use of open-source software provided by the Python, Astropy, SciPy, and Matplotlib communities.

Appendix A: Correction for Dilution by Stellar Companions

To correct for the dilution of our lightcurves by the nearby stellar companions WASP-12BC, we apply the dilution factors from Stevenson et al. (2014a): $\alpha_{\rm comp}=0.1149\pm 0.0039$ and $0.1196\pm 0.0042$ for 3.6 $\mu$m and 4.5 $\mu$m respectively. Since our phase curve amplitudes are normalized by the eclipse depth, no corrections need to be made to $C_{1}$, $D_{1}$, $C_{2}$, or $D_{2}$. Additionally, while the planetary radii need to be corrected for dilution from WASP-12BC, the ratio $R_{p,2}/R_{p}$ remains the same for models with ellipsoidal variations. Following Stevenson et al. (2014b); Stevenson et al. (2014a), the multiplicative correction factor is

[TABLE]

where $g(\beta,\lambda)$ is the fraction of WASP-12BC’s flux which falls within our aperture of size $\beta$. We estimated $g(\beta,\lambda)$ using STINYTIM222http://irsa.ipac.caltech.edu/data/SPITZER/docs/dataanalysistools/tools/contributed/general/stinytim/, the point response function modelling software for Spitzer. We made $10\times$ oversampled point response functions calculated at the pixel position (25,25) assuming a $T=3660$ blackbody source (the effective temperature of WASP-12BC; Stevenson et al., 2014b). We found $g(2.5,4.5~{}\mu$m) = 0.8147, $g(3.2,4.5~{}\mu$m) = 0.8608, $g(4.3,3.6~{}\mu$m) = 0.9089, and $g(2.9,3.6~{}\mu$m) = 0.8580. For the $3\times 3$ pixel stamp used in M. Zhang’s PLD analyses, we find $g(3{\times}3,3.6~{}\mu$m) = 0.6518 and $g(3{\times}3,4.5~{}\mu$m) = 0.6291. For P. Cubillos’ analyses, we find $g(3.0,4.5~{}\mu$m) = 0.8533, $g(2.5,4.5~{}\mu$m) = 0.6957, $g(2.5,3.6~{}\mu$m) = 0.8254, and $g(4.0,3.6~{}\mu$m) = 0.9015. We also checked $g(2.25,3.6~{}\mu$m$)$ and $g(2.25,4.5~{}\mu$m$)$ to compare our calculation to that of Stevenson et al. (2014b); we find values of 0.8007 and 0.7586, where Stevenson et al. (2014b) found 0.7116 and 0.6931. This discrepancy is likely caused by an incorrect angular separation used in the previous work’s calculation.

The planet’s radius was then corrected using

[TABLE]

with the elongated axis, $R_{p,2}$, in bi-axial ellipsoid models corrected similarly. The dayside flux was corrected using

[TABLE]

with the white noise amplitude, $\sigma_{F}$, corrected similarly.

Appendix B: Computing Astrophysical Parameters

Tables A2–A9 present many of the fitted astrophysical values from all models run in all three independent analyses. $T_{b,\rm day}$ and $T_{b,\rm night}$ are the apparent brightness temperatures of the planet’s day and night hemispheres which we calculate using only the contribution from the first order sinusoid. In doing so, we are assuming that the second order sinusoidal variations are attributable to something other than the planet, although the second-order sinusoidal variations end up having negligible contributions during transit and eclipse anyway. These brightness temperatures are calculated by inverting the Planck function (Cowan & Agol, 2011), using

[TABLE]

where $h$ is Planck’s constant, $c$ is the speed of light, $k_{B}$ is the Boltzmann constant, $\lambda$ is the wavelength. For $T_{b,\rm day}$, $\psi=(F_{\rm day}/F_{*})(R_{p}/R_{*})^{-2}$, and for $T_{b,\rm night}$, $\psi=(F_{\rm day}/F_{*})(1-2\,C_{1})(R_{p}/R_{*})^{-2}$. The stellar brightness temperature, $T_{*,\rm b}$ was calculated by fitting blackbodies to the relevant wavelengths from a PHOENIX stellar model (Husser et al., 2013) with previously measured (Hebb et al., 2009) values of $T_{\rm*,eff}=6300$ K and $\log(g)=4.5$. We find $T_{*,\rm b}=6000$ K for 4.5 $\mu$m and 5800 K for 3.6 $\mu$m. The tabulated first and second order offsets are measured in degrees after the secondary eclipse and are calculated using:

[TABLE]

[TABLE]

Appendix C: Tidal Distortion Calculations

To assess the impact of stellar and planetary tidal distortion, we model the stellar/planetary surfaces using the dimensionless Roche potential, defined by

[TABLE]

where $r$ is the distance from the host star, $\theta$ is the polar angle, $\phi$ is the azimuthal angle, and $q$ is the mass ratio, $M_{*}/M_{p}$. We find that the star’s radius should be 0.0085% longer along the star–planet axis compared to the perpendicular equatorial axis, while the planet’s radius should be 5.5% longer along the star–planet axis compared to the dawn–dusk axis seen at transit.

We first assume that the planet and star have a constant temperature of 3000 K and 6300 K, respectively, and then perturb these temperatures to account for gravity darkening using the $T_{\rm eff}\propto g_{\rm eff}^{\beta}$ model (Espinosa Lara & Rieutord, 2011) where $\beta$ is 0.24 for the appreciably distorted planet and 0.25 for the more spherical host star. Next, we convert these temperature maps into flux maps using the Planck blackbody function. We then compute disk-integrated phase curves (Cowan et al., 2013) while also accounting for the variations in apparent areas of the two objects. Our calculations show that the planet’s expected variations are only ${\sim}3.5$ times stronger than that of the host star at Spitzer/IRAC wavelengths (see Figure A3 for a depiction).

Our predicted ellipsoidal and gravity darkening variations are consistent with past predictions (Budaj, 2011) and with the amplitude of the Zhang PLD model with second order sinusoidal variations fitted to the 3.6 $\mu$m data collected in 2010 (we set the offset to zero as there is no significant detection of an offset in this phase curve). However, the expected ellipsoidal and gravity darkening variations are highly discrepant with the observed amplitude at 4.5 $\mu$m (see Figure A4). Running simulations where the planet fills its Roche lobe ($R_{p,2}/R_{p}\approx 1.4$), our ellipsoidal variations and gravity darkening model would be able to explain the full amplitude of the 4.5 $\mu$m phase curve, but the variations remain mostly monochromatic and the model drastically over predicts the variations in the 3.6 $\mu$m phase curve.

Appendix D: Red Noise Tests

The bottom rows of Figures A5–A10 show the observed standard deviation in the residuals versus the number of data cubes binned together for each lightcurve made using the binrms routine from the Multi-Core Markov-Chain Monte Carlo (MC3)333http://pcubillos.github.io/MCcubed/ package (Cubillos et al., 2017); this allows us to test for any red noise remaining in our residuals (Winn et al., 2007; Cowan et al., 2012). These figures show that minimal red noise remains after our fiducial models have been subtracted from the data (the photometric uncertainty decays roughly as $\sqrt{N_{\rm binned}}$). There is, however, some lower-frequency noise in the 2010 4.5 $\mu$m observations between the first eclipse and the transit that cannot be modelled by any of the three decorrelation pipelines.

Appendix E: NUV Evidence for Mass Loss

Across the NUV, WASP-12b appears to be larger than the planet’s Roche radius, implying significant mass loss (Fossati et al., 2010; Haswell et al., 2012; Nichols et al., 2015). The first Hubble Space Telescope, Cosmic Origins Spectrograph transit observation of WASP-12b (Fossati et al., 2010) also detected an early ingress in the NUV; this suggests the presence of a stream of gas stripped from the planet flowing in toward the star (Lai et al., 2010; Bisikalo et al., 2013b; Matsakos et al., 2015) which forms a bow shock ahead of the planet (Vidotto et al., 2010; Llama et al., 2011; Bisikalo et al., 2013b; Cherenkov et al., 2014; Matsakos et al., 2015; Turner et al., 2016), although the position of this shock can vary (Vidotto et al., 2011; Llama et al., 2013). There is also evidence for variable NUV ingress times (Haswell et al., 2012; Nichols et al., 2015) which suggests variable mass-loss rates and/or a variations in the planet–shock distance (Vidotto et al., 2011). The non-detection of stellar activity indicators from WASP-12A (Knutson et al., 2010; Fossati et al., 2013) may also suggest that WASP-12b is undergoing mass loss. The final resting place of the gas stripped from WASP-12b is debated, with some suggesting an accretion disk interior to the planet’s orbit (Lai et al., 2010; Li et al., 2010) and others suggesting an extended circumstellar torus of gas with the planet embedded inside (Debrecht et al., 2018).

Appendix F: Discussion of Variability

To date, no Spitzer phase curve observation has shown variability in the phase curve offset of an exoplanet, although significant near-infrared variability has been seen for brown dwarfs and isolated planetary mass objects (Artigau et al., 2009; Radigan et al., 2012), and Kepler phase curves of the hot Jupiter HAT-P-7b have been reported to vary (Armstrong et al., 2016). Variability is expected for WASP-12b due to coupling between the planet’s partially ionized atmosphere and the planet’s magnetic field (Rogers, 2017). The timescale of this variability is set by the Alfvén timescale (${\sim}115$ days for WASP-12b assuming magnetic effects occur on the dayside where the atmosphere is dominated by atomic hydrogen; Rogers, 2017; Dang et al., 2018). Variability may also arise in the presence of time-variable cloud coverage, although optically reflective clouds on the planet’s dayside were stringently rejected using Hubble/STIS optical eclipse spectroscopy of WASP-12b (Bell et al., 2017). Any time-variability in the gas streaming from the planet could also obscure different portions of the planet over time and lead to an apparent variation in the 3.6 $\mu$m phase curve.

Appendix G: Model Selection

The preferred model for each phase curve was chosen to be the model with the lowest Bayesian Information Criterion (BIC), defined as

[TABLE]

where $N_{\rm par}$ is the number of model parameters and $N_{\rm dat}$ is the number of data. The log-likelihood is

[TABLE]

where $\sigma_{F}$ is the fitted photometric uncertainty (assumed to be constant throughout the observation) and

[TABLE]

is a measure of the badness-of-fit, where $F_{{\rm obs},i}$ are the observed flux measurements. We adopt the threshold that models with a $\Delta{\rm BIC}\leq 5$ with respect to the favoured model cannot be strongly ruled out.