Characterizing the WASP-4 system with TESS and radial velocity data: Constraints on the cause of the hot Jupiter's changing orbit and evidence of an outer planet

Orbital dynamics provide valuable insights into the evolution and diversity of exoplanetary systems. Currently, only one hot Jupiter, WASP-12b, is confirmed to have a decaying orbit. Another, WASP-4b, exhibits hints of a changing orbital period that could be caused by orbital decay, apsidal precession, or the acceleration of the system towards the Earth. We have analyzed all data sectors from NASA's Transiting Exoplanet Survey Satellite together with all radial velocity (RV) and transit data in the literature to characterize WASP-4b's orbit. Our analysis shows that the full RV data set is consistent with no acceleration towards the Earth. Instead, we find evidence of a possible additional planet in the WASP-4 system, with an orbital period of ~7000 days and $M_{c}sin(i)$ of $5.47^{+0.44}_{-0.43} M_{Jup}$. Additionally, we find that the transit timing variations of all of the WASP-4b transits cannot be explained by the second planet but can be explained with either a decaying orbit or apsidal precession, with a slight preference for orbital decay. Assuming the decay model is correct, we find an updated period of 1.338231587$\pm$0.000000022 days, a decay rate of -7.33$\pm$0.71 msec/year, and an orbital decay timescale of 15.77$\pm$1.57 Myr. If the observed decay results from tidal dissipation, we derive a modified tidal quality factor of $Q^{'}_{*}$ = 5.1$\pm$0.9$\times10^4$, which is an order of magnitude lower than values derived for other hot Jupiter systems. However, more observations are needed to determine conclusively the cause of WASP-4b's changing orbit and confirm the existence of an outer companion.


INTRODUCTION
One of the most powerful tools used to study exoplanets is observing them while they transit across the disk of their stars. The transit method can be used to search for temporal variations in the planetary orbital parameters and allows for a direct measurement of the planetary radius and thus its atmosphere (e.g. Charbonneau et al. 2000Charbonneau et al. , 2007. Specifically, light curves of transiting planets can be used to search for transit timing variations (TTVs; Agol et al. 2005;Agol & Fabrycky 2018), transit duration variations (Agol & Fabrycky 2018), and impact parameter variations (Herman et al. 2018;Szabó et al. 2020). The presence of TTVs may indicate addi-tional bodies in the system, a decaying planetary orbit, precession in the orbit, or a variety of other effects (e.g., Applegate 1992, Dobrovolskis & Borucki 1996, Miralda-Escudé 2002Schneider 2003;Agol et al. 2005;Mazeh et al. 2013;Agol & Fabrycky 2018).
It is theorized that orbital decay might occur on short-period massive planets orbiting stars with surface convective zones due to exchange of energy with their host stars through tidal interactions (e.g., Rasio et al. 1996;Lin et al. 1996;Chambers 2009;Lai 2012;Penev et al. 2014;Barker 2020). Measurements of orbital decay would expand our understanding of the hot Jupiter population and its evolution (e.g. Jackson et al. 2008;Hamer & Schlaufman 2019). Even though such decay may occur over millions of years, it is possible to search for small changes in the orbital period of hot Jupiter sys-tems since many of these planets have been monitored for decades.
Currently, WASP-12b is one of the few hot Jupiters confirmed to have a varying period. It is an ultra-hot planet around a G0 star with an orbital period of 1.09 days (Hebb et al. 2009). WASP-12b is believed to have an escaping atmosphere (e.g., Lai et al. 2010;Bisikalo et al. 2013;Turner et al. 2016a) as suggested by Hubble Space Telescope near-ultraviolet observations (Fossati et al. 2010;Haswell et al. 2012;Nichols et al. 2015). Maciejewski et al. (2016) were the first to detect its decreasing orbital period, and subsequent studies have confirmed the period change (Patra et al. 2017;Maciejewski et al. 2018;Bailey & Goodman 2019;Baluev et al. 2019) and established orbital decay as its cause (Yee et al. 2020). The decaying orbit of WASP-12b was confirmed also using transit and occultation observations with NASA's Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2015) (Turner et al. 2021, see also Owens et al. 2021). The decay rate of WASP-12b was found to be 32.53±1.62 msec yr −1 corresponding to an orbital decay timescale of τ = P/|Ṗ | = 2.90 ± 0.14 Myr (Turner et al. 2021), shorter than the estimated mass-loss timescale of 300 Myr (Lai et al. 2010;Jackson et al. 2017). Assuming the observed decay results from tidal dissipation, Turner et al. (2021) derived a modified tidal quality factor, a dimensionless quantity that describes the efficiency of tidal dissipation, of Q = 1.39±0.15 ×10 5 , which falls at the lower end of values derived for binary star systems (Meibom et al. 2015) and hot Jupiters (Jackson et al. 2008;Husnoo et al. 2012;Barker 2020).
Besides WASP-12b, WASP-4b is a well-studied system with hints of a changing period. WASP-4b is a typical hot Jupiter with a short orbital period of 1.34 days and orbits a G7V star . The planet's atmosphere (e.g. Cáceres et al. 2011;Ranjan et al. 2014;Bixel et al. 2019) and orbital parameters (e.g. Southworth 2009;Winn et al. 2009a;Hoyer et al. 2013;Baluev et al. 2020) have been studied extensively since its discovery in 2008. Bouma et al. (2019) first reported an orbital period variation of WASP-4b using TESS and ground-based observations. Southworth et al. (2019) confirmed that the period was decaying with additional ground-based observations and found a decay rate of 9.2 msec yr −1 . They found that orbital decay and apsidal precession could explain the TTVs after ruling out instrumental issues, stellar activity, the Applegate mechanism, and light-time effect. Bouma et al. (2020) obtained additional radial-velocity (RV) data on WASP-4b using HIRES on Keck and found that the observed orbital period variation could be explained by the sys-tem accelerating toward the Sun at a rate of -0.0422 m s −1 day −1 . Recently, Baluev et al. (2020) analyzed a comprehensive set of 129 transits and additional RV data (presented in Baluev et al. 2019(presented in Baluev et al. ) from 2007(presented in Baluev et al. -2014 (mostly in years not covered by the RV data presented in Bouma et al. 2020). They also confirmed a period change in WASP-4b's orbit but do not confirm the RV acceleration found by Bouma et al. (2020). However, Baluev et al. (2020) did not include the new RV HIRES data from Bouma et al. (2020) in their analysis. Therefore, the cause of the period variation in WASP-4b is still an open question.
Motivated by the possible changing period of WASP-4b, we analyze all three sectors (Figure 1) from TESS and combine the results with all transit, occultation, and RV measurements from the literature to verify its changing period and derive updated planetary properties. TESS is well suited for our study because it provides high-precision time-series data, ideal for searching for TTVs (e.g. Hadden et al. 2019;Pearson 2019;von Essen et al. 2020;Ridden-Harper et al. 2020;Turner et al. 2021). Altogether, the transit data span 13 years from 2007-2020 and the RV data span 12 years from 2007-2019. Using the combined data set we hope to shed light on the cause of WASP-4b's period variation.

OBSERVATIONS AND DATA REDUCTION
TESS observed WASP-4b in Sector 2 (2018-Aug-22 to 2018-Sep-20), Sector 28 (2020-Jul-30 to 2020-Aug-26), and Sector 29 (2020-Aug-26 to 2020. The TESS observations were processed by the Science Processing Operations Center (SPOC) pipeline 1 (Jenkins et al. 2016). The SPOC pipeline produces light curves ideal for characterizing transiting planets since they are corrected for systematics. SPOC produces Presearch Data Conditioning (PDC) and Data Validation Timeseries (DVT) light curves. The PDC light curves are corrected for instrumental systematics (pointing or focus related), discontinuities resulting from radiation events in the CCD detectors, outliers, and flux contamination. The DVT light curves are created by using a running median filter on the PDC light curves to remove any long-period systematics. We use the DVT light curves ( Figure 1) for our analysis because they have less scatter in their out-of-transit (OoT) baseline. As shown in Ridden-Harper et al. (2020) for XO-6b TESS data, the DVT and PDC light curves produce similar results on the timing of the transits. For the Sector 2 data, the light curves produced from the SPOC pipeline have a known issue 2 that overestimates their uncertainties. Therefore, we estimated the uncertainties using the scatter in the OoT baseline as we did in our previous study of WASP-12b (Turner et al. 2021) and as recommended (Barclay, T., private communication). The Sector 28 and 29 data are unaffected by this problem. Therefore, we used the uncertainties provided by the SPOC pipeline.

Transit Modeling
All the TESS transits of WASP-4b were modeled with the EXOplanet MOdeling Package (EXOMOP; Pearson et al. 2014;Turner et al. 2016bTurner et al. , 2017 3 to find a bestfit. EXOMOP creates a model transit using the analytic equations of Mandel & Agol (2002)  Each TESS transit ( Figure 1) was modeled with EXOMOP independently. We used 20 6 links and 20 chains for the DE-MCMC model and use the Gelman-Rubin statistic (Gelman & Rubin 1992) to ensure chain convergence (Ford 2006). The mid-transit time (T c ), scaled semi-major axis (a/R * ), planet-to-star radius (R p /R * ), and inclination (i) are set as free parameters for every transit. The linear and quadratic limb darkening coefficients and period (P orb ) are fixed during the analysis. The linear and quadratic limb darkening coefficients are taken from Claret (2017) and are set to 0.382 and 0.210, respectively.
The parameters derived for every TESS transit event can be found in Tables A1-A3 in Appendix A. The modeled light curves for each individual transit can be found in Figures A1-A7 in Appendix A. All parameters for each transit event are consistent within 2σ of every other transit. The Sector 2 data for WASP-4b was analyzed in Bouma et al. (2019) and our timing analysis for each individual transit is consistent within 1σ of their findings (see Figure B8 in Appendix B).
The light curve of WASP-4b was phase-folded at each derived mid-transit time and modeled with EXOMOP to find the final fitted parameters. The phase-folded light curve and model fit can be found in Figure 2. We use the light curve model results combined with literature values to calculate the planetary mass (M p ; Winn 2010), radius (R p ), density (ρ p ), equilibrium temperature (T eq ), surface gravity (log g p ; Southworth et al. 2007), or-  Table 4 from the best-fit 2-planet RV model (Model #3). The period of WASP-4b was taken from the orbital decay model in Table 6. All the physical parameters for WASP-4c were taken from the best-fit 2-planet RV model (Model #3).

Occultation
We created an occultation light curve by phase-folding all the data about the secondary eclipse using the first TESS transit as the reference transit time. As shown in Figure 3, we do not see an occultation of WASP-4b in the TESS data. We only show the PDC light curve but this is also the case for the DVT light curve. We find a 3-σ upper limit on the occultation depth (δ occ ) to be 1.34 ×10 −5 . The geometric albedo (A g,occ ) of WASP-4b 4 The Safronov number is a qualitative measure of the potential of a planet to capture or gravitationally scatter other objects in close by orbits can be calculated assuming no thermal contribution (1) We find a 3-σ upper limit on A g,occ of 0.017 using Equation (1). Our upper limit on A g,occ is consistent with the overall trend that hot Jupiters are very dark (e.g. Kipping & Bakos 2011;Močnik et al. 2018;Kane et al. 2020).

Modeling Radial Velocities
We combined all available radial velocity (RV) data of WASP-4 in the literature in our analysis. The complete RV data set is given in Table 2. The set includes CORALIE and HARPS data from Baluev et al. (2019), HARPS data from Pont et al. (2011), and HIRES data from Bouma et al. (2020). We did not include the RV data from Triaud et al. (2010) as these data were taken for Rossiter-McLaughlin measurements and were reduced with a nonstandard pipeline. Note-This table is available in its entirety in machine-readable form. Note-The number of data points (N data ) for all the models was 74.
priors for the orbital period and time of inferior conjunction that were set to the values derived by Bouma et al. (2019). We set the eccentricity of the orbit to zero as indicated by previous upper limit studies (Beerer et al. 2011;Knutson et al. 2014;Bonomo et al. 2017. As done in previous studies, we first modeled the data with only one planet. The free parameters in the model were the orbital velocity semi-amplitude (K b ), the instrument zero-points, white noise instrument jitter for each instrument (σ, added in quadrature to its uncertainties), and the linear acceleration (γ). We ran models with and withoutγ to check if an acceleration is needed to fit the data (Models #1 and #2). The results of the 1-planet models of WASP-4b can be found in Table 4 and Figure 4. We use the Bayesian Information Criterion (BIC) to assess the preferred model. The BIC is defined as where k is the number of free parameters in the model fit and N pts is the number of data points. The power of the BIC is the penalty for a higher number of fitted model parameters, making it a robust way to compare different best-fit models. The preferred model is the one that produces the lowest BIC value. We find a BIC value of 675.08 and 682.76 for the model without (Model #1) and with fittingγ (Model #2), respectively. For two generic models i and j, we can relate the difference in the BIC values between models, ∆(BIC j,i ) = BIC j -BIC i , and the Bayes factor B i,j , the ratio of the likelihood between models i and j, assuming a Gaussian distribution for the posteriors (e.g. Faulkenberry 2018): Therefore, Model #1 without fitting forγ is the preferred model with a ∆(BIC 1,2 ) = BIC 1 -BIC 2 = -7.68 and Bayes factor, B 1,2 , of 46.5. When fitting forγ, we find a linear acceleration term that is positive but is consistent with zero within the uncertainties (γ = 0.0001 +0.0034 −0.0036 m s −1 day −1 ). Our results are in conflict with the findings by Bouma et al. (2020) that find an acceleration along our line of sight at a rate oḟ γ = −0.0422 +0.0028 −0.0027 m s −1 day −1 . We can reproduce the results of Bouma et al. (2020) by modeling only their data (See Figure C13 and Table C5 in Appendix C). Based on this test, we conclude that the acceleration found by Bouma et al. (2020) was caused by modeling only part of the full RV data set. Therefore, we conclude that the changing orbital period detected using the transit data is not caused by the WASP-4 system accelerating towards Earth.
Next, we fit the RV data with a 2-planet model because the residuals of the one-planet fit showed some sinusoidal structure (panel b in Figures 4I-II). This sinusoidal trend is not caused by stellar activity as the S-index time series from the HIRES data shows no signs of secular or sinusoidal trends (Bouma et al. 2020). We performed several different models where we fit for K b , σ,γ, the orbital velocity semi-amplitude of the 2nd body (K c ). and an eccentricity of the 2nd planet (e c ). The results of Models #3 and #4 are summarized in Table 3. We find that Model #3 with K b and K c set as free parameters and e b , e c , andγ fixed to zero finds the best fit with a BIC of 628.34. The derived orbital parameters of this model can be found in Table 4 and the two-planet RV fit can be found in Figures 4III-IV. The two-planet model (Model #3) is highly preferred over the one-planet model (Model #1) with a ∆(BIC 3,1 ) = -46.74 and a Bayes factor, B 3,1 , of 1.41×10 10 . For the second body, we find a period (P c ) of 7001.0 +6.0 −6.6 days, a semi-major axis of 6.82 +0.23 −0.25 AU, and a M c sin(i) of 5.47 +0.44 −0.43 M Jup (Table 1, Table 4). The companion is expected to be much fainter than the planet host star because Wilson et al. (2008) found no evidence for changing spectral line bisectors in their spectroscopic observation. Becker et al. (2017) found that distant exterior companions to hot Jupiters around cool stars (T star < 6200 K) are typically coplanar within 20-30 degrees. Therefore, we find that the mass of WASP-4c is between 5.47-6.50 M Jup assuming that its inclination is within 30 degrees of WASP-4b's inclination. However, more RV measurements are needed to verify the existence of this second planet around WASP-4.

Transit Timing Variations
For the timing analysis, we combined the TESS transit data with all the prior transit and occultation times. All the transit and occultation times used in this analysis can be found in Table 5. In the table we give the  Table 4 are the median values of the posterior distributions. The thin blue line is the best fit model. We add in quadrature the RV jitter terms listed in Table 4 with the measurement uncertainties for all RVs. b panels) Residuals to the best fit model. The error bars of the residuals reflect both the measurement error and the jitter from the MCMC fit. The larger error bars for the 1 planet fits reflect the larger amounts of jitter needed to fit the data with only 1 planet. c panels) RVs phase-folded to the ephemeris of WASP-4b. The phase-folded model for WASP-4b is shown as the blue line. The Keplerian orbital models for all other planets (if modeled) have been subtracted. d panels) RVs phase-folded to the ephemeris of the second planet WASP-4c. The Keplerian orbital model for the other planet (panels c, WASP-4b) has been subtracted. Red circles are the velocities binned in 0.08 units of orbital phase.     Note-This table is available in its entirety in machine-readable form.  original reference in which the data were reported and the reference for the timing if different from the original source. We combine transit data as tabulated by Bouma et al. (2020) and additional transits reported by amateur observers from Baluev et al. (2020). This table is available in its entirety in machine-readable form online. Similar to what was done in Yee et al. (2020), Patra et al. (2017), andTurner et al. (2021), we fit the timing data to three different models. The first model is the standard constant period formalization:

References-
where T 0 is the reference transit time, P orb is the orbital period, E is the transit epoch, and T tra (E) is the calculated transit time at epoch E.
The second model assumes that the orbital period is changing uniformly over time: where dP orb /dE is the decay rate. The third model assumes the planet is precessing uniformly (Giménez & Bastero 1995): where e is a nonzero eccentricity, ω is the argument of pericenter, P s is the sidereal period and P a is the anomalistic period. For all three models, we found the best-fitting model parameters using a DE-MCMC analysis. We used 20 chains and 20 6 links in the model and again we ensure chain convergence using the Gelman-Rubin statistic.
The results of timing model fits can be found in Table 6. Figure 5 shows the transit and occultation timing data fit with the orbital decay and apsidal precession models. In this figure, the best-fit constant-period model has been subtracted from the timing data. The orbital decay model fits the transit and occultation data Orbital Decay Apsidal Precession Figure 5. WASP-4b transit (panels a, b, and c) and occultation (panel d) timing variations after subtracting the data with a constant-period model. The filled black triangles are the data points from the literature and the square orange points are from the TESS data in this paper. All the transit and occultation times can be found in Table 5. The orbital decay and apsidal precession models are shown as the blue and red lines, respectively.
slightly better than the precession model (Table 6, Figure 5).
Our finding that the constant-period model does not fit the data well is consistent with previous studies (Bouma et al. 2019;Southworth et al. 2019;Baluev et al. 2020). The orbital decay and apsidal precession models fit the data with a minimum chi-squared (χ 2 min ) of 276.35 and 270.42, respectively. We find that the orbital decay model is the preferred model with a ∆(BIC) = -5.93 and a Bayes factor of 9.3. Therefore, based on our analysis of the observed timing residuals, the orbital decay model is only slightly preferred over apsidal precession.
Due to the RV measurements showing evidence of a possible second planet, we modeled the two-planet system to see if they could reproduce the TTVs. For this analysis, we used the publicly available TTV analysis package, OCFit 5 (Gajdoš & Parimucha 2019). Specifi-5 https://github.com/pavolgaj/OCFit cally, we used the AgolExPlanet function which is an implementation of Equation 25 in Agol et al. (2005). The priors in OCFit were set to the values given for Model #4 in Table 4 for both planets where the outer planet has an eccentricity of 0.094. We were not able to fit the TTVs well with an outer planet consistent with the RV constraints. We also produced several forward models with OCFit that show that the expected TTV signal from the outer body is less than 2 seconds (dependent on the real mass of the body) over the full observational period (Figure 6). We did not use the preferred two-planet model (Model #3) because this model did not produce a detectable signal. The two objects are assumed to be co-planar but relaxing this condition will only decrease the TTV signal.
We can also analytically calculate the expected TTV signal on WASP-4b (δt b ) using the following equation from Agol et al. (2005) assuming non-resonant planets with large period-ratios on circular orbits:  Table 4 where the outer planet has an eccentricity of 0.094 (Model #4). The two objects are assumed to be co-planar.
where M c and P c are the mass and orbital period of the outer companion, M star is the mass of WASP-4, and P b is the period of WASP-4b. The assumptions of equation (12) are all satisfied within the constraints of the best-fit RV model parameters found by RadVel (Table 4). For a M c between 5.47 M Jup and 2 M sun , we find a δt b using equation (12) to be between 1.9×10 −10 secs to 7.4×10 −8 secs. Hence, the expected TTV signal is many orders magnitude below the observed TTVs regardless of the mass of the outer companion. Our results are expected as resonant perturbations between close planets is the main cause of large (∼>mins) TTVs (e.g. Agol et al. 2005;Steffen et al. 2012;Nesvorný et al. 2013;Dawson et al. 2019). Therefore, we conclude that the observed TTVs are caused by orbital decay or apsidal precession and not gravitational perturbations from the outer body.

DISCUSSION
From our analysis, we derived updated planetary parameters and find that the orbital decay model is slightly preferred to explain the observed TTVs. We conclusively rule out linear acceleration as the cause of the period change. We also show that TTVs can not be caused by the second body orbiting WASP-4 regardless of its mass. However, apsidal precession is not ruled out. Further transit and occultation measurements of WASP-4b are needed to disentangle the cause of the variation. The orbital decay and apsidal precession models exhibit very different timing variations in the mid-2020s ( Figure  7). Therefore, we predict that it will be possible to conclusively determine the cause of WASP-4b's changing orbit by then.  Figure 7. WASP-4b transit (panels a) and occultation (panel b) timing variations predictions using the best-fit models to the timing data presented in this paper. The lines are 100 random draws from the posteriors of the orbital decay (blue) and apsidal precession (red) models. On both panels, we show the average error of the previous observations.
Assuming the TTVs can be explained entirely by orbital decay, we derive an updated period of 1.338231587±0.000000022 days and a decay rate of -7.33±0.71 msec year −1 . Our results indicate an orbital decay timescale of τ = P/|Ṗ | = 15.77 ± 1.57 Myr, slightly longer than the value derived by Bouma et al. (2019) of 9.2 Myr. Assuming that the planet mass is constant, the rate of change of WASP-4b's orbital period,Ṗ , can be related to its host star's modified tidal quality factor by the constant-phase-lag model of Goldreich & Soter (1966), defined aṡ where M p is the mass of the planet, M is the mass of the host star, R is the radius of the host star and a is the semi-major axis of the planet. Using our derived value ofṖ and planetary parameters from Table 1, we find a modified tidal quality factor of Q = 5.1±0.9 ×10 4 . This value is higher than the value found by The cause of this small discrepancy (all Q values are within 2σ) is due to our including more transit data in our analysis than previous studies. It is currently not clear how to account theoretically for all the observed low values of Q . Our value is an order of magnitude lower than the observed values for binary star systems (10 5 -10 7 ; Meibom et al. 2015) and hot Jupiters (10 5 − 10 6.5 ; e.g., Jackson et al. 2008;Husnoo et al. 2012;Barker 2020) (2019) who found that hot Jupiter host stars tend to be young requiring Q 10 7 . Some possible causes of a large tidal dissipation rate might be that WASP-4 is turning off the main sequence (as suspected for WASP-12; Weinberg et al. 2017) or that an exterior planet could be trapping WASP-4b's spin vector in a high-obliquity state (also predicted for WASP-12b; Millholland & Laughlin 2018). The latter theory is intriguing as we now have evidence for an additional body in the system (Figure 4). More investigation is needed to understand the low value of Q .
If confirmed, WASP-4b would be the second exoplanet after WASP-12b (Yee et al. 2020;Turner et al. 2021) to show evidence of tidal orbital decay. Future observations of WASP-4b are needed to verify this possibility. Other planets such as WASP-103b, KELT-16b, and WASP-18b are also predicted to exhibit large rates of tidal decay (Patra et al. 2020). Additionally, the planets HATS-2b and WASP-64b are also ideal candidates to search for orbital decay because they are in systems similar to WASP-4 (Southworth et al. 2019). To understand the formation and evolution of the hot Jupiter population, timing observations of additional systems are needed.
The discovery of WASP-4c makes the WASP-4 system unique, as it would be the most widely separated companion of a transiting hot-Jupiter known to date. However, our discovery of WASP-4c is not surprising because Knutson et al. (2014) estimated that 27% ± 6% of hot Jupiters have a planetary companion at a distance of 1-10 AU and a mass of 1-13 M Jup . Similar planets may be common in other systems but could only be detected with long time-baseline RV data sets. Therefore, we expect that the unique status of the WASP-4 system is a result of observational bias rather than an intrinsic rarity of such systems.
Future observations of the WASP-4 system can help us put the system in context with the overall hot Jupiter population and shed light on the possible formation scenarios of the system. Specifically, stronger constraints on the obliquity, mutual inclinations, and full orbital parameters of WASP-4c will help us understand planetary migration in this system.

CONCLUSIONS
We analyzed all available sectors of TESS data of WASP-4b to investigate its possible changing orbit. Our TESS transit timing investigations confirm that the planet's orbit is changing ( Figure 5). We conclude that the acceleration of the WASP-4 system towards Earth is not the cause of the period variation after analyzing all available RV data (Figure 4; Table 4). From the RV analysis, we also find evidence of a possible second planet orbiting WASP-4 with a period (P c ) of 7001.0 +6.0 −6.6 days, semi-major axis of 6.82 +0.23 −0.25 AU, and a M c sin(i) of 5.47 +0.44 −0.43 M Jup (Figure 4, Table 4). WASP-4c is the most widely separated companion of a transiting hot-Jupiter discovered to date. This outer planet is not the cause of the observed TTVs ( Figure 6). Our timing analysis slightly favors the orbital decay scenario over apsidal precession as the cause of the TTVs (Table 6, Figure  5). However, apsidal precession cannot be ruled out. We find an updated period of 1.338231587±0.000000022 days, a decay rate of -7.33±0.71 msec year −1 , and an orbital decay lifetime of 15.77±1.57 Myr assuming the system is undergoing orbital decay. The planetary physical parameters are also updated with greater precision than previous studies. More transit, occultation, and RV data are needed over the next few years to determine conclusively the cause of WASP-4b's changing orbit and help place the system in context with the overall hot Jupiter population. ACKNOWLEDGMENTS Support for this work was provided by NASA through the NASA Hubble Fellowship grant #HST-HF2-51495.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555.
We are grateful to Benjamin (BJ) Fulton for his help with RadVel. We thank Dong Lai and Maryame El Moutamid for useful discussions. We were inspired to pursue this project after attending a TESS hackathon hosted by the Carl Sagan Institute.
We also thank the anonymous referee for the useful comments.
This paper includes data collected by the TESS mission, which are publicly available from the Mikulski Archive for Space Telescopes (MAST    Figure A2. Individual TESS transit events (8-16) from Sector 20 of WASP-4b. Other comments are the same as Figure A1.  Figure A3. Individual TESS transit events (17-18) from Sector 2 of WASP-4b. Other comments are the same as Figure A1.  Figure A4. Individual TESS transit events (1-8) from Sector 28 of WASP-4b. Other comments are the same as Figure A1.  Figure A5. Individual TESS transit events (9-15) from Sector 28 of WASP-4b. Other comments are the same as Figure A1.  Figure A6. Individual TESS transit events (1-8) from Sector 29 of WASP-4b. Other comments are the same as Figure A1.  Figure A7. Individual TESS transit events (9-14) from Sector 29 of WASP-4b. Other comments are the same as Figure A1.

B. DIFFERENCE IN MID-TRANSIT TIMES FOR THE SECTOR 2 TESS DATA
The comparison in the derived mid-transit times between our analysis and Bouma et al. (2020) for the Sector 2 TESS data is found in Figure B8. Our results are consistent within 1σ.  Figure B8. Difference in mid-transit times between our analysis and Bouma et al. (2020) for the Sector 2 TESS data.

C. RADIAL-VELOCITY MODELS
We performed an RV analysis on all the RV data in the literature with RadVel. Table 3 compares all the one-planet and two-planet RV models carried out in our analysis. The priors used in RadVel can be found in Table C4. The best fit of the data is the two-planet model (Model # 3 in Table 3) with a BIC of 628.34. In this model, the eccentricities of both bodies (e b and e c ) are fixed to zero and there is no linear acceleration. The planetary parameters derived of all models can be found in Table 4. The posterior distributions for all free parameters of Model # 1 -4 can be found in Figures  C11-C12. The best-fit Keplarian orbital models compared to the RV data can be found in Figure 4.
We also performed an RV analysis only on the data presented in Bouma et al. (2020). We used priors listed in Table C4. The results of that analysis can be found in Table C5 and Figure C13. We find aγ = −0.0400 +0.0037     Figure C10. Posterior distributions for all free parameters for the one-planet RV analysis without (Model # 2) fitting for the linear acceleration. The results of this model can be found in Table 4 and Figure 4.  Figure C11. Posterior distributions for all free parameters for the best-fit two-planet RV analysis (Model # 3). The results of this model can be found in Table 4 and Figure 4.  Figure C12. Posterior distributions for all free parameters for the two-planet RV analysis (Model # 4) fitting for ec. The results of this model can be found in Table 4 and Figure 4.   Table C5 are the median values of the posterior distributions. The thin blue line is the best fit 1-planet model. We add in quadrature the RV jitter terms listed in Table C5 with the measurement uncertainties for all RVs. b) Residuals to the best fit 1-planet model. c) RVs phase-folded to the ephemeris of WASP-4b. The small point colors and symbols are the same as in panel a. Red circles are the same velocities binned in 0.08 units of orbital phase.
The phase-folded model for WASP-4b is shown as the blue line.