WASP-4b Arrived Early for the TESS Mission

WASP-4 was observed by TESS with Camera 2 from 2018 August 23 to September 20, within the second "sector" of science operations. The star is designated as TIC 402026209 in the TESS Input Catalog (Stassun et al. 2018). The pixel data for an 11 × 11 array surrounding WASP-4 were averaged into 2 min stacks by the onboard computer. The data were downlinked via the Deep Space Network,¹⁴ and the spacecraft time stamps were calibrated against the ground-station clocks. The spacecraft clock times were then transformed by the Payload Operations Center into the Temps Dynamique Barycentrique (TDB) reference system. The images were then reduced to light curves by the Science Processing Operations Center (SPOC) at NASA Ames (Jenkins et al. 2016). During this processing, the SPOC used the known spacecraft trajectory to compute the barycentric time corrections on a target-by-target basis and expressed the time stamps as Barycentric TESS Julian Dates (BTJD), which is simply the Barycentric Julian Date minus 2,457,000. Light curves that were flagged by the SPOC pipeline as crossing a transit detection threshold were then vetted and released by the MIT TESS Science Office to the Mikulski Archive for Space Telescopes on 2018 November 29 (Ricker & Vanderspek 2018).

We began our analysis with the Presearch Data Conditioning (PDC) light curve, which has had nonastrophysical variability removed through the methods discussed by Smith et al. (2017a, 2017b). We then processed the light curve as follows. First, we removed all points with nonzero quality flags. This removed data that might have been adversely affected by "momentum dumps," the firing of thrusters and resetting of reaction wheels¹⁵ that took place every 2.5 days during sector 2.   The data during these events were assigned quality flags corresponding to "Reaction Wheel Desaturation Event" and "Manual Exclude" (Tenenbaum & Jenkins 2018, Table 28). For WASP-4, these flags were simultaneously set for 54 distinct cadences, and there were 10 momentum dumps, averaging about 10 min of flagged data per dump. Out of caution, we clipped out an additional 10 min before and after every momentum dump. We also removed the data within the first and last hours of both orbits, because of correlated red noise that appears during those time ranges.

All told, we removed 8% of the original data points and were left with 18,165 measurements of the relative flux of WASP-4. We normalized the data by dividing out the median flux, yielding the measurements shown in Figure 1. We converted the time stamps from BTJD_TDB into BJD_TDB by adding the appropriate 2,457,000 day offset (Tenenbaum & Jenkins 2018). Many of these and subsequent processing steps were performed using astrobase (Bhatti et al. 2018). We did not "flatten" the light curves, as is often done with splines, polynomials, or Gaussian processes. Instead, we modeled the out-of-transit flux variations simultaneously with the transit parameters, as described below.

Using the cleaned PDC light curve, we applied the Box Least Squares algorithm (Kovács et al. 2002) to estimate the orbital period, transit duration, and a reference epoch using the TESS data alone. Based on the results, we isolated the data within four transit durations of each transit midpoint. To find the transit parameters that best fit the data, we first fitted a line to the out-of-transit flux measurements surrounding each transit and divided it out. We then created a phase-folded light curve from all 18 transits and fitted a standard transit model using the analytic formulae given by Mandel & Agol (2002) and implemented by Kreidberg (2015, BATMAN). We assumed the orbit to be circular, consistent with the limits from radial velocities and occultation timing (Beerer et al. 2011; Knutson et al. 2014; Bonomo et al. 2017). The free parameters were the reference epoch, the planet-to-star radius ratio R_p/R_⋆, the orbital distance to stellar radius ratio a/R_⋆, the inclination i, two quadratic limb-darkening coefficients $({u}_{\mathrm{linear}},{u}_{\mathrm{quad}})$, and the orbital period P.

We sampled the posterior probability distribution for all of the parameters using the algorithm proposed by Goodman & Weare (2010) and implemented by Foreman-Mackey et al. (2013, emcee). Table 1 gives the results, which are in reasonable agreement with the parameters reported by, for example, Southworth et al. (2009) and Huitson et al. (2017). Figure 2 shows the phase-folded light curve.

Table 1.  Selected System Parameters of WASP-4b

Parameter	Value	68% Confidence Interval	Comment
Transit/RV parameters:
Rp/R⋆	0.15201	$+0.00040$, −0.00033	A
i (deg)	89.06	+0.65, −0.84	A
a/R⋆	5.451	+0.023, −0.052	A
${u}_{\mathrm{linear}}$	0.382	⋯	A
${u}_{\mathrm{quad}}$	0.210	⋯	A
K (m s−1)	241.1	+2.8, −3.1	B
Stellar parameters:
${T}_{\mathrm{eff}}$ (K)	5400	±90	C
$\mathrm{log}{g}_{\star }$ (cgs)	4.47	±0.11	C
$[\mathrm{Fe}/{\rm{H}}]$	−0.07	±0.19	C
${F}_{\mathrm{bol}}$ (erg cm−2 s−1)	$2.802\times {10}^{-10}$	$\pm 0.076\times {10}^{-10}$	D
AV (mag)	0.03	+0.02, −0.01	D
π (mas)	3.7145	0.0517	F
R⋆ (R⊙)	0.893	±0.034	E
${\rho }_{\star }$ (g cm−3)	1.711	$+0.022$, −0.048	E
${M}_{\star }$ (M⊙)	0.864	$+0.084$, −0.090	E
T magnitude	11.778	±0.018	G
Planetary parameters:
a (au)	0.0226	$+0.0007$, −0.0008	E
Mp (MJup)	1.186	$+0.090$, −0.098	E
Rp (RJup)	1.321	±0.039	E

Note. (A) From phase-folded TESS light curve (Section 2.2). Orbital periods are in Table 4. The limb-darkening parameters were allowed to float around the Claret (2017) prediction, but were unconstrained. (B) Triaud et al. (2010). (C) From HARPS spectra (Doyle et al. 2013). (D) Stassun et al. (2017). (E) This work, see Section 2.3. (F) Gaia Collaboration et al. (2018). (G) Stassun et al. (2018).

Download table as:  ASCIITypeset image

To measure the transit times, we returned to the "cleaned" PDC time series and fitted the data within four transit durations of each transit separately. We used four free parameters, the time of midtransit t_tra, the planet-to-star radius ratio, and the slope and intercept of a linear trend, to account for any slow variations unrelated to the transit. We fixed the remaining parameters at the values that had been determined from the phase-folded TESS light curve. The uncertainty in each photometric data point was set equal to the rms level of the out-of-transit data.

To verify that the measured uncertainties are estimated accurately, we computed the χ² value for a linear ephemeris fit to the measured TESS midtransit times. We found that ${\chi }^{2}=9.2$, with n = 16 degrees of freedom. The variance of the χ² distribution is 2n, so we would expect χ² = 16 ± 5.7. Visually inspecting the residuals showed that the error variance had been overestimated, so we multiplied the measured TESS errors by a factor f = 0.76, forcing a reduced χ² of unity. This lowered the mean uncertainty of the transit midtimes from 29.8 to 22.6 s. We verified that omitting this step did not appreciably alter any of our conclusions.

Figure 1 shows the light curve of each individual transit, the best-fit models, and the residuals. Table 2 reports the midtransit times and their uncertainties. After binning the residuals to 1 hr windows, the light curves have an rms scatter of 586 ppm. The prelaunch TESS noise model¹⁶ (Winn 2013; Sullivan et al. 2015, Section 6.4) would have predicted an error budget consisting of the following terms added in quadrature: 410 ppm from photon-counting noise, 202 ppm from detector read noise, and 673 ppm from the zodiacal background light. The level of background light appears to have been overestimated in the model.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2.  WASP-4b Transit Times, Uncertainties, and References

ttra (BJDTDB)	${\sigma }_{{t}_{\mathrm{tra}}}$ (days)	Epoch	H13?	References
2454368.59279	0.00033	−1073	1	Wilson et al. (2008)
2454396.69576	0.00012	−1052	1	Gillon et al. (2009b)
2454697.79817	0.00009	−827	1	Winn et al. (2009)
2454701.81303	0.00018	−824	1	Hoyer et al. (2013)
2454701.81280	0.00022	−824	1	Hoyer et al. (2013)
2454705.82715	0.00029	−821	1	Hoyer et al. (2013)
2454728.57767	0.00042	−804	1	Hoyer et al. (2013)
2454732.59197	0.00050	−801	1	Hoyer et al. (2013)
2454740.62125	0.00035	−795	1	Hoyer et al. (2013)
2454748.65111	0.00007	−789	1	Winn et al. (2009)
2454752.66576	0.00069	−786	1	Dragomir et al. (2011)
2455041.72377	0.00018	−570	1	Hoyer et al. (2013)
2455045.73853	0.00008	−567	1	Sanchis-Ojeda et al. (2011)
2455049.75325	0.00007	−564	1	Sanchis-Ojeda et al. (2011)
2455053.76774	0.00009	−561	1	Sanchis-Ojeda et al. (2011)
2455069.82661	0.00029	−549	1	Nikolov et al. (2012)
2455069.82617	0.00038	−549	1	Nikolov et al. (2012)
2455069.82670	0.00028	−549	1	Nikolov et al. (2012)
2455069.82676	0.00031	−549	1	Nikolov et al. (2012)
2455073.84108	0.00029	−546	1	Nikolov et al. (2012)
2455073.84128	0.00026	−546	1	Nikolov et al. (2012)
2455073.84111	0.00023	−546	1	Nikolov et al. (2012)
2455073.84114	0.00018	−546	1	Nikolov et al. (2012)
2455085.88418	0.00086	−537	1	Dragomir et al. (2011)
2455096.59148	0.00022	−529	1	Hoyer et al. (2013)
2455100.60595	0.00012	−526	1	Sanchis-Ojeda et al. (2011)
2455112.64986	0.00039	−517	1	Nikolov et al. (2012)
2455112.65009	0.00033	−517	1	Nikolov et al. (2012)
2455112.65005	0.00031	−517	1	Nikolov et al. (2012)
2455112.65005	0.00049	−517	1	Nikolov et al. (2012)
2455132.72310	0.00041	−502	1	Hoyer et al. (2013)
2455468.61943	0.00046	−251	1	Hoyer et al. (2013)
2455526.16356	0.00008	−208	0	Ranjan et al. (2014)
2455828.60375	0.00041	18	1	Hoyer et al. (2013)
2455832.61815	0.00041	21	1	Hoyer et al. (2013)
2455844.66287	0.00009	30	0	Huitson et al. (2017)
2456216.69123	0.00006	308	0	Huitson et al. (2017)
2456288.95622	0.00015	362	0	C. Baxter et al. (2019, in preparation)
2456292.97025	0.00019	365	0	C. Baxter et al. (2019, in preparation)
2456576.67556	0.00005	577	0	Huitson et al. (2017)
2456924.61561	0.00006	837	0	Huitson et al. (2017)
2458355.18490	0.00025	1906	0	This work
2458356.52251	0.00027	1907	0	This work
2458357.86105	0.00026	1908	0	This work
2458359.19946	0.00026	1909	0	This work
2458360.53707	0.00028	1910	0	This work
2458361.87538	0.00025	1911	0	This work
2458363.21411	0.00027	1912	0	This work
2458364.55193	0.00025	1913	0	This work
2458365.89057	0.00026	1914	0	This work
2458369.90506	0.00028	1917	0	This work
2458371.24298	0.00026	1918	0	This work
2458372.58124	0.00026	1919	0	This work
2458373.91981	0.00028	1920	0	This work
2458375.25792	0.00025	1921	0	This work
2458376.59623	0.00024	1922	0	This work
2458377.93434	0.00026	1923	0	This work
2458379.27319	0.00025	1924	0	This work
2458380.61098	0.00028	1925	0	This work

Note. t_tra is the measured transit midtime, and ${\sigma }_{{t}_{\mathrm{tra}}}$ is its 1σ uncertainty. ${\sigma }_{{t}_{\mathrm{tra}}}$ was evaluated from the sampled posteriors by taking the maximum of the difference between the 84th percentile minus the median, and the median minus the 16th percentile. The resulting error variances then appeared to have been overestimated, so we lowered the uncertainties as described in Section 2.2. The "Reference" column refers to the work describing the original observations. The "H13?" column is 1 if the midtime value was taken from Hoyer et al. (2013). Otherwise, the midtime came from the column listed in "Reference." The Hoyer et al. (2013) BJD_TT times are equal to BJD_TDB for our purposes (Urban & Seidelmann 2012). We omitted the timing measurements from Southworth et al. (2009) because there were technical problems with the computer clock at the time of observation (Nikolov et al. 2012). The two C. Baxter et al. (2019, in preparation) times were obtained from Spitzer/IRAC transit light curves in the 3.6 and 4.5 μm channels.

Download table as:  ASCIITypeset images: 1 2

We calculated the stellar and planetary parameters in the following way. We computed the star's spectral energy distribution based on the Gaia DR2 parallax (after making the small correction advocated by Stassun & Torres 2018) and the broadband magnitudes from the available all-sky catalogs: G from Gaia DR2, B_T and V_T from Tycho-2, BVgri from APASS, JHK_S from 2MASS, and the WISE 1–4 passbands, thus spanning the wavelength range 0.4–22  μm. We adopted the effective temperature from the work by Doyle et al. (2013), who determined the spectroscopic parameters of WASP-4 using high signal-to-noise ratio observations with the High Accuracy Radial-velocity Planet Searcher (HARPS). Then, we determined the stellar radius through the combination of the bolometric luminosity and the effective temperature, using the Stefan–Boltzmann law. To determine the stellar mass, we first computed the mean stellar density based on the value of a/R_⋆ that gave the best fit to the phase-folded TESS light curve (for the relevant equation, see Seager & Mallén-Ornelas 2003 or Winn 2010). The mass was calculated from the radius and density, and the orbital distance was also calculated from the radius and a/R_⋆. The planetary radius was calculated as the product of R_⋆ and R_p/R_⋆. Finally, the planet mass was calculated based on the stellar mass, the radial-velocity amplitude observed by Triaud et al. (2010), and the orbital inclination.

Table 1 gives the resulting parameters, which we adopted for the remaining analysis. The uncertainties in our derived stellar and planetary parameters are propagated according to standard analytic formulae, under the assumption that the variables are uncorrelated and normally distributed. Our system parameters are in agreement with those of previous investigators, but have the benefit of incorporating the Gaia parallax (Wilson et al. 2008; Gillon et al. 2009a; Winn et al. 2009; Southworth 2011; Petrucci et al. 2013; Huitson et al. 2017). By comparing the star's luminosity and spectroscopic parameters with the outputs of the Yonsei–Yale stellar-evolutionary models, we found that WASP-4 is a main-sequence star with an age of approximately 7 Gyr (Demarque et al. 2004).

Table 2 gives the transit times we used in our analysis. We included data from peer-reviewed literature for which the analysis was based on observations of a single transit, and for which the midpoint was allowed to be a free parameter. We also required that the time system be clearly documented. Many of the times were previously compiled by Hoyer et al. (2013). We confirmed that the times in that paper were in agreement with the original sources and that barycentric corrections had been performed when needed.

The earliest epoch is from EulerCam on the 1.2 m Euler telescope (Wilson et al. 2008). The second epoch is based on z-band photometry acquired by Gillon et al. (2009b) at the VLT 8.2 m with FORS2. Subsequent observations were performed by Winn et al. (2009), Dragomir et al. (2011), Sanchis-Ojeda et al. (2011), Nikolov et al. (2012), Hoyer et al. (2013), and Ranjan et al. (2014). Finally, Huitson et al. (2017) acquired optical transit spectra with the 8.1 m Gemini South telescope between 2011 and 2014, one transit per season. The per-point standard deviation of their light curves was a few hundred parts per million. The average precision in their reported transit times was 5.6 s. Since these data points carry significant weight in the analysis, we corresponded with the authors to confirm that the time stamps in their data represent midexposure times, that the barycentric correction was performed correctly, and that the time system of the final results was BJD_TDB. These same authors also used the same instrument and method to analyze other hot Jupiters, none of which showed a departure from a constant-period model. Finally, these authors also measured two transit midpoints using Spitzer (C. Baxter et al. 2019, in preparation); these times agree with the Gemini South results.

We also compiled the available occultation times, which are given in Table 3. The tabulated values have been corrected for the light travel time across the diameter of the orbit by subtracting 2a/c = 22.8 s from the observed time. Beerer et al. (2011) observed two occultations of WASP-4b using warm Spitzer in the 3.6 and 4.5 μm bands. Cáceres et al. (2011) detected an occultation from the ground in the K_S band and gave a time in HJD, without specifying the time standard. We assumed the standard was UTC, and we performed the appropriate corrections to convert to BJD_TDB. We verified that none of our conclusions would be changed if this assumption was mistaken. Finally, Zhou et al. (2015) observed an occultation with the Anglo-Australian Telescope. They did not report the observed midpoint, but they did report a result for $e\cos \omega$ based upon the observed midpoint. We calculated the implied midpoint using the formula (e.g., Winn 2010)

$$\begin{eqnarray}&&{t}_{\mathrm{occ}}(E)={t}_{0}+{PE}+\displaystyle \frac{P}{2}\left(1+\displaystyle \frac{4}{\pi }e\cos \omega \right),\end{eqnarray} \tag{ 1 }$$

for E the transit number, t₀ the reference epoch, e the eccentricity, and ω the argument of pericenter. In total, there are four available occultation times.

Table 3.  WASP-4b Occultation Times, Uncertainties, and References

${t}_{\mathrm{occ}}$ (BJDTDB)	${\sigma }_{{t}_{\mathrm{occ}}}$ (days)	Epoch	References
2455102.61210	0.00074	−511	Cáceres et al. (2011)a
2455172.20159	0.00130	−459	Beerer et al. (2011)
2455174.87780	0.00087	−457	Beerer et al. (2011)
2456907.88714	0.00290	838	Zhou et al. (2015)b

Notes. ${t}_{\mathrm{occ}}$ is the measured occultation midtime, minus the 2a/c = 22.8 s light travel time; ${\sigma }_{{t}_{\mathrm{occ}}}$ is the 1σ uncertainty on the occultation time.

^aCáceres et al. (2011) reported this time in "HJD," with an unspecified time standard. We assumed the time was originally in ${\mathrm{HJD}}_{\mathrm{UTC}}$ and converted to ${\mathrm{BJD}}_{\mathrm{TDB}}$ for the tabulated time. ^bZhou et al. (2015) fixed the epoch and let $e\cos \omega$ float. Using the reported dates of observation, we converted their $e\cos \omega$ values into an occultation time using Equation (1) of the text.

Download table as:  ASCIITypeset image

First, we performed a weighted least-squares fit of a constant-period model (a "linear ephemeris") to the pre-TESS data, and we used it to extrapolate to the epochs of the TESS observations. The residuals of the best-fitting model are shown in Figure 3. The transits observed by TESS occurred earlier than expected. Because the TESS mission is still in an early stage, we were concerned about a possible offset in the TESS time stamps due to an error with the TESS clock or the data processing pipeline. The appendix describes some tests that convinced us that a simple offset is unlikely. Assuming that the observed timing variation is astrophysical, we proceeded by exploring three models for the timing data in a manner identical to the study by Patra et al. (2017).

Zoom In Zoom Out Reset image size

The first model assumes a constant orbital period on a circular orbit:

$$\begin{eqnarray}&&{t}_{\mathrm{tra}}(E)={t}_{0}+{PE},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{t}_{\mathrm{occ}}(E)={t}_{0}+\displaystyle \frac{P}{2}+{PE},\end{eqnarray} \tag{ 3 }$$

where E is the epoch number. We defined the epoch numbers such that E = 0 is near the weighted average of the observed times. This helps to reduce the covariance between t₀ and P.

The second model assumes the period is changing at a steady rate:

$$\begin{eqnarray}&&{t}_{\mathrm{tra}}(E)={t}_{0}+{PE}+\displaystyle \frac{1}{2}\displaystyle \frac{{dP}}{{dE}}{E}^{2},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{t}_{\mathrm{occ}}(E)={t}_{0}+\displaystyle \frac{P}{2}+{PE}+\displaystyle \frac{1}{2}\displaystyle \frac{{dP}}{{dE}}{E}^{2}.\end{eqnarray} \tag{ 5 }$$

The three free parameters are the reference epoch t₀, the period at the reference epoch, and the period derivative, ${dP}/{dt}=(1/P){dP}/{dE}$.

The third model assumes the planet has a slightly eccentric orbit and that the line of apsides is rotating (Giménez & Bastero 1995):

$$\begin{eqnarray}&&{t}_{\mathrm{tra}}(E)={t}_{0}+{P}_{{\rm{s}}}E-\displaystyle \frac{{{eP}}_{{\rm{a}}}}{\pi }\cos \omega ,\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{t}_{\mathrm{occ}}(E)={t}_{0}+\displaystyle \frac{{P}_{{\rm{a}}}}{2}+{P}_{{\rm{s}}}E+\displaystyle \frac{{{eP}}_{{\rm{a}}}}{\pi }\cos \omega ,\end{eqnarray} \tag{ 7 }$$

where P_s is the sidereal period, e is the eccentricity, ${P}_{{\rm{a}}}$ is the anomalistic period, and ω is the argument of pericenter. In this model, the angular velocity of the line of apsides $d\omega /{dE}$ is constant,

$$\begin{eqnarray}&&\omega (E)={\omega }_{0}+\displaystyle \frac{d\omega }{{dE}}E,\end{eqnarray} \tag{ 8 }$$

and the sidereal and anomalistic periods are connected through the equation

$$\begin{eqnarray}&&{P}_{{\rm{s}}}={P}_{{\rm{a}}}\left(1-\displaystyle \frac{1}{2\pi }\displaystyle \frac{d\omega }{{dE}}\right).\end{eqnarray} \tag{ 9 }$$

The sidereal period is the duration required to return to the same orientation with respect to the stars; the slightly longer anomalistic period is the duration required to reach a fixed longitude with respect to the rotating line of apsides. The five free parameters of this model are $({t}_{0},{P}_{{\rm{s}}},e,{\omega }_{0},d\omega /{dE})$, denoting the reference epoch, the sidereal period, the eccentricity, the argument of pericenter at the reference epoch, and the angular velocity of the line of apsides.

We fitted each model by assuming a Gaussian likelihood and sampling over the posterior probability distributions. The prior for the quadratic model allowed the period derivative to have any sign. We considered two possible priors for the precession model: the first is a wide prior that allows nonphysical values of the planetary Love number (Equation (18)). The second prior requires that the planetary Love number is less than that of a sphere of constant density.

Figure 4 shows the residuals with respect to the constant-period model. The best-fitting constant-period model has χ² = 174 and 61 degrees of freedom. The best-fitting quadratic model has χ² = 62.6 and 60 degrees of freedom. The best-fitting precession model under the wide prior has χ² = 64.3 and 58 degrees of freedom. Under the physical prior, the best fit is slightly worse, with χ² = 64.8. The precession and quadratic models both provide much better fits than the constant-period model. The difference in χ² between the linear and quadratic models corresponds to $p\approx {10}^{-26}$.

Zoom In Zoom Out Reset image size

The quadratic model provides a marginally better fit to the data than the precession model. It is favored by Δχ² = 1.7 and has two fewer free parameters. A useful heuristic for model comparison is the Bayesian information criterion,

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

where k is the number of free parameters, and n is the number of data points. In this case, n = 62. The difference in the BIC between the precession and decay models is ${\rm{\Delta }}\mathrm{BIC}\,={\mathrm{BIC}}_{\mathrm{prec}}-{\mathrm{BIC}}_{\mathrm{quad}}=10$,  assuming a wide prior for the precession model. This corresponds to a Bayes factor of $1.1\times {10}^{4}$. Likewise, the Akaike information criterion favors the constant-period-derivative model by ΔAIC = 5.8. Differences of this magnitude are traditionally deemed "strong evidence" that one model is a better description of the data than the other (Kass & Raftery 1995), although we prefer to reserve judgment until more data can be obtained.

In the quadratic model, the period derivative is

$$\begin{eqnarray}&&\dot{P}=-(4.00\pm 0.38)\times {10}^{-10}=-12.6\pm 1.2\,\mathrm{ms}\,{\mathrm{yr}}^{-1}.\end{eqnarray} \tag{ 11 }$$

For comparison, the best-fitting period derivative of WASP-12b is $\dot{P}=-29\pm 3\,\mathrm{ms}\,{\mathrm{yr}}^{-1}$ (Maciejewski et al. 2016; Patra et al. 2017). If both planets are truly falling onto their stars, then WASP-4b is falling at about half the rate of WASP-12b.

In the precession model, the best-fit eccentricity is

$$\begin{eqnarray}&&e=({1.92}_{-0.76}^{+1.93})\times {10}^{-3}.\end{eqnarray} \tag{ 12 }$$

The longitude of periastron advances by $\dot{\omega }\,={13.6}_{-3.6}^{+4.7}\ \mathrm{degrees}\,{\mathrm{yr}}^{-1}$, and the precession period is ${27}_{-7}^{+10}$ yr. All of the best-fitting parameters (the medians of the posterior distributions) and the 68% credible intervals are reported in Table 4.

Table 4.  Best-fit Transit Timing Model Parameters

Parameter	Median Value (Unc.)a
Constant period
t0 (${\mathrm{BJD}}_{\mathrm{TBD}})$	2455804.515752(+19)(−19)
P (days)	1.338231466(+23)(−22)
Constant period derivative
t0 (${\mathrm{BJD}}_{\mathrm{TBD}})$	2455804.515918(+24)(−24)
P (days)	1.338231679(+31)(−31)
dP/dt	$-4.00(+37)(-38)\times {10}^{-10}$
Apsidal precession (wide prior)
t0 (${\mathrm{BJD}}_{\mathrm{TBD}})$	2455804.51530(+25)(−31)
Ps (days)	1.33823127(+20)(−48)
e	${1.92}_{-0.76}^{+1.93}\times {10}^{-3}$
${\omega }_{0}$ (rad)	2.40(+38)(−34)
$d\omega /{dE}$ (rad epoch−1)	${8.70}_{-2.30}^{+3.01}\times {10}^{-4}$

Note.

^aThe numbers in parentheses give the 68% confidence interval for the final two digits, where appropriate.

Download table as:  ASCIITypeset image

To assess the overall robustness of our results, we considered a few possible systematic effects in the timing data set.

Suspect pre-TESS light curves—Some of the pre-TESS light curves have incomplete phase coverage: a handful of the transit times in Table 2 are from light curves with gaps. A separate issue is the effect of spot-crossing anomalies on transit timings. To address these concerns, we repeated the model fitting described above, but omitted epochs −827, −804, −537, and −208 because of gaps in their coverage. We also omitted epochs −526 and −561 because of visible spot anomalies during the transits. (All epoch numbers are as in Table 2.) The resulting best-fit transit timing model parameters were all within 1σ of the values quoted in Table 4. The uncertainties, goodness-of-fit statistics, and model comparison statistics did not appreciably change.

Spot-crossing events in TESS data—To explore the effect of possible spot-crossing events on the TESS transit time measurements, we performed a separate test. We injected triangular spot anomalies with amplitude 0.03% and duration 30 min at random phases into each transit. The amplitude was chosen to be larger than the spot-crossing anomalies observed by Southworth et al. (2009) and Sanchis-Ojeda et al. (2011), and the duration was chosen to be comparable to those of previously observed events. Spots of these amplitudes resemble the possible anomalies present in transits "1360.54" and "1372.58" of Figure 1.

With spots injected, we repeated our measurement of the transit times. On average, the measured transit times did not change after injecting spots, because the flux deviations are equally likely to occur in the first and second halves of the transit. For individual transits, there were no cases for which the timing deviation was larger than one minute. The largest shifts occur when the spot anomaly occurs during transit ingress or egress, in which case the measured midtime is shifted either late or early by between 30 and 50 s.

The TESS observations could therefore be skewed early if there were spot-crossing events during every egress. Two arguments rule out this possibility. (1) The light-curve residuals do not show evidence for these events. (2) The stellar rotation period is between 20 and 40 days, and the sky-projected stellar obliquity is less than 10° (Triaud et al. 2010; Sanchis-Ojeda et al. 2011; Hoyer et al. 2013). Since the planet orbits every 1.3 days, requiring that spot anomalies always occur during egress would be equivalent to requiring a stellar spot distribution that is exquisitely (and thus implausibly) distributed to match the planet egress times.

Detrending choices in pre-TESS data—There is a final concern that is difficult to address. We collected the midtransit time values derived by different authors, who used heterogeneous methods to fit and detrend their light curves. We have also assumed that these authors have correctly documented the time systems in which the data are reported. Further, though many choices in transit fitting (e.g., parameterization of limb darkening and eccentricity) do not affect transit midtime measurements, different detrending approaches can asymmetrically warp transits and shift midtransit times. The magnitude of this systematic effect is hard to quantify, but the situation is fairly clear from Figure 4. Many independent authors provided transit measurements shortly after WASP-4b's discovery, and the data are consistent with each other. Huitson et al. (2017) provided the most important data from epochs 0 to 1000. If their data were systematically affected by detrending choices or time-system confusion at the level of several times their reported uncertainties, then it is possible that the orbital period is constant despite the evidence in the TESS data. For this reason, we paid careful attention to the Huitson et al. (2017) data set, and we corresponded with the authors to confirm that their results are not affected by systematic effects of the required amplitude.

None of the concerns mentioned in this subsection seem likely to explain the observed timing variations. We proceed by considering possible astrophysical explanations.

If the timing variation is caused entirely by orbital decay, then the best-fit model parameters yield a characteristic decay timescale of

$$\begin{eqnarray}&&\displaystyle \frac{P}{{dP}/{dt}}=9.2\,\mathrm{Myr}.\end{eqnarray} \tag{ 13 }$$

For comparison, the corresponding time for WASP-12b is 3.2 Myr (Patra et al. 2017).

If WASP-4 really is undergoing rapid orbital decay, then how many of the other known hot Jupiters should have orbits that are decaying at detectable (or nearly detectable) rates? Figure 5 compares some key properties of WASP-4 with those of a larger ensemble of hot Jupiters. The middle panel displays two parameters that strongly affect the expected orbital-decay timescale, ${{PM}}_{\star }/{M}_{{\rm{p}}}$ and a/R_⋆. WASP-4 has one of the shortest theoretical timescales for orbital decay. Figure 5 shows about 20 hot Jupiters (including WASP-12) for which the theoretical timescale is shorter. In almost all of those cases, though, the planet was discovered more recently than WASP-4, and a decade-long baseline of observations is not yet available. A separate consideration not shown in Figure 5 is that the hot Jupiter host stars have a variety of different structures, from being fully convective to nearly fully radiative, which may lead to widely divergent tidal dissipation timescales.

Zoom In Zoom Out Reset image size

In the simple "constant phase lag" model for tidal interaction (Zahn 1977), the rate of dissipation can be parameterized by a modified¹⁷ quality factor, ${Q}_{\star }^{{\prime} }=3{Q}_{\star }/(2{k}_{\star })$. Here, Q_⋆ is the ratio between the energy stored in the equilibrium deformation of the star and the energy lost to heat per tidal period (e.g., Goldreich & Soter 1966). A larger Q_⋆ implies less efficient tidal dissipation. The dimensionless number k_⋆ is the stellar Love number, which is smaller when the star's density distribution is more centrally concentrated. In this model, once the planet's spin and orbit are synchronized, then the semimajor axis and eccentricity evolve as (Appendix B of Metzger et al. 2012)

$$\begin{eqnarray}&&\displaystyle \frac{1}{{\tau }_{{\rm{e}}}}=\displaystyle \frac{| \dot{e}| }{e}=\displaystyle \frac{63\pi }{2{Q}_{{\rm{p}}}^{{\prime} }}{\left(\displaystyle \frac{{R}_{{\rm{p}}}}{a}\right)}^{5}\left(\displaystyle \frac{{M}_{\star }}{{M}_{{\rm{p}}}}\right)\left(\displaystyle \frac{1}{P}\right)\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&\displaystyle \frac{1}{{\tau }_{{\rm{a}}}}=\displaystyle \frac{| \dot{a}| }{a}=\displaystyle \frac{9\pi }{{Q}_{\star }^{{\prime} }}{\left(\displaystyle \frac{{R}_{\star }}{a}\right)}^{5}\left(\displaystyle \frac{{M}_{{\rm{p}}}}{{M}_{\star }}\right)\left(\displaystyle \frac{1}{P}\right).\end{eqnarray} \tag{ 15 }$$

The orbital period evolves as

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

The modified quality factor of WASP-4 corresponding to the observed value of $\dot{P}$ is

$$\begin{eqnarray}&&{Q}_{\star }^{{\prime} }=(2.9\pm 0.3)\times {10}^{4}.\end{eqnarray} \tag{ 17 }$$

This is about an order of magnitude lower than the value that was inferred for WASP-12b. It is also smaller than most theoreticians would have expected. The ${Q}_{\mathrm{Jup}}^{{\prime} }$ value of Jupiter is estimated to be ≈$1.4\times {10}^{5}$, based on the observed motions of the Galilean moons (Lainey et al. 2009). For stars, studies of the binary eccentricity distribution have been interpreted with tidal models, giving ${Q}_{\star }^{{\prime} }\approx {10}^{5}\mbox{--}{10}^{7}$ (e.g., Meibom & Mathieu 2005; Belczynski et al. 2008; Geller et al. 2013; Milliman et al. 2014). Population studies of hot Jupiter systems have also been undertaken, generally finding ${Q}_{\star }^{{\prime} }\approx {10}^{5}\mbox{--}{10}^{8}$ using different models (Jackson et al. 2009; Hansen 2010; Penev et al. 2012, 2018; Collier Cameron & Jardine 2018). For instance, motivated by the rapid rotation of some hot Jupiter hosts (Pont 2009; Ciceri et al. 2016a; Penev et al. 2016), Penev et al. (2018) modeled the evolution of hot Jupiter systems under the influence of a magnetized wind and a constant phase-lag tide. For WASP-4, their method gave ${Q}_{\star }^{{\prime} }\approx ({1.2}_{-0.5}^{+1.0})\times {10}^{7}$, which would correspond to $\dot{P}\approx -30$ μs per year. This strongly disagrees with the period change that we have observed.

Essick & Weinberg (2016) studied the problem of the orbital decay of hot Jupiters using a theory in which gravity modes are excited at the base of the stellar convective zone, propagate inward through the radiative core, and break near the stellar core, leading to energy dissipation. They predicted the stellar quality factors in hot Jupiter systems to vary in ${Q}_{\star }^{{\prime} }\approx {10}^{5}\mbox{--}{10}^{6}$. From their Equation (26), the prediction for WASP-4 is ${Q}_{\star }^{{\prime} }=7\times {10}^{5}$, which is an order of magnitude larger than implied by the observed period change.

The applicability of the Essick & Weinberg (2016) model depends on the evolutionary state of the star. Weinberg et al. (2017) showed that more rapid dissipation—enough to account for the period change of WASP-12b—could exist in stars that have begun evolving into red giants. The bottom panel of Figure 5 shows a Hertzsprung–Russell diagram of hot Jupiter hosts, including WASP-4. On the y axis is $G=g-\mu$, for g the apparent Gaia-band magnitude and μ the distance modulus reported by Gaia Collaboration et al. (2018). The x axis is the effective temperature from Bonomo et al. (2017), which for WASP-4 agrees within 1σ of that from Table 1. Inspecting the HR diagram, we see WASP-4 shows little evidence of being evolved, in agreement with our analysis from Section 2.3.

To summarize, if the observed period change is caused entirely by tidal orbital decay, then the constant-phase-lag tidal model implies a stellar tidal dissipation rate that is higher than expected by at least an order of magnitude. It might be possible that we are observing at a special time, shortly after the planet's inward migration, or when the planet is near resonance with a stellar oscillation mode. Tidal dissipation rates might also be increased if the star is just turning off the main sequence. Another hypothesis, recently advanced by Millholland & Laughlin (2018) for the case of WASP-12b, is that an exterior planet could be trapping WASP-4b's spin vector in a high-obliquity state, leading to rapid dissipation through planetary obliquity tides.

If instead the observed timing variation is just a small portion of an apsidal precession cycle, then the orbital eccentricity is a few times 10⁻³, and the full precession period is about 27 yr. Ragozzine & Wolf (2009) calculated apsidal precession periods for hot Jupiters, finding them to range between about 10 and 100 yr. They highlighted that for many hot Jupiters, including WASP-4b, the theoretical precession rate is dominated by the non-Keplerian force that is due to the planet's tidal bulge. Precession from general relativity, the planet's rotational bulge, and the star's rotational and tidal bulges contribute at the 10% level at most. Thus, a measurement of the precession rate can be used to determine the planet's Love number. From their Equation (14), the implied Love number for WASP-4b is

$$\begin{eqnarray}&&{k}_{2,{\rm{p}}}={1.59}_{-0.47}^{+0.70}\quad (\mathrm{wide}\,\mathrm{prior}),\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&{k}_{2,{\rm{p}}}={1.20}_{-0.26}^{+0.20}\quad ({k}_{2,{\rm{p}}}\sim { \mathcal U }[0.015,1.5]),\end{eqnarray} \tag{ 19 }$$

where the uncertainties in the semimajor axis and planet radius have been propagated through our Markov chains, and ${ \mathcal U }$ denotes the uniform distribution. For comparison, the Love number of Jupiter is about 0.55 (Wahl et al. 2016; Ni 2018), and a uniform-density sphere has ${k}_{2}=1.5$. The uncertainty in ${k}_{2,{\rm{p}}}$ for WASP-4b is large because the eccentricity, reference time, and $d\omega /{dE}$ have strongly correlated errors, and the four available occultation times only weakly constrain these parameters (Figure 7).

The results under both priors are larger than the Love number of Jupiter, at weak statistical significance. Imposing the physically motivated prior truncates the allowed values of ${k}_{2,{\rm{p}}}$ and drives the precession period toward larger values. Different priors therefore change the predicted evolution of the orbit. While the precession model matches the orbital decay model for the next two decades in transits (Figure 7), in occultations the precession model can deviate from the orbital decay model, particularly if one requires the planetary Love number to be smaller. This suggests that promptly obtaining occultation time measurements for the system could help rule out the apsidal precession model. For the time being, the data are insufficient to make a stronger judgment.

Regardless of the detailed model assumptions, the main physical problem with the apsidal precession hypothesis is explaining why the eccentricity would be as large as ∼10⁻³ despite rapid tidal circularization. For WASP-4, Equation (14) gives ${\tau }_{e}\,=0.29({Q}_{{\rm{p}}}^{{\prime} }/{10}^{5})\,\mathrm{Myr}$. The star is several billion years old, so unless the planet arrived very recently, any initial eccentricity should have been lowered well below ∼10⁻³. The top panel of Figure 5 compares the expected eccentricity damping time of WASP-4b with that of other transiting giant planets. WASP-4b has one of the shortest known eccentricity damping times.

Neighboring companion—One way to maintain a significant eccentricity is through the gravitational perturbations from another planet. Mardling (2007) considered the long-term tidal evolution of hot Jupiters with companions. The companion in their model is coplanar and can have a mass down to an Earth mass; the main requirement is that both the hot Jupiter and the outer companion start on eccentric orbits. They found that although the early phases of the two-planet eccentricity evolution occur quickly, the final phase of the joint eccentricity evolution toward circularity would occur on timescales several orders of magnitude longer than the circularization time of an isolated hot Jupiter (see their Figures 4 and 5).

A separate way a neighboring companion could excite the hot Jupiter's eccentricity is through the Kozai–Lidov mechanism (Kozai 1962; Lidov 1962). In this case, the orbital plane of the outer companion, "c," would need to be inclined relative to that of the hot Jupiter, "b," by at least ${\sin }^{-1}\sqrt{2/5}\approx 39^\circ$. For the Kozai–Lidov mechanism to operate at maximum efficiency, we need (Bailey & Goodman 2019, Equation (20))

$$\begin{eqnarray}&&{M}_{{\rm{c}}}\gt 7.5\,{M}_{\oplus }\times {\left(\displaystyle \frac{{a}_{{\rm{c}}}}{{a}_{{\rm{b}}}}\right)}^{3/2},\end{eqnarray} \tag{ 20 }$$

where ${a}_{{\rm{b}}}$ is the semimajor axis of WASP-4b, ${M}_{{\rm{c}}}$ is the mass of the hypothetical WASP-4c, and ${a}_{{\rm{c}}}$ is WASP-4c's semimajor axis. Owing to our imprecise measurement,  in Equation (20) we have assumed WASP-4b's Love number is ${k}_{2,{\rm{b}}}\approx 0.6$,  similar to Jupiter. For the radial velocity (RV) signal of the companion to remain undetected, it would need to be in the residual ${(O-C)}_{\mathrm{RV}}=15.2\,{\rm{m}}\ {{\rm{s}}}^{-1}$ reported by Triaud et al. (2010). Again following Bailey & Goodman (2019), this implies

$$\begin{eqnarray}\begin{array}{rcl}{M}_{{\rm{c}}} & \lt & {(O-C)}_{\mathrm{RV}}{\left(\displaystyle \frac{{M}_{\star }{a}_{{\rm{b}}}}{G}\right)}^{1/2}{\left(\displaystyle \frac{{a}_{{\rm{c}}}}{{a}_{{\rm{b}}}}\right)}^{1/2}{f}^{-1/2}\\ {M}_{{\rm{c}}} & \lt & 23.7\,{M}_{\oplus }\times {\left(\displaystyle \frac{{a}_{{\rm{c}}}}{{a}_{{\rm{b}}}}\right)}^{1/2}{f}^{-1/2},\end{array}\end{eqnarray} \tag{ 21 }$$

for $f({e}_{{\rm{c}}},{\omega }_{{\rm{c}}},{i}_{{\rm{c}}})\propto {\sin }^{2}{i}_{{\rm{c}}}$, a geometric prefactor that depends on the argument of periastron ${\omega }_{{\rm{c}}}$ and inclination ${i}_{{\rm{c}}}$ of the exterior companion (Bailey & Goodman 2019, Equation (23)). Since WASP-4b is transiting, f can be arbitrarily small, and both of the preceding limits can be satisfied.

Fluctuations in the gravitational potential from convection—An independent mechanism to pump the eccentricity invokes the gravitational fluctuations from stellar convection (Phinney 1992, Section 7). From Equation (7.33) of that work, the mean-squared eccentricity of the orbit is

$$\begin{eqnarray}&&\langle {e}^{2}\rangle =\displaystyle \frac{2\langle {E}_{{\rm{e}}}\rangle }{\mu {n}^{2}{a}^{2}}=6.8\times {10}^{-5}\displaystyle \frac{{({L}^{2}{R}_{\mathrm{conv}}^{2}{M}_{\mathrm{conv}}^{2})}^{1/3}}{\mu {n}^{2}{a}^{2}},\end{eqnarray} \tag{ 22 }$$

where L is the stellar luminosity, ${R}_{\mathrm{conv}}$ and ${M}_{\mathrm{conv}}$ are the width and mass of the convective region, μ is the reduced mass, n is the orbital frequency, and a is the semimajor axis. For the luminosity, reduced mass, and semimajor axis, we used values from Table 1, combined with the Stefan–Boltzmann law and standard definitions. To estimate the width and mass of the convective region, we ran the MESA code for a star with mass and metallicity matched to WASP-4 and the input physics detailed in the MIST isochrones project (Paxton et al. 2011, 2013, 2015; Choi et al. 2016; Dotter 2016). We identified the tachocline boundary using the mixing types specified in the resulting radial profiles and found ${R}_{\mathrm{conv}}\approx 0.33{R}_{\odot }$ and ${M}_{\mathrm{conv}}\approx 9\times {10}^{-4}{M}_{\odot }$. For WASP-4, this implies $\langle {e}^{2}{\rangle }^{1/2}\lesssim {10}^{-5}$. Hence, this mechanism does not seem capable of producing the required eccentricity of ∼10⁻³.

An acceleration of the center of mass of the system toward our line of sight could cause a decrease in the apparent orbital period. The period derivative would be

$$\begin{eqnarray}&&\dot{P}=\displaystyle \frac{{\dot{v}}_{{\rm{r}}}P}{c},\end{eqnarray} \tag{ 23 }$$

where ${\dot{v}}_{{\rm{r}}}$ is the time derivative of the radial velocity.

Knutson et al. (2014), combining radial-velocity data from their own program with those of Wilson et al. (2008), Pont et al. (2011), and Husnoo et al. (2012), found evidence for a long-term trend in WASP-4, with low statistical significance:

$$\begin{eqnarray}&&{\dot{v}}_{{\rm{r}}}=-{0.0099}_{-0.0054}^{+0.0052}\,{\rm{m}}\,{{\rm{s}}}^{-1}\,{\mathrm{day}}^{-1}.\end{eqnarray} \tag{ 24 }$$

This acceleration translates to an expected $\dot{P}=-(4.4\pm 2.4)\,\times {10}^{-11}$. The period decrease from the observed RV trend is an order of magnitude smaller than the observed $\dot{P}$ from transit times, $(4.0\pm 0.4)\times {10}^{-10}$ (Table 4).

Nonetheless, it is intriguing that the host star shows a weakly significant acceleration, and with the correct sign needed to explain the transit timing variations. Given the potential importance of the radial velocities for interpreting this system, we performed an independent analysis, as follows.

First, we collected the usable RV measurements from CORALIE, HARPS, and HIRES. We included the CORALIE measurements from Wilson et al. (2008) and Triaud et al. (2010), using the homogeneous radial velocities calculated by the latter authors. We included the HARPS values reported by Pont et al. (2011), which are identical to those from Husnoo et al. (2012). We omitted the HARPS data points taken over three nights by Triaud et al. (2010) for Rossiter–McLaughlin observations because they were calculated using a different pipeline than the longer-baseline Pont et al. HARPS measurements, and necessary inclusion of an extra offset term would nullify their statistical value. Finally, we included the five HIRES measurements taken over many years by Knutson et al.(2014).

We then fitted a single Keplerian orbit, plus instrument offsets, jitters, and a long-term trend (Fulton et al. 2018, radvel). We set Gaussian priors on the period and time of inferior conjunction using the values from Table 4, and we fixed the eccentricity to zero, consistent with results from Beerer et al. (2011), Knutson et al. (2014), and Bonomo et al. (2017). The remaining free parameters were the velocity semiamplitude, the instrument zero points, the instrument jitters (an additive white noise term for each instrument), and linear ($\dot{{v}_{{\rm{r}}}}$) and optionally second-order ($\ddot{{v}_{{\rm{r}}}}$) acceleration terms.

We found that the best-fitting model with both linear and quadratic radial velocity terms was marginally preferred (by ${\rm{\Delta }}\mathrm{BIC}=5.8$) over the best-fitting model with only a linear term. Regardless, for consistency with Knutson et al. (2014), who fixed the quadratic component of the long-term trend to zero, in Figure 6 we show best-fitting models for the linear-trend case. The best-fit value for the line-of-sight acceleration,

$$\begin{eqnarray}&&{\dot{v}}_{{\rm{r}}}=-{0.0077}_{-0.0047}^{+0.0052}\,{\rm{m}}\,{{\rm{s}}}^{-1}\,{\mathrm{day}}^{-1},\end{eqnarray} \tag{ 25 }$$

is within 1σ of the value reported by Knutson et al. (2014). (Note that the CORALIE data used in our and their analyses differ, as we included additional measurements reported by Triaud et al. 2010.) The implied period derivative is still therefore about an order of magnitude smaller than our observed $\dot{P}$ from transit timing.

Zoom In Zoom Out Reset image size

To summarize, only about one-tenth of the observed period decrease can be explained through a constant acceleration of the WASP-4 system's center of mass. However, given the limited amount and uneven time coverage of the existing radial-velocity data (Figure 6), it remains possible that the center of mass has a more complex motion, perhaps due to a companion on an eccentric orbit (e.g., WASP-53 or WASP-81, Triaud et al. 2017). It would be useful to gather more radial-velocity data to confirm or refute this possibility.

A separate candidate explanation for the timing deviations is the Applegate (1992) effect. Some eclipsing binaries exhibit period modulations with amplitudes of ≲0.05 day over timescales of decades (e.g., Söderhjelm 1980; Hall 1989). The Applegate mechanism explains these modulations by positing that the internal structure of a magnetically active star changes shape via cyclic exchange of angular momentum between the inner and outer zones of the star. This model could also apply to a hot Jupiter orbiting a star with a convective zone. The changing gravitational quadrupole of the star would cause the orbit of the planet to precess on the timescale of the stellar activity cycle. An essential difference between this process and apsidal precession is that Applegate timing variations need not be strictly periodic (e.g., Söderhjelm 1980, Figure 12). This mechanism would also produce transit and occultation timing deviations of the same sign, while for apsidal precession they would have opposite signs. For WASP-4, Watson & Marsh (2010) estimated that the Applegate effect could produce timing deviations of up to 15 s, depending on the modulation period of the stellar dynamo. If this analysis is accurate, then the Applegate mechanism cannot explain the majority of our observed 82 s variation.

There are two other small effects worth noting. The first is the Shklovskii (1970) effect, due to the star's proper motion, which leads to an apparent period change of $P{\mu }^{2}d/c$, which is only 6 × 10⁻¹³ for the case of WASP-4. The second effect, described by Rafikov (2009), comes from the star's on-sky motion altering our viewing angle and leads to an observed apsidal precession. The corresponding period change is $\dot{P}\sim {(P\mu )}^{2}/2\pi$, which is on the order of 10⁻²¹ for WASP-4, too small to be of any consequence.