The Orbit of WASP-12b Is Decaying

We observed four occultations of WASP-12b with the Spitzer Space Telescope in 2019 January and February. The first and last event were separated by 16 planetary orbits. All of the data were obtained with the 4.5 μm channel, in 32 × 32 pixel subarray mode with 2 s exposures. For each event, approximately 11,000 exposures were obtained over a timespan of 7 hr bracketing the 3 hr duration of the occultation.

To reduce the data, we first determined the background level in each exposure by calculating the median flux in the image after excluding the pixels associated with the host star. We subtracted this background level from each image. Then, to measure the pixel location of the centroid of the stellar image, we fitted a two-dimensional Gaussian function to the central 25 pixels of each exposure. Using these centroid positions, we performed circular-aperture photometry, with trial aperture radii ranging from 1.6 to 3.5 pixels in 0.1 pixel increments. We identified a few outliers based on an unusually large deviation in the centroid time series; specifically, we flagged any exposures for which the centroid coordinates were more than 5σ away from the median of the surrounding 10 exposures. The flux values for the offending exposures were replaced by the median flux value within that same 10-exposure window. We also removed from consideration a few data points from the orbit 2026 data set that had obvious image artifacts.

To correct for the well-known effects of intra-pixel sensitivity variations, we used the Pixel Level Decorrelation technique of Deming et al. (2015). We selected a grid of pixels surrounding the stellar image and divided each pixel value by the total flux in that exposure. The intention of this normalization procedure is to eliminate the information from the astrophysical signal (which affects all the pixels), leaving behind only the changes due to pointing fluctuations and differences in pixel sensitivity. Following Equation (4) of Deming et al. (2015), we modeled the light curve as a linear combination of the normalized pixel values ${\hat{P}}_{i}(t)$, a time-dependent trend ft + gt² that accounts for any phase curve variation or long-term instrumental artifacts, a constant offset h, and a geometric eclipse model E(t) with depth D:

$$\begin{eqnarray}&&{F}_{\mathrm{calc}}(t)=\displaystyle \sum _{i=1}^{N}{c}_{i}{\hat{P}}_{i}(t)+{ft}+{{gt}}^{2}+h+{DE}(t).\end{eqnarray} \tag{ 1 }$$

This model could be extended to include cross-terms between the ${\hat{P}}_{i}$ terms and the eclipse model, or higher-order terms in the pixel fluxes (Luger et al. 2016), but we chose not to do so, given the small values of $\sum {c}_{i}{\hat{P}}_{i}$ (<0.01) and the eclipse depth (∼0.005). For a given set of eclipse parameters, we used the batman code (Kreidberg 2015) to calculate E(t). We used linear regression to solve Equation (1) for the coefficients c_i, f, g, and D that provide the best fit to the data F_obs(t).

To speed up computations, it is helpful to reduce the data volume by binning the data in time. Deming et al. (2015) found that the optimized values of the coefficients c_i sometimes depend on the size of the time bins, which they attributed to time-correlated noise. They recommended choosing a bin size that minimizes the amplitude of correlated noise on the timescale of the eclipse. We determined this optimal bin size as follows. We obtained an initial estimate for the occultation time by fitting the unbinned data with a model in which all of the eclipse parameters (apart from the occultation time) were fixed to the values found by Collins et al. (2017). Then, using a fixed eclipse model with this occultation time, we determined the coefficients c_i, f, g, h, D for time-binned light curves, with bin sizes ranging from 2 to 60 exposures (4–120 s). In each case, we examined the residuals by binning them and computing the standard deviation. For uncorrelated noise, the residuals should scale approximately as N^−1/2 where N is the number of exposures per bin. We identified the optimal case as the one for which the residuals best match this expectation. As a further degree of optimization, we repeated this procedure for each choice of aperture radius, and selected the radius that led to the smallest standard deviation of the residuals. Table 2 gives the optimal set of photometric parameters for each observation, while the light curves are provided in Table 1.

Table 2.  New Occultation Midpoints and Depths

Source	Spitzer				WIRC	H19a	V19b
UT Date	2019 Jan 16	2019 Jan 20	2019 Jan 24	2019 Feb 2	2017 Mar 18	2017 Jan 16	2019 Jan 2
Orbit Number	2010	2014	2018	2026	1397	1341	1997
Midpointc	58499.75572	58504.11988	58508.48459	58517.21641	57830.7139	57769.5957	58485.5642
Timing Uncertainty	0.00077	0.00087	0.00091	0.00074	0.0011	0.0014	0.0014
Eclipse Depth (ppm)	${4720}_{-279}^{+289}$	${4243}_{-265}^{+270}$	${3601}_{-262}^{+261}$	${4632}_{-258}^{+266}$	${3232}_{-112}^{+110}$	${1089}_{-71}^{+72}$	${1095}_{-176}^{+175}$
Photometric band	Spitzer 4.5 μm				Ks	i'	V
Aperture radius (pixels)	2.3	2.4	2.1	2.0	⋯	⋯	⋯
Bin size (exposures)	22	22	14	14	⋯	⋯	⋯

Notes.

^aEclipse observed by Hooton et al. (2019). ^bEclipse observed by von Essen et al. (2019). ^cTimes given in ${\mathrm{BJD}}_{\mathrm{TDB}}\mbox{--}{\rm{2,400,000}}$.

^dAperture radius and bin sizes reported in this table only for the Spitzer observations.

Download table as:  ASCIITypeset image

After adopting the optimal aperture and bin size for each of the four observations, we jointly fitted all of the Spitzer data using a single eclipse model. This time, all of the eclipse parameters were allowed to vary, subject to prior constraints. We placed Gaussian priors on the orbital inclination I, planet-to-star radius ratio ${R}_{{\rm{P}}}$/${R}_{\star }$, and orbit-to-star radius ratio a/${R}_{\star }$, based on the best-fit values and uncertainties reported by Collins et al. (2017), shown in Table 3. These parameters are sufficient to describe the loss of light as a function of the planet's position on the stellar disk. To specify the loss of light as a function of time for each event, the timescale ${R}_{\star }P/a$ must also be specified (see Equation (19) of Winn 2010). We did so by holding P fixed at the value 1.09142 days, but importantly, we did not require the interval between occultations to be equal to 1.09142 days. We allowed the occultation midpoints and depths to be freely varying parameters.

Table 3.  WASP-12 System Parameters based on Spitzer Occultations

Parameter	Priora	Best-fit
Period (days)	1.09142	fixed
${R}_{{\rm{P}}}$/${R}_{\star }$	${ \mathcal G }$(0.11785, 0.00054)	0.11786 ± 0.00027
$a/{R}_{\star }$	${ \mathcal G }$(3.039, 0.034)	3.036 ± 0.014
I (deg)	${ \mathcal G }$(83.37, 0.7)	83.38 ± 0.3
e	0.0	fixed
ω (deg)	0.0	fixed

Note.

^aBased on values from Collins et al. (2017). We used Gaussian priors denoted by ${ \mathcal G }(\mu ,\sigma )$.

Download table as:  ASCIITypeset image

Table 2 gives the final fit eclipse times and depths, while Table 3 gives the remaining fit results. The timing precision ranged from 1.0 to 1.3 minutes, which is similar to the results that were achieved by Deming et al. (2015) and Patra et al. (2017) for the same star. Figure 2 shows all four detrended Spitzer light curves, along with the best-fitting eclipse model curves.

Zoom In Zoom Out Reset image size

An occultation of WASP-12b was also observed with the Wide-Field Infrared Camera (WIRC, Wilson et al. 2003) on the Hale 200 inch telescope at Palomar Observatory on 2017 March 18. This observation was made in the K_s band using a new Hawaii-II detector installed on WIRC in 2017 January (Tinyanont et al. 2019) and a near-infrared Engineered Diffuser (Stefansson et al. 2017). We obtained 1828 images with an exposure time of 2 s, spanning 5 hr.

Each image was corrected for dark current, flat field, and bad pixels with the WIRC Data Reduction Pipeline (Tinyanont et al. 2019). We performed circular-aperture photometry of WASP-12 and five comparison stars following the procedure described by Vissapragada et al. (2019). A global background was subtracted from each image using iterative 3σ clipping of the flux values, while a local background was determined for each star using annuli with inner and outer radii of 20 and 50 pixels. Optimizing the pipeline over various aperture sizes found a best circular-aperture radius of 9 pixels. The resulting light curve is provided in Table 1. We modeled the flux time series for WASP-12 as the product of an eclipse model E(t) and a model of systematic effects, consisting of a linear function of time and a linear combination of the fluxes from the five comparison stars:

$$\begin{eqnarray}&&M(t)=\left({ft}+\displaystyle \sum _{i=1}^{5}{c}_{i}{S}_{i}(t\right)\times E(t).\end{eqnarray} \tag{ 2 }$$

We used linear regression to solve for the coefficients in the systematics model, given a choice of parameters for the eclipse model (Vissapragada et al. 2019).

Figure 3 shows the resulting light curve after removing the best-fitting model for the systematic effects. The latter part of the observation was affected by intermittent cirrus clouds, causing sudden and large-amplitude fluctuations in the measured fluxes and leading to larger scatter in the detrended light curve. For this reason, we chose not to fit the data that were obtained during that time period (the gray region in Figure 3). This meant that the egress time and the total transit duration could not be determined from the WIRC data alone. Instead, we held fixed the geometric eclipse parameters at the values taken from the best-fitting model of the Spitzer data (Table 3). We allowed the eclipse depth and midpoint to vary freely. The results of this fit are given in Table 2, and plotted as a red curve in Figure 3. When the entire light curve is fitted, instead of masking out the latter part of the transit, the derived mideclipse time shifts by 0.5σ and has a formal uncertainty that is a factor of 2 smaller.

Zoom In Zoom Out Reset image size

Two other groups have recently reported on observations of occultations of WASP-12b. Hooton et al. (2019) detected the occultation at the 7σ level using the Isaac Newton Telescope on La Palma, in 2017 January. Separately, von Essen et al. (2019) observed three occultations of WASP-12b in 2019 January with the 2.5 m Nordic Optical Telescope (NOT). While neither of these authors published the mideclipse times of their observations, they kindly provided us the light curves. For the observations by von Essen et al. (2019), only the first occultation was securely detected. We only reanalyzed the data from this event. We followed the same detrending procedures that are described in their papers to refit the light curves, holding the eclipse parameters fixed at the values from the best-fitting model to the Spitzer data (apart from the eclipse depth and midpoint). We were able to reproduce their results for the eclipse depths. Table 2 gives the corresponding mideclipse times.

Following Patra et al. (2017), we fitted three models to the timing data. The first model assumes the orbital period to be constant:

$$\begin{eqnarray}\begin{array}{rcl}{t}_{\mathrm{tra}}(N) & = & {t}_{0}+{NP}\\ {t}_{\mathrm{occ}}(N) & = & {t}_{0}+\displaystyle \frac{P}{2}+{NP},\end{array}\end{eqnarray} \tag{ 3 }$$

where N is the number of orbits from a fixed reference orbit,¹⁰ while t₀ is the midtransit time of this reference orbit.

The second model assumes the orbital period to be changing uniformly with time:

$$\begin{eqnarray}\begin{array}{rcl}{t}_{\mathrm{tra}}(N) & = & {t}_{0}+{NP}+\displaystyle \frac{1}{2}\displaystyle \frac{{dP}}{{dN}}{N}^{2}\\ {t}_{\mathrm{occ}}(N) & = & {t}_{0}+\displaystyle \frac{P}{2}+{NP}+\displaystyle \frac{1}{2}\displaystyle \frac{{dP}}{{dN}}{N}^{2}.\end{array}\end{eqnarray} \tag{ 4 }$$

The third model assumes the planet has a nonzero eccentricity e and its argument of pericenter ω is precessing uniformly: (Giménez & Bastero 1995):

$$\begin{eqnarray}\begin{array}{rcl}{t}_{\mathrm{tra}}(N) & = & {t}_{0}+{{NP}}_{s}-\displaystyle \frac{{{eP}}_{a}}{\pi }\cos \omega (N)\\ {t}_{\mathrm{occ}}(N) & = & {t}_{0}+\displaystyle \frac{{P}_{a}}{2}+{{NP}}_{s}+\displaystyle \frac{{{eP}}_{a}}{\pi }\cos \omega (N)\\ \omega (N) & = & {\omega }_{0}+\displaystyle \frac{d\omega }{{dN}}N\\ {P}_{s} & = & {P}_{a}\left(1-\displaystyle \frac{1}{2\pi }\displaystyle \frac{d\omega }{{dN}}\right),\end{array}\end{eqnarray} \tag{ 5 }$$

where P_s is the sidereal period and P_a is the anomalistic period.

In all three cases, we found the best-fitting model parameters by minimizing χ². We again used the emcee code to perform an MCMC sampling of the posterior distribution in parameter space, given broad uniform priors on all the parameters. We ran the MCMC with 100 walkers, discarding the first 40% of the steps as burn-in and running the code for >10 autocorrelation times. We also double-checked convergence by inspecting the posteriors and computing the Geweke scores for each chain (Geweke 1992). Table 6 gives the fit results.

Table 6.  Timing Model Fit Parameters

Parameter	Value (Uncertainty)
Constant Period Model
Period, P (days)	1.091419649(25)
Midtransit Time of Reference Orbit, t0	2456305.455521(26)
Ndof	156
${\chi }_{\min }^{2}$	380.745
BIC	390.871
Orbital Decay Model
Period, P (days)	1.091420107(42)
Midtransit Time of Reference Orbit, t0	2456305.455809(32)
Decay Rate, dP/dN (days/orbit)	−10:04(69)×10−10
Ndof	155
${\chi }_{\min }^{2}$	167.566
BIC	182.754
Apsidal Precession Model
Sidereal Period, Ps (days)	1.091419633(81)
Midtransit Time of Reference Orbit, t0	2456305.45488(12)
Eccentricity, e	0.00310(35)
Argument of Periastron, ω0 (rad)	2.62(10)
Precession Rate, dω/dN (rad/orbit)	0.000984(+70, −61)
Ndof	153
${\chi }_{\min }^{2}$	179.700
BIC	205.013

Note. Uncertainties in parentheses are the 1σ confidence intervals in the last two digits.

Download table as:  ASCIITypeset image

As was already shown by Maciejewski et al. (2016) and Patra et al. (2017), the constant period model does not fit the data. The minimum value of χ² is 380.7 with 156 degrees of freedom. Figure 4 shows the residuals. Also plotted are the best-fitting model curves for the orbital decay and apsidal precession models. The best-fitting orbital decay model has ${\chi }_{\min }^{2}=167.6$, while the best-fitting apsidal precession model has ${\chi }_{\min }^{2}=179.7$. Thus, while both models fit the data much better than the constant-period model, the orbital decay model is preferred. The difference in χ² is 12.1. Patra et al. (2017) also found a preference for orbital decay, but with a weaker statistical significance (${\rm{\Delta }}{\chi }^{2}=5.5$). Most of the increase in ${\rm{\Delta }}{\chi }^{2}$ is from the newest Spitzer observations of occultations, for which the midpoints are consistent with the predictions of the orbital decay model, but occurred earlier than would be expected based on the apsidal precession model.

Zoom In Zoom Out Reset image size

Our confidence that orbital decay is a better description of the data is enhanced by the fact that the orbital decay model has only three free parameters while the apsidal precession model has five free parameters. A commonly used way to reward a model for fitting the data with fewer free parameters is to compare models with the Bayesian Information Criterion (BIC; Schwarz 1978):

$$\begin{eqnarray}&&\mathrm{BIC}={\chi }^{2}+k\mathrm{log}n,\end{eqnarray} \tag{ 6 }$$

where k is the number of free parameters, and n is the number of data points. In this case, the BIC favors the orbital decay model by Δ(BIC) = 22.3. The interpretation of this number is not completely straightforward, but if we assume the posterior distribution of all the parameters to be a multivariate Gaussian function, then there is a simple relation between ΔBIC and the Bayes factor B:

$$\begin{eqnarray}&&B=\exp \left[-{\rm{\Delta }}(\mathrm{BIC})/2\right]={\rm{70,000}},\end{eqnarray} \tag{ 7 }$$

representing an overwhelming preference for the orbital decay model.

5.2.1. Rømer Effect

If the center of mass of the star–planet system is accelerating along the line of sight with a magnitude ${\dot{v}}_{r}$, then the apparent period of the hot Jupiter would be observed to change due to the Rømer effect:

$$\begin{eqnarray*}&&\displaystyle \frac{\dot{P}}{P}=\displaystyle \frac{{\dot{v}}_{r}}{c}.\end{eqnarray*}$$

If we insert the measured values of P and $\dot{P}$ for WASP-12 into this equation, the implied acceleration is ${\dot{v}}_{r}=0.25$ m s⁻¹ day⁻¹. This is more than an order of magnitude larger than the 2σ upper limit of 0.019 m s⁻¹ day⁻¹ that Patra et al. (2017) obtained by fitting the previously available radial-velocity data. Here, we use the newly obtained radial-velocity data to place an even more stringent upper limit.

In addition to the HIRES radial velocities from Knutson et al. (2014) and described in Section 4, we analyzed the data obtained with the SOPHIE spectrography by Hebb et al. (2009) and Husnoo et al. (2011), as well as data obtained with the HARPS-N spectrograph by Bonomo et al. (2017). We allowed for a constant velocity offset and a "jitter" value specific to each spectrograph. We also excluded the data points obtained during transits, to avoid having to model the Rossiter–McLaughlin effect. We used the radvel code (Fulton et al. 2018) to fit a model consisting of the sum of a sinusoidal function (representing the circular orbit of the planet) and a linear function of time (representing a constant radial acceleration). We fixed the period and time of conjunction to the values from the best-fitting constant period timing model (Table 6).

Figure 5 shows the residuals, after subtracting the best-fitting model. Any line-of-sight acceleration must have an amplitude $| \dot{{v}_{r}}| \lt 0.005$ ${\rm{m}}\,{{\rm{s}}}^{-1}\,{\mathrm{day}}^{-1}$, at the 2σ level. This is a factor of four improvement over the constraints reported by Patra et al. (2017), and two orders of magnitude smaller than the value that would be observed if the observed period derivative were entirely due to the Rømer effect. Thus, we can dismiss the Rømer effects as a significant contributor to the observed period change.¹¹

Zoom In Zoom Out Reset image size

5.2.2. Orbital Eccentricity

5.2.3. Joint RV-timing Fit

In principle, the radial-velocity signal should also be affected by either orbital decay or apsidal precession. A change in the orbital period would affect the RV signal via the true anomaly ν. We computed the radial-velocity signal of a planet with a constant and decaying period using the parameters in Table 6. Over the 9 yr timespan of the HIRES radial-velocity observations, we found that the maximum deviation between the two models is $\sim 10\,{\rm{m}}\,{{\rm{s}}}^{-1}$.

As for apsidal precession, Csizmadia et al. (2019) presented a formula for the associated RV signal:

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{RV} & = & K\left[e\cos \omega (t)+\cos \left(\nu +\omega (t\right)\right.\\ & & +\left.\displaystyle \frac{\dot{\omega }}{n}\displaystyle \frac{{\left(1-{e}^{2}\right)}^{3/2}\cos \left(\nu +\omega (t\right)}{1+e\cos \nu }\right]\\ \omega (t) & = & {\omega }_{0}+\dot{\omega }\left(t-{t}_{0}\right),\end{array}\end{eqnarray} \tag{ 8 }$$

where $n\equiv 2\pi /P$ is the mean motion. Based on this equation, along with the precession period from the best-fitting apsidal precession model, we found that the maximum deviation between the RV signal of a precessing orbit with parameters in Table 6 and a circular model would be $\sim 6\,{\rm{m}}\,{{\rm{s}}}^{-1}$ over a decade. We note that Csizmadia et al. (2019) claimed that apsidal precession would result in residuals on the order of K_orb; however, that only holds if the anomalistic period could be determined independently. In reality, the anomalistic period would need to be determined using the same data, and the residuals would be on the order of e K_orb.

In both cases, there should be small but potentially measurable effects on the RV data. We tried fitting the timing and RV data sets jointly, but the results were essentially unchanged. In particular, the RV data did not alter the ${\rm{\Delta }}\mathrm{BIC}$ between the orbital decay and apsidal precession models. This is likely because the RV observations at later times are sparsely sampled, preventing these small deviations from being measured effectively. Furthermore, the RV jitter of WASP-12 in the HIRES observations is ${\sigma }_{\mathrm{jit}}\sim 10\,{\rm{m}}\,{{\rm{s}}}^{-1}$, a similar magnitude to the decay or precession RV signal. Future RV measurements could be helpful, but only if the RV systematic effects are better understood, including the tidal effect described in Section 5.2.2.

The possibility discussed in the Introduction is that the angular momentum is being transferred directly to the star through the gravitational torque on the star's tidal bulge, and the energy is being dissipated inside the star as the tidal oscillations are converted into heat. In the "constant phase lag" model of Goldreich & Soter (1966), assuming that the planet's mass stays constant, the decay rate is

$$\begin{eqnarray*}&&\dot{P}=-\displaystyle \frac{27\pi }{2{Q}_{\star }^{{\rm{{\prime} }}}}\left(\displaystyle \frac{{M}_{p}}{{M}_{\star }}\right){\left(\displaystyle \frac{{R}_{\star }}{a}\right)}^{5},\end{eqnarray*}$$

where ${Q}_{\star }^{{\prime} }$ is the star's "modified tidal quality factor," defined as the quality factor ${Q}_{\star }$ divided by 2/3 of the Love number k₂. Inserting the measured decay rate, we obtain for WASP-12 a tidal quality factor of

$$\begin{eqnarray*}&&{Q}_{\star }^{{\prime} }={1.75}_{-0.11}^{+0.13}\times {10}^{5}.\end{eqnarray*}$$

This result for ${Q}_{\star }^{{\prime} }$ is lower than most of the estimates in the literature, which are based on less direct observations. By analyzing the eccentricity distribution of stellar binaries, Meibom & Mathieu (2005) found ${Q}_{\star }^{{\prime} }$ to be in the range from 10⁵ to 10⁷, while similar studies applied to hot Jupiter systems by Jackson et al. (2008) and Husnoo et al. (2012) found ${Q}_{\star }^{{\prime} }={10}^{5.5}\mbox{--}{10}^{6.5}$. Hamer & Schlaufman (2019) found that hot Jupiter host stars are kinematically younger than similar stars without hot Jupiters, and interpreted the result as evidence for the tidal destruction of hot Jupiters, finding ${Q}_{\star }^{{\prime} }$ $={10}^{6}-{10}^{6.5}$. Collier Cameron & Jardine (2018) modeled the orbital period distribution of the hot Jupiter population and found ${Q}_{\star }^{{\prime} }={10}^{7}\mbox{--}{10}^{8}$. Penev et al. (2012) modeled the star/planet tidal interactions in selected systems and also found ${Q}_{\star }^{{\prime} }\gt {10}^{7}$. However, these studies assumed Qs to be a universal constant, even though one would expect it to depend on forcing frequency and perhaps many other parameters. Penev et al. (2018) updated and expanded the approach of system-by-system modeling to allow for frequency dependence, finding that ${Q}_{\star }^{{\prime} }$ ranges from 10⁵ to 10⁷ for orbital periods of 0.5–2 days (Penev et al. 2018).

If tidal dissipation is responsible for the orbital decay of WASP-12, then the dissipation rate is higher than would have been expected according to these earlier studies. The physical mechanism for such rapid dissipation is unclear. Attempts to compute the tidal quality factor from physical principles for either the equilibrium or dynamical tide generally find larger values of ${Q}_{\star }^{{\prime} }$, from 10⁷ to 10¹⁰ (as reviewed by Ogilvie 2014). Weinberg et al. (2017) argued that WASP-12 could be a subgiant star, in which case nonlinear wave-breaking of the dynamical tide near the stellar core would lead to efficient dissipation and ${Q}_{\star }^{{\prime} }\sim 2\times {10}^{5}$. However, Bailey & Goodman (2019) examined this possibility and found that the observed properties of WASP-12 are more compatible with the expected properties of a main-sequence star rather than a subgiant.

Millholland & Laughlin (2018) proposed an alternate hypothesis, in which WASP-12b is in a spin-orbit resonance with an external perturber, allowing it to maintain a large obliquity and giving rise to obliquity tides. The ongoing dissipation of obliquity tides might be strong enough to explain the observed decay rate. While such a hypothetical perturber cannot yet be ruled out, this scenario calls for a specific system architecture with a misaligned perturber just below the limits of detection.

The Roche limit for a close-orbiting planet can be expressed as a minimum orbital period depending on the density distribution of the planet (Rappaport et al. 2013). For WASP-12b, with a mean density of 0.46 g cm⁻³, the Roche-limiting orbital period is 14.2 hr assuming the mass of the planet to be concentrated near the center, and 18.6 hr in the opposite limit of a spherical and incompressible planet. The true orbital period of 26 hr is longer than either of these limits, implying that the planet is not filling its Roche lobe. This simple comparison does not take into account the planet's tidal distortion, which would lead to a longer period for the Roche limit, but probably not as long as 26 hr.

Nevertheless, there may exist an optically thin exosphere or wind that does fill the Roche lobe. Ultraviolet observations of the WASP-12 system indicate that the planet is surrounded by a cloud of absorbing material that is larger than the Roche lobe (Fossati et al. 2010; Haswell et al. 2012; Nichols et al. 2015). Recent infrared phase curve observations of the WASP-12 system also corroborate the idea that gas is being stripped from the planet (Bell et al. 2019). The mass-loss rate is not measured directly (Haswell 2018), but models of this process suggest it could be as high as $\sim 3\times {10}^{14}$ g s⁻¹, corresponding to ${M}_{{\rm{p}}}/\dot{M}\approx 300\,\mathrm{Myr}$ (Lai et al. 2010; Jackson et al. 2017). This mass-loss timescale, while longer than the tidal decay timescale, is still short relative to the age of the star. Putting this evidence together with the observed changes in orbital period, it seems that tidal orbital decay has brought the planet close enough to initiate Roche lobe overflow of the planet's tenuous outer atmosphere.

The escaping mass bears energy and angular momentum, which can lead to changes in the orbital period separate from the effects of the tidal bulge of the star. To assess whether or not the escaping mass bears enough angular momentum to be relevant for period changes, we need to know more about the flow of mass away from the planet. We would expect the gas to flow out from the inner Lagrange point, with lower specific angular momentum than the rest of the planet. As a result, the planet's specific angular momentum would rise, causing its orbit to widen. In this scenario, the rate of period decrease that we have measured would be the result of a competition between tidally driven decay and mass-loss-driven growth. Jia & Spruit (2017) presented an expression relating the change in semimajor axis to the mass-loss rate:

$$\begin{eqnarray}&&\displaystyle \frac{\dot{a}}{a}=\displaystyle \frac{2}{1+q}\displaystyle \frac{{\dot{M}}_{{\rm{p}}}}{{M}_{{\rm{p}}}}\left[\displaystyle \frac{{x}_{{\rm{L}}1}^{2}}{{x}_{{\rm{p}}}^{2}}(1-\epsilon )+{q}^{2}-1\right]\gt 0,\end{eqnarray} \tag{ 11 }$$

where q is the planet–star mass ratio, ${x}_{{\rm{L}}1},{x}_{{\rm{p}}}$ are the distances from the L1 point and planet to the system center of mass, and parameterizes the fraction of angular momentum that is transferred back to the planet from the outflowing gas. Using the mass-loss rate estimated by Lai et al. (2010) and Jackson et al. (2017), and assuming no angular momentum transfer back to the planet ($\epsilon =0$), this yields

$$\begin{eqnarray}&&\displaystyle \frac{\dot{P}}{P}\approx 3\displaystyle \frac{{\dot{M}}_{{\rm{p}}}}{{M}_{{\rm{p}}}}\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&\dot{P}\approx 1\,\mathrm{ms}\,{\mathrm{yr}}^{-1}.\end{eqnarray} \tag{ 13 }$$

Under these assumptions, the mass loss from the planet causes a period increase about 30 times smaller than the observed period change for WASP-12b, and is unlikely to be slowing down the planet's orbital decay.

On the other hand, if gas flows out of the L2 point, there could be a net loss in angular momentum from the planet, hastening any tidally driven decay of the orbit. Valsecchi et al. (2015) examined this process and found that the orbital period begins to shrink only in the final stages of mass loss, when the remaining mass of the planet is dominated by the dense core. This is not the current situation of WASP-12b, given its Jovian mass and low mean density. Furthermore, Jia & Spruit (2017) found that even in such a scenario, the mass loss from L1 continues to dominate over the outflow from L2, leading to continued expansion of the orbit. In either case, it seems unlikely that the mass loss from WASP-12b is contributing significantly to its orbital evolution.