A Lack of Variability Between Repeated Spitzer Phase Curves of WASP-43b

Though the global atmospheres of hot Jupiters have been extensively studied using phase curve observations, the level of time variability in these data is not well constrained. To investigate possible time variability in a planetary phase curve, we observed two full-orbit phase curves of the hot Jupiter WASP-43b at 4.5 microns using the Spitzer Space Telescope, and reanalyzed a previous 4.5 micron phase curve from Stevenson et al. (2017). We find no significant time variability between these three phase curves, which span timescales of weeks to years. The three observations are best fit by a single phase curve with an eclipse depth of 3907 +- 85 ppm, a dayside-integrated brightness temperature of 1479 +- 13 K, a nightside-integrated brightness temperature of 755 +- 46 K, and an eastward-shifted peak of 10.4 +- 1.8 degrees. To model our observations, we performed 3D general circulation model simulations of WASP-43b with simple cloud models of various vertical extents. In comparing these simulations to our observations, we find that WASP-43b likely has a cloudy nightside that transitions to a relatively cloud-free dayside. We estimate that any change in WASP-43bs vertical cloud thickness of more than three pressure scale heights is inconsistent with our observed upper limit on variation. These observations, therefore, indicate that WASP-43bs clouds are stable in their vertical and spatial extent over timescales up to several years. These results strongly suggest that atmospheric properties derived from previous, single Spitzer phase curve observations of hot Jupiters likely show us the equilibrium properties of these atmospheres.


INTRODUCTION
Weather and climate are intrinsically linked to interpretations of planetary properties deduced from remote observations. The history of observations within our own Solar System tells cautionary tales for neglecting weather phenomena. For example, early observers of Mars noticed seasonal variations in dark spots on Mars' surface, and interpreted them as the coming and going of vast patches of vegetation (see, e.g., Sinton 1957Sinton , 1959. This variation was later revealed to just be the seasonal obscuration and uncovering of the rocky, inanimate surface by dust storms. Exoplanetary science is entering a revolutionary era for high-precision spectroscopy of exoplanet atmospheres, with the potential for similar difficulties in interpreting newfound phenomena.
Our understanding of an exoplanet's atmosphere is unavoidably linked to the effects of weather and climate in that atmosphere. This is because variability in atmospheric mixing, cloud formation, and circulation can drive the planet's observable properties to change in time. As a result, revealing the true nature of the atmosphere may require repeated observations. In an atmosphere that varies in time, properties derived from only a single observation provide only a snapshot state arXiv:2212.03240v3 [astro-ph.EP] 15 Feb 2023 that does not necessarily represent the general nature of a planet.
Some early simulations of circulation in hot Jupiter atmospheres predicted that atmospheric dynamics would naturally produce features like polar vortices, zonal jets, and an asymmetric temperature field. These weather features were predicted to cause variations in a planet's emission properties, depending on the planet's wind speeds and irradiation environment (Cho et al. 2003(Cho et al. , 2008. Rauscher et al. (2008) simulated how these weather features, particularly polar vortices, could cause a planet's phase curve to change over time. They predicted that weather would cause the phase curve to exhibit amplitude modulations upwards of 25% over the course of consecutive orbits.
Other theoretical work, however, has not predicted strong variations driven by the planet's climate (Showman & Guillot 2002;Cooper & Showman 2005;Dobbs-Dixon & Lin 2008). Some models further predicted only minor variations in observable quantities over time. For example, Showman et al. (2009)'s 3D, radiativedynamical models of the hot Jupiters HD 189733b and HD 209458b exhibited only ∼1% variability in the eclipse depth as measured in all Spitzer channels. Dobbs-Dixon & Agol (2013)'s simulations of HD 189733b exhibited very stable dayside emission, varying by less than 0.1% over 30 orbits for all Spitzer channels. Similarly, Menou (2020)'s high-resolution models of HD 209458b exhibited less than 2% variability in the infrared dayside-integrated flux. Finally, Komacek & Showman (2020)'s simulations of phase curve observations of HD 209458b predicted upper limits on eclipse depth variability of 2% and on phase curve amplitude of 1%.
It is worth noting that all of the above models were cloud-free. Adding clouds to a planet's atmosphere is known to change its emission characteristics significantly compared to a cloud-free atmosphere (e.g. Mendonça et al. 2018;Venot et al. 2020;Roman et al. 2021), including inducing variability in the observed emission (e.g. Parmentier et al. 2013;Line & Parmentier 2016;Lines et al. 2018;Powell et al. 2018;Jackson et al. 2019). In addition to clouds, magnetic interactions between a planet's magnetic field and ions embedded in a planet's atmospheric winds may also cause variability. This would be true for hot Jupiters warm enough that they experience significant thermal ionization on their daysides (Rogers 2017;Hindle et al. 2019Hindle et al. , 2021. Early infrared observations of hot Jupiter atmospheres agreed with the lack of variability predicted by cloud-free modeling. Agol et al. (2010)'s repeated secondary eclipse observations of HD 189733b using Spitzer /InfraRed Array Camera (IRAC) at 8 µm placed an upper limit of 2.7% on the variability of HD 189733b's eclipse depth at that wavelength. Crossfield et al. (2012)'s repeated transit and eclipse observations of HD 209458b using Spitzer /MIPS at 24 µm found no significant variations in either the transit or eclipse depth. Kilpatrick et al. (2020) found no variations in the eclipse depths of HD 189733b or HD 209458b at both 3.6 µm and 4.5 µm using Spitzer /IRAC. Similarly, analyses by Wong et al. (2014), Ingalls et al. (2016), and Bell et al. (2021) find no significant evidence for variability in XO-3b's eclipse depth at 4.5 µm .
On the other hand, optical observations of hot Jupiters have seen significant variability. The Kepler phase curve of HAT-P-7b exhibits large variations in its shape and offset over tens to hundreds of days (Armstrong et al. 2016). Similarly, the Kepler phase curve of Kepler-76b shows variation in its eclipse depth of ∼40% on timescales of tens of days (Jackson et al. 2019). Bell et al. (2019) and von Essen et al. (2019) also find variation in WASP-12b's eclipse depth and phase curve offset, but the authors suggest that these changes may be the result of WASP-12b transferring mass to its star, rather than intrinsic atmospheric variations. Long-term TESS observations of WASP-12b have not found significant variability in the planet's eclipse depth over a timescale of couple of years (Wong et al. 2022).
It is unclear whether the observed variation at optical wavelengths is actually due to weather in the planet's atmosphere. Optical observations are more vulnerable to stellar contamination than infrared observations, which can introduce time variability into long-duration time series. For example, Lally & Vanderburg (2022) showed that the observed variability in HAT-P-7b's phase curve can also be explained by supergranulation on the host star.
To date, there have been no published observing programs designed to search for phase curve variability at medium infrared wavelengths, such as those probed by Spitzer /IRAC. Repeated phase curve observations have been taken using Spitzer to help correct the strong systematics that are often present in 3.6 µm data (e.g., WASP-12b (Bell et al. 2019), WASP-76b (May et al. 2021), and WASP-52b (May et al. 2022)) and to measure low-amplitude phase curve signals (e.g., K2-141b (Zieba et al. 2022)), but it is difficult to probe variability with these data sets. Repeat phase curve observations at medium infrared wavelengths present a bridge between existing constraints on variability, and a new test on the variable nature of hot Jupiter atmospheres. We set out to do this by observing two new phase curves of the hot Jupiter WASP-43b using Spitzer /IRAC at 4.5 µm . Together with a previous phase curve from Stevenson et al. (2017), these observations provide three epochs over which to sample WASP-43b's global weather.
WASP-43b is a hot Jupiter orbiting the bright (K=9.2) K7 star WASP-43, discovered in 2011 by Hellier et al. (2011). The planet is on a short 19.5 hour orbital period (Hellier et al. 2011) which drives a dayside brightness temperature of 1485 ± 41 K (May & Stevenson 2020). Previous observations of WASP-43b suggest it has a cloudy nightside (Stevenson et al. 2017;Kataria et al. 2015), and models further predict WASP-43b should have clouds on its dayside as well (Helling et al. 2020). This expected cloudiness, in addition to the planet's short period and high signal-to-noise ratio in phase curve observations (Stevenson et al. 2017), make WASP-43b an ideal target for our investigation into potential infrared variability.

OBSERVATIONS & DATA REDUCTION
We obtained two full-orbit photometric phase curves of WASP-43b at 4.5 µm using Spitzer /IRAC (Fazio et al. 2004) separated by two weeks, on 2019 September 11 and 2019 September 26 (PID 14297). Hereafter, we will refer to these visits as "2019a" and "2019b", respectively. Each visit used the subarray mode with 2.0 second exposures and lasted 25.2 hours. We employed a 30-minute pre-observation before each visit's science observations to allow any detector persistence to build up to a steady state, and excluded this time from our data analysis. The science visits were split into two consecutive 12.5 hour astronomical observing requests (AORs), so that Spitzer would repoint and position WASP-43 onto the IRAC photometric "sweet spot" at the start of each AOR. This re-pointing process added a gap of 6 minutes between each AOR. To stabilize the pointing in each AOR, we used the PCRS peak-up mode targeted on WASP-43 itself.
In addition to our observations, we reanalyzed the Spitzer/IRAC 4.5 µm full-orbit phase curve curve data taken as part of observing program PID 10169 on 2014 August 27, described by Stevenson et al. (2017). Hereafter, we will refer to this visit as "2014". This visit was split into three consecutive AORs. Table 1 provides specific details of each observing visit.
Our data reduction process follows that of Beatty et al. (2018). We started with the basic calibrated data (BCD) images from the Spitzer pipeline. First, we performed background subtraction by masking out the middle of each image, 10 pixels in from each edge, to remove star light, as well as the bottom row of each image since it systematically contains several nan values, and then computed the median value of the unmasked pixels. We then subtracted this median value from each pixel in the image. Then, we corrected for bad pixels by performing three iterations of 5-σ clipping on each pixel's time series and replacing the values of identified bad pixels with the median value of its time series. Using these background-subtracted and bad-pixel-corrected images, we then determined the position of WASP-43 in each image using the Howell centroiding technique (Howell 2006).
We used the background-subtracted images, prior to bad-pixel correction, for the photometric extraction. We extracted the flux within a circular aperture centred on WASP-43's position in each image. Figure 1 shows histograms of WASP-43b's position on the detector for each visit. The aperture size used depended on the data set. We chose the optimal aperture for each visit by sampling an astrophysical model on the data, which is described in Section 3.2, identifying which aperture size gave both the largest Bayesian likelihood and smallest data-model residual. The 2014 and 2019b visits both preferred an aperture radius of 2.4 pixels, while the 2019a visit preferred a radius of 2.1 pixels. For reference, the median FWHMs were 0.50 pixels, 0.75 pixels, and 0.51 pixels for the 2014, 2019a, and 2019b visits, respectively. We then identified and removed any outlier flux measurements in each visit's normalized flux time series by clipping any flux value more than 5-σ from the best-fit model found during aperture optimization. The 2019b light curve exhibited a nonlinear decreasing trend which, in lieu of fitting for systematics at this stage, we fitted with a third-order polynomial in this outlier-cutting model. This removed 59 of 44928 exposures from the 2014 data, 57 of 44928 exposures from the 2019a data, and 132 of 44928 exposures from the 2019b data.
Spitzer time series observations are susceptible to detector persistence, which causes a notable ramp effect at the start of each AOR. Each of our visits showed this initial ramp effect. This effect was identifiable by small increases in the background flux level over the first 5-10 minutes of each AOR, after which the background level reached a constant saturation level throughout the entire visit. The exact duration of this initial persistence ramp was slightly different for each AOR and each visit. We determined the end of each ramp period as the time when the background flux level appeared to flatten off at a constant median value, then cut out all data points prior to this time. This process amounted to cutting 150 points, 174 points, and 73 points from the three AORS of the 2014 visit, 250 points and 206 points from the two AORS of the 2019a visit, and 200 points and 186 points from the two AORs of the 2019b visit. After these two clipping routines, we were left with 44719 points in the We determined the time of each exposure by assuming that each 64 frame image stack began at the MJD OBS header time and the image frames were evenly spaced between the AINTBEG and AINTEND header times. We then converted these spacecraft times to Barycentric Julian Dates in the Barycentric Dynamical Time standard (BJD TDB).

LIGHT CURVE MODEL SELECTION
We modeled the observed light curve as a combination of a time-varying astrophysical signal and systematic effects introduced by the instrument. We denote the astrophysical signal as A(t). The systematic error is a function of both time and spatial position (x,y) on the detector, and we denote it as S(x, y, t). The full light curve model, which we call F obs (t), was then F obs = A(t) S(x, y, t). (1)

Spitzer Systematics
We found that our data displayed the usual systematic effects seen in Spitzer time series. We divided these into intrapixel sensitive variations B(x,y), a linear trend in time f t (t), and linear trends in the x-and y-FWHMs, f x, FWHM and f y, FWHM , respectively. All together, the model for systematics is Intrapixel sensitivity variations, which cause variations in flux at the subpixel scale, are the dominant source of systematic error in Spitzer's 4.5 µm channel. We corrected for these sensitivity variations by using Bilinearly-Interpolated Subpixel Sensitivity (BLISS) mapping (Stevenson et al. 2012) to fit for the relative sensitivity B(x, y) of the detector as a function of position. May & Stevenson (2020) have suggested that a fixed BLISS map, determined from multiple observations, can return the most precise fit, and have published a fixed sensitivity map. However, this map does not cover the full extent of our observed stellar positions. We therefore fit a BLISS map simultaneously with the other model parameters.
In addition to BLISS mapping, previous analyses have shown correlations between the measured flux and the FWHM of the stellar point response function (PRF) (e.g. Mendonça et al. 2018). To measure the FWHM in each frame, we used the Photometry for Orbits, Eclipses, and Transits (POET) pipeline (Stevenson et al. 2012;Cubillos et al. 2013;Stevenson et al. 2017) which fits the stellar point spread function of each image with a Gaussian function. We used these measured FWHMs to fit a linear detrending function that depended on the xand y-width of the star where s F W HM is the slope and w is either the x-or y-FWHM. We used separate slopes for the x-and y-width models and left them as free parameters in the sampling. Lastly, we fit a visit-long linear trend in time of where s t is the slope, which we left as a free parameter in our sampling, and t median is the median time of the visit. We tested this visit-long detrending against one where each AOR was given its temporal trend, and found that the visit-long detrending was preferred in terms of the Bayesian Information Criterion (BIC). Similarly, we tested several other systematic correction methods, including not using FWHM detrending f F W HM , not using flux detrending f t , using a quadratic form of f t , and detrending using noise pixel values, and found that our methods described here were preferred in terms of the BIC.

Astrophysical Model
We modeled the astrophysical signal from the WASP-43 system as the combination of the star's signal (F ), which we treat as just a transit model, and the planet's signal (F planet ), which is a coupled eclipse and phase curve model. Altogether, our astrophysical model was then Here, θ represents the astrophysical parameters for which we fit: the transit center time, orbital period, visit shows considerable overlap, but there is a notable trend in centroids toward the lower right likely due to re-pointing inconsistency between its three AORs. The 2019a visit has full overlap between its two AORs. The 2019b visit has two distinct centroid regions. The top right region consists frames from the first AOR and the bottom left region consists of frames from the second AOR. This suggests the spacecraft was unable to return to its original pointing when starting the second AOR. ratio of the semi-major axis to the stellar radius, inclination, ratio of the planetary and stellar radius, and all phase curve parameters described below. We did not fit for WASP-43b's eccentricity or longitude of periastron. Previous works have shown that WASP-43b's eccentricity is consistent with zero (Hellier et al. 2011;Gillon et al. 2012;Wang et al. 2013;Blecic et al. 2014), so we fixed the eccentricity to zero and the longitude of periastron to ω = π/2.
We modeled the planetary transit F (θ, t) using the BATMAN python package (Kreidberg 2015). We used a quadratic law for stellar limb darkening and fixed the coefficients to u 1 = 0.0796 and u 2 = 0.1621. We estimated these limb darkening coefficients by matching WASP-43's stellar properties in the table of modelled limb darkening coefficients by Claret & Bloemen (2011), using those computed via the least-squares method.
We did not consider any variability on part of the star in our model, which is justified by a lack of stellar variability seen by previous observations of the system as well as the fact that we do not end up finding variability between our observations.
To model the relative flux from the planet F planet (θ, t), we used a combined geometric eclipse and phase curve model The geometric eclipse model, f eclipse , models the relative change in flux during the planet's occultation. We normalize it to be equal to one while out of eclipse and have a depth of unity, so that it physically represents the fraction of the planet's disk which is visible as a function of time. Rather than having to fit for the eclipse depth itself, this method lets the eclipse depth be naturally set by the phase curve parameters. Also, when fitting the eclipse time, we account for WASP-43b's Roemer delay of 15.1 seconds. Our phase curve model represents the emission from the planet as a function of its orbital phase. We use a double harmonic sinusoidal model of the form where Here, f 0 is the vertical phase offset, P is the orbital period, t c is the transit midpoint time, c 1 and c 3 are the phase amplitudes, and c 2 and c 4 are the phase offsets in units of degrees. We add a factor of π in the argument of each cosine term to center the phase curve maximum at the secondary eclipse midpoint for an offset of zero degrees. Stevenson et al. (2017) used a similar formula for their phase curve model and found a nonzero second harmonic in their fit of the 2014 phase curve. We later test whether our fits also show a significant second harmonic.

Fitting Scenarios
Previous observations have ruled out any decay of WASP-43b's orbit (Blecic et al. 2014;Murgas et al. 2014;Hoyer et al. 2016;Patra et al. 2020;Garai et al. 2021), which would naturally lead to changes in the shape of WASP-43b's phase curve over time. In addition, previous observations of this system at shorter wavelengths have not seen significant stellar activity, so we do not expect stellar contamination to be a major concern (Hellier et al. 2011;Gillon et al. 2012;Czesla et al. 2013;Murgas et al. 2014). Therefore, we expect any variability between the phase curves to be the result of intrinsic atmospheric variability. If the planet's climate or weather is causing significant changes to the hemisphere-integrated emission at each phase over time, we should see this as changes in the shape and amplitude of the phase curve. This reasoning motivated us to use three fitting scenarios that differ in which fitting parameters are allowed to vary between visits.
In the first scenario, which we refer to as the "individually fit" case, we fit each visit's light curve individually. We gave each light curve its own set of planetary, orbital, and phase curve parameters, and its own BLISS map. In the second scenario, which we refer to as the "semishared fit" case, we simultaneously fit a shared set of planetary and orbital parameters, and BLISS map, but gave each light curve its own phase curve parameters. In the third scenario, which we refer to as the "fully shared fit" case, we fit all parameters simultaneously so that each light curve used a single, shared set of planetary, orbital, and phase curve parameters, and a shared BLISS map. In all scenarios, we gave each light curve its own set of systematic parameters for f t , f x,F W HM , and f y,F W HM (see Eqns. 3 and 4). As we will show in Section 4, we found no significant variation between visits so this third, "fully shared fit" scenario is preferred. The other two fitting scenarios then served as checks that our individual data sets did not have any uncorrected systematic effects present.
To fit the data in all three scenarios, we employed Markov Chain Monte Carlo (MCMC) sampling using the emcee python package (Foreman-Mackey et al. 2013). We enforced Gaussian priors on only four parameters: the transit center time, orbital period, semi-major axis, and inclination. Within the MCMC sampling, we parametrized the orbital period as log 10 (P ), the semimajor axis as log 10 (a/R ), and the inclination as cos (i), following Eastman et al. (2013). We give these prior values in Table 2. We did not enforce a radius prior since the only radius measurement at this wavelength was by Stevenson et al. (2017), whose data we are reanalyzing in this work. We sampled each scenario for 20,000 steps, which was over 20 times their corresponding autocorrelation times as estimated by emcee, using a number of walkers that was twice the number of free parameters in each case. The 2014 and 2019b individual fits needed a 1,000 step burn-in, the 2019a individual fit needed a 1,500 step burn-in, and the semi-shared and fully shared fits each needed a 5,000 step burn-in.

RESULTS
We found no significant variation between the three observed phase curves. The fully shared fitting case was preferred of our three fitting scenarios, meaning the three visits are best described by a single, shared phase curve rather than three different phase curves.

Fully Shared Fits
As mentioned, in the fully shared fit case, we fit the three visits simultaneously with a single, shared set of planetary, orbital, and phase curve parameters, and a single BLISS map. The best-fit parameter values and relevant derived properties are given in Table 3. The full light curve is plotted in Figure 2, and just the phase curve component is plotted as the orange dotted line in the bottom panel of Figure 4. Also, the raw data and FWHMs from each visit, along with each fitted model, is plotted in Figure 8 of the Appendix.
The fully-shared fit phase curve exhibits a small eastward hotspot offset of 10.4 • ± 1.7 • , which is close to the 12.3 • ± 1.0 • eastward offset of WASP-43b's Hubble Space Telescope (HST)/WFC3 phase curve (Stevenson et al. 2014). We find a maximum phase curve value of 3936 ± 84 ppm, a nonzero minimum phase curve value of 425 ± 105 ppm, a peak-to-peak amplitude of 3511 ± 89 ppm, and an eclipse depth of 3907 ± 85 ppm. Using a BT-Settl model for the star (Allard et al. 2003;Barber et al. 2006;Allard et al. 2007Allard et al. , 2011Caffau et al. 2011;Allard et al. 2012Allard et al. , 2013) and a blackbody model for the planet, we convert these phase curve extrema into integrated brightness temperatures of 1479 ± 13 K on the dayside and 755 ± 46 K on the nightside. The phase curve's second harmonic term does not strongly contribute, as c 3 is close to zero and c 4 has over 50% uncertainty and is not well constrained.  Figure 2. Emission phase curve of WASP-43b at 4.5 µm , which is the fully shared fit in our analysis. The black data points are a combination of the data from the three visits analysed in this work, binned in 10 minute intervals. This binning is only done for presentation purposes, and only unbinned data were used during the analysis. The bottom panel shows the full light curve, while the top panel highlights the phase curve.

Individual and Semi-Shared Fits
We found that the phase curves in the individual fit case and the semi-shared fit case were very consistent with one another. The best-fit parameter values and relevant derived properties are given in Table 4 for the individual fit case, and in Table 5 for the semi-shared fit case. Between each individually fit visit, the primary phase curve parameters (c 1 and c 2 ) all agree within 2.18σ. The 2014 individually fit phase curve does exhibit a non-negligible second harmonic term, with an amplitude that is approximately 11% that of the fundamental term. This is similar to what Stevenson et al. (2017) found for the 2014 data, as their second harmonic amplitude is 14% that of the fundamental term. The 2019a and 2019b phase curves, however, do not show significant second harmonic terms. The semi-shared phase curves are in similar consistency with one another. Between each individual phase curve and its corresponding semi-shared fit phase curve, the parameters are all consistent within 1σ.
We plot the full individually fit light curves in the lefthand column of Figure 3 and show just their phase curve components in the top panel of Figure 4. Likewise, the Phase curve vertical offset 0.00210 ± 0.00008 Note-This fully-shared fitting scenario was our preferred method, and thus these values represent the best-fit values of this work. Values given only in one column with arrows in adjacent columns represent parameters that were shared between the three sets. These values represent the median of the resulting posterior distribution determined from our MCMC sampling. The upper and lower errors given represent the 16th and 84th percentiles of the posterior distribution, respectively.
semi-shared fit light curves are plotted in the middle column of Figure 3 and their phase curve components in the middle panel of Figure 4. The best-fit R p /R values were all consistent within 1.14σ, while the other planetary and orbital parameters were consistent within 0.3σ between each individually fit visit. In addition, they are all within 1.09σ of the semi-shared fit planetary and orbital parameters. As mentioned, all three visits used the same set of these parameters in the semi-shared fit case. The fact that the orbital parameters are all consistent between visits, even when fit entirely individually, indicates that our fits are self-consistent.
When focusing on the phase curve, our fits of the 2014 visit do seem different from the 2019 visits at first glance. The individually fit 2014 visit's minimum phase curve value is -77 ± 200 ppm, which is lower than the 2019 visits' values of 568 ± 153 ppm and 583 ± 330, respectively. The 2014 visit also exhibits a much larger hotspot offset of -19.6 • ± 2.5 • compared to -4.2 • ± 2.5 • for the 2019a visit and -11.4 • ± 6.8 • for the 2019b visit, in the individually fit cases. The discrepancy between the minimum phase curve values is not statistically significant, as they only differ to 2.6σ. The hotspot offsets differ by 4.3σ, which is significant. However, we believe that these discrepancies are the result of strong intrapixel sensitivity effects in the 2014 data, rather than astrophysical variability. We checked for excess red noise left over in the detrended data by plotting the standard deviations of residuals as a function of temporal bin sizes for each visit's light curve, shown in Figure 9 of the Appendix. We find that each visit's residuals roughly follow the 1/ √ N relation expected from pure white noise, and that the 2014 visit does not appear to have significant excess noise compared to our 2019 visits. As we will discuss further in Section 7.1, other works have reanalyzed the 2014 data with different approaches to correcting for the subpixel sensitivity variations and found values that more closely match our results for the 2019 data (Mendonça et al. 2018;Morello et al. 2019;May & Stevenson 2020;Bell et al. 2021).

Comparison of Fits
The three fitting scenarios used span the range of potential variation between visits. If there is significant variation, then one of the two scenarios where the three visits had separate phase curve parameters ought to be preferred. If there is no variation, then the fully shared scenario ought to be preferred.
We compare the goodness-of-fit of each scenario using the BIC, and find that the BIC decisively prefers the fully shared fit. The semi-shared fit came second, with ∆|BIC| = −80 relative to the absolute BIC of the fully shared fit, and the individual fits were least preferred with ∆|BIC| = −112. The fully shared fit phase curve, therefore, represents a single phase curve of the planet, whose emission does not strongly change with time. We plot the corresponding full astrophysical model along with the combined data from all three visits in Figure  2.

GCM Setup
To interpret the observed phase curves, we simulated atmospheric conditions on WASP-43b using a 3D general circulation model (GCM). Our GCM used a doublegray radiative transfer scheme. Specifically, we used the model of Rauscher & Menou (2010, adapted from Hoskins & Simmons (1975), with recent updates to account for radiative feedback and scattering from clouds (Roman & Rauscher 2017Roman et al. 2021). In particular, we used the cloud model of Roman et al. (2021) that includes up to 13 different cloud species for various vertical distributions.
In the Roman et al. (2021) model, clouds are treated as purely temperature-dependent sources of opacity, with constant mixing ratios set by the assumed solar elemental abundances. The optical thicknesses of the clouds are determined by converting the relative molecular abundances (or partial pressures) of each species into particles with prescribed densities and radii, following Roman et al. (2021).
We ran simulations for both clear atmospheres and two different sets of cloud compositions. Following Ro-man et al. (2021), one cloud composition case includes 13 different species: KCl, ZnS, Na 2 S, MnS, Cr 2 O 3 , SiO 2 , Mg 2 SiO 4 , VO, Ni, Fe, Ca 2 SiO 4 , CaTiO 2 , and Al 2 O 3 . We refer to this case using all 13 species as the "all species" case. We also modeled a case in which only eight of these species are included, based on considerations of nucleation efficiency (Gao et al. 2020). This compositional case lacks iron, nickel, and three sulfides, so the clouds are limited to KCl, Cr 2 O 3 , SiO 2 , Mg 2 SiO 4 , VO, Ca 2 SiO 4 , CaTiO 2 , and Al 2 O 3 . We refer to this case with only eight species as the "nucleation limited" case. We refer the reader to (Roman et al. 2021) for more details on the cloud properties and parameters for each species. In general, the silicates form thick, conservatively scattering clouds, while the Fe, Al 2 O 3 , and Cr 2 O 3 form clouds of lower albedos that absorb starlight. Places where clouds overlap have mixed properties, weighted by the optical thickness of each species.
For both cloud composition scenarios, we explored the observational consequences of variations in the cloud deck's vertical thickness. While we expect the vertical extent of the clouds to be determined by a balance between vertical mixing and settling, along with other limiting processes, these processes are computationally expensive and poorly constrained both observationally and theoretically. Rather than attempt to model these cloud physics explicitly, we therefore defined the maximum vertical thickness, relative to the cloud base pressure, after which the cloud is truncated, and run simulations for a range of values. We choose values at fivelayer intervals, ranging from 5 to 45 layers of our 50-layer model. Each cloud species is allowed to have its own individual vertical extent, though all species are subject to the same maximum limit. This effectively mimics a range of vertical mixing strengths and amounts to nine simulations for each of our two cloud composition scenarios. Each 5-layer interval roughly corresponds to one scale height, and we describe each model's maximum vertical thickness in terms of the number of scale heights hereafter.
Our model assumes a constant optical depth per bar of pressure, which results in a constant aerosol mixing ratio for each species. As a result, vertically thicker clouds are correspondingly optically thicker in our model. This directly ties the emission characteristics in our cloudy atmosphere to the vertical extent of the clouds, so that the vertical thickness directly traces the optical thickness of the cloud above the thermal photosphere. In reality, however, optical depth can also be set by the concentration of cloud particles.
Using the GCM and planetary parameters listed in Table 6, we began our simulations with clear skies, no Note-These values represent the median of the resulting posterior distribution determined from our MCMC sampling. The upper and lower errors given represent the 16th and 84th percentiles of the posterior distribution, respectively.
winds, and no horizontal temperature gradients. We ran the simulations for 3500 planetary orbits, assuming tidal synchronization. We allowed the winds, temperatures, and clouds to evolve over the course of the simulation. The results from the GCMs were recorded at 100 day intervals to track possible temporal variability. Additional details of the GCM dynamics can be found in Rauscher & Menou (2012), Roman & Rauscher (2019), and Roman et al. (2021). Once the GCM runs completed, we reduced them to the corresponding emission phase curves at 4.5 µm.

Model Phase Curve Extraction
In order to generate phase curves from the GCM, we followed the routine outlined by Malsky et al. (2021). First, we remapped the GCM from constant pressure levels to 250 constant altitude layers by using vertical hydrostatic equilibrium, as assumed in the GCM. Next, we recovered the locations and properties of each cloud species, for each layer, by using the same criteria as within the GCM. We calculated the particle size, wavelength-dependent extinction efficiency, asymmetry parameter, and single scattering albedo for each cloud species on a 50x50 wavelength and pressure grid following the Mie scattering code from Wolf & Voshchinnikov (2004) and the Rosseland Clouds software package 1 . Lastly, we used the vertical gradient from Roman et al. (2021) to relate the atmospheric pressure to the predicted particle sizes for the radiatively active clouds. At pressures of less than 10 mbar we assumed particles with radii of 0.1 µm, with an exponential increase to particle sizes of approximately 80 µm at 100 bar.
We calculated simulated phase curves from the GCM models using the output 3D atmospheric structure. Using a line-by-line ray-tracing radiative transfer routine that is set up to calculate the intensity toward the observer from each line-of-sight path through the atmosphere (e.g Zhang et al. 2017), we calculated the emission spectrum from the planet, for viewing geometries Note-Values given only in one column with arrows in adjacent columns represent parameters that were shared between the three sets. These values represent the median of the resulting posterior distribution determined from our MCMC sampling. The upper and lower errors given represent the 16th and 84th percentiles of the posterior distribution, respectively. throughout its orbit. We assumed solar abundance (consistent with the GCMs) and recovered the cloud properties as described above. We then integrated the spectrum at each phase over the Spitzer 4.5 micron filter response function to calculate the photometric flux from the planet at each phase. To convert these integrated emission phase curves into relative flux units, we divided by the band-integrated stellar flux computed from a BT-Settl model stellar spectrum of the star WASP-43.

Variability Analysis
Variability due to cloud-driven weather may manifest in multiple ways. For example, the cloud thickness may remain constant but the spatial (i.e. latitudinal and longitudinal) distribution of the clouds at any given phase may change (e.g. Lines et al. 2018). Alternatively, the vertical cloud thickness may change, but the spatial distribution remains roughly constant. Since vertical cloud thickness is a freely prescribed parameter in our models, we are able to investigate this second scenario directly.
We remind the reader that by the "thickness" of each model, we are referring to the maximum vertical extent beyond which the clouds were truncated. By comparing our observations to the extracted phase curves from models of different thicknesses, we seek to constrain the extent to which the thickness of WASP-43b's clouds may have changed between observations while still remaining consistent with the data.
We first checked for temporal variability within each GCM model, by comparing the extracted phase curve at different output times after the initial settling period. Such variability could be caused by changes in the spatial distribution of clouds. We do not see any significant temporal variations in these phase curves, implying that there is likely no significant spatial variation in the cloud coverage at scales larger than the model's spatial resolution of T31 (Table 6).
We also compared our observations to various phase curves extracted from GCMs with different cloud thick- ness values. We computed the reduced-χ 2 statistic between each model phase curve and the detrended data for each individual visit, as well as the three visits combined -excluding data points within transit and eclipse. Figure 5 shows the reduced-χ 2 values for each comparison as a function of the model's maximum cloud thickness.
As shown in Figure 5, no single GCM output perfectly fits the data. We can, however, rule out a wide range of possible cloud thickness models, as all models with thicknesses of 3 to 6 scale heights are strongly rejected with reduced-χ 2 values ranging from 2.6 -12. The remaining models with very thin clouds (1 to 2 scale heights) and very thick clouds (7 to 9 scale heights) are marginally consistent with the detrended observations. Figure 6 shows these remaining model phase curves overplotted with our observations, with all three visits combined, and the fully shared phase curve from section 4.1 that was best fit to these data.
Our observations strongly reject intermediate thicknesses. But, our observations cannot distinguish between the very thin models and the very thick models, as the data are statistically fit equally well by both. This statistical degeneracy is, in part, due to the inability to collect phase curve measurements at transit or eclipse, when we would see only the nightside or dayside, respectively. Rather, the majority of our points of comparison span phases where we see the terminator regions, where the model phase curves are much more similar than they are at day or nightside phases. It does, however, seem like our thin and thick models make fairly similar predictions for the dayside emission. For example, the peak emissions of the one scale height phase curve and the seven scale height phase curve only differ by ∼ 1 data error bar. Cloud feedback provides an explanation for why such different thicknesses could lead to similar phase curves. The thin cloud model has clouds condensing at or below the photosphere leading to a clear, hot atmosphere with strong emission. At intermediate thicknesses, reflective silicate clouds begin condensing above the photosphere which reduces absorption, cools the atmosphere, and lowers emission. The thickest models, however, have built up sufficient vertical extent that highly absorbing clouds with high condensation temperatures finally peak above the photosphere, which reincreases absorption on the dayside, warms the local atmosphere, and re-increases dayside emission. This effect of increased emission due to cloud feedback has also been Figure 4. Overlays of the fitted phase curves from this work. Solid lines are the best fit, while the transparent lines represent the 1σ uncertainty associated with each curve. Top: the individually fit phase curves. Middle: the semi-shared fit phase curves overlayed with the individually fit phase curves. Bottom: all fit phase curves, including the fully shared fit phase curve which we deem representative of all three visits. seen by Roman & Rauscher (2019). It is likely also true that, since our models' clouds are allowed to thin below their maximum prescribed thickness due to evaporation, some of the dayside clouds in our thickest models have evaporated, which will also affect the dayside emission. The lower dayside emission predicted by the thin one scale height model of 4256 ppm is closest to our observed maximum relative flux of 3936 ± 84 ppm. Conversely, the lower nightside emission predicted by the thick models of ∼300 ppm best match our observed minimum relative flux of 425 ± 105 ppm. We believe that, in reality, WASP-43b has a complicated combination of thick nightside clouds and relatively thinner dayside clouds which our models cannot accurately disentangle.
Since neither of the seven, eight, and nine scale height models can be strongly rejected from one another, this thicknesses range of three scale heights places an upper limit on how much the maximum cloud thickness on WASP-43b may have varied between visits. Combining this with the lack of spatial variability seen, our observations suggest that WASP-43b's clouds are stable in terms of their vertical and spatial distributions over timescales as long as several years.

PREDICTIONS FOR JWST OBSERVATIONS
JWST has successfully launched, and has shown great promise to deliver years of exciting infrared science as a successor to the Spitzer Space Telescope. In particular, JWST is well equipped to study atmospheric variability on exoplanets. JWST 's unprecedented precision should be able to probe smaller phase curve amplitude variations, and its wider wavelength range may be able to probe silicate cloud spectral features in the 5 -10 µm range. During the review process of this work, JWST observed a single spectroscopic phase curve of WASP-43b using JWST /MIRI LRS for the JWST Transiting Exoplanet Community Early Release Science Program (hereafter ERS) (Bean et al. 2018). Here, we predict how well such an observation could further constrain the nightside cloud thickness, and the level of thickness variations in WASP-43b's atmosphere.
MIRI LRS covers a wide wavelength range of 5 -12 µm, but achieves its highest precision around 5 -6 µm (Glasse et al. 2015;Kendrew et al. 2015, and our Pan-dExo calculations). For simplicity, we only consider this 5 -6 µm bandpass. Following the same method as in Section 5.2, we extracted cloudy phase curves in this bandpass from our GCMs.
We then simulated MIRI LRS observations of each model phase curve. We used 7172 integrations of 9.86 seconds each, as computed by PandExo (Batalha et al. 2017) for a full-orbit observation of WASP-43b. This PandExo estimate is slightly shorter than the ERS plan to use 8595 integrations of 10.3 seconds each, since they include a longer baseline to test for systematics. We scattered data points around each model phase curve at a 9.86 second cadence, assuming Gaussian noise. We set the 1σ noise level equal to the per-integration flux uncertainty of 603 ppm, which we calculated from our PandExo simulation.
First, we tested whether these faux MIRI LRS observations could significantly distinguish between the thin and thick models, which our Spitzer data could not. For this test, we focused just on the one and seven scale height models with all cloud species. The top panel of Figure 7 shows these one and seven scale height model phase curves and their simulated observations binned into three minute intervals. To test how well we could rule one model out given the other as a truth, we computed the reduced-χ 2 statistic between the simulated observed data of one model and the other model's phase curve. We find that we could distinguish each model from the other at high significance with χ 2 reduced > 2.2. As mentioned in Section 5.3, a primary difference between these thin and thick models is in their prediction of the nightside emission. We repeated the above reduced-χ 2 comparison using only points at phases within ± 0.1 of phase 0 or phase 1 (phase 0 and 1 are equivalent). With this subset of points, we find that we could distinguish the nightside predictions at much stronger significance compared to the full-orbit data, of χ 2 reduced > 3.5. This suggests that the JWST /MIRI LRS phase curves of WASP-43b being taken during the upcoming ERS observations will be able to place a strong constraint on WASP-43b's nightside cloud thickness. We note the caveat, however, that in our models there are other factors in addition to just the prescribed cloud thickness that contribute to the exact amount of emission at any given phase. Then, we tested whether these MIRI LRS observations, if repeated, could place a better limit on thickness variations than the three scale height upper limit made from our Spitzer observations. The bottom panel of Figure 7 shows the 7 -9 scale height model phase curves and their simulated observations, again binned to three minute intervals. We find that, in this 5-6 µm integrated bandpass, there is little difference between the models at any phase and the simulated MIRI LRS observations also cannot distinguish between them. Therefore, repeated MIRI LRS observations would likely not be able to place a better limit on the potential cloud thickness variation in WASP-43b's atmosphere. We find that no model fits the data extremely well, as the minimum reduced-χ 2 = 1.3 for the seven scale height nucleation limited model. However, we find that our observations are reasonably consistent with both the thinnest models (of 1 -2 scale height maximum thicknesses) and the thickest models (7 -9 scale height layer maximum thicknesses). Figure 6. Plot of the cloudy GCM phase curves that are consistent with our combined observations by χ 2 reduced < 2.2. Models with all cloud species are shown as colored solid lines, while those with nucleation limited species are the colored dotted lines. Also shown are the detrended observed data, combined from the 2014, 2019a, and 2019b visits, as fit by our fully shared model. Points within transit and eclipse are removed. The solid black line is the fully shared fit phase curve, which was best fit to the data as described in Section 4.  (Bean et al. 2018) may achieve. Top: here, we focus only on the one and seven scale height models as these represent the two thickness regimes that our Spitzer observations could not distinguish. Assuming one of these models as truth, these MIRI LRS observations would be able to rule out the other model at very high significance. This suggests that MIRI LRS observations will be able to measure accurately WASP-43b's nightside emission and approximate cloud thickness. Bottom: here, we show the 7 -9 scale height models, since our Spitzer observations could not significantly rule out any of these models, thus placing an upper limit on the possible thickness variation. We find that, in this 5 -6 µm bandpass, there is no significant difference between the three models at any phase. Therefore, repeated MIRI LRS observations would likely not be able to place an improved upper limit on the thickness variations.

Reanalyses of WASP-43b's 4.5 micron phase curve
The results of Spitzer /IRAC phase curve observations are generally sensitive to where the star landed on the detector, and the data reduction method used, due primarily to complicated intrapixel sensitivity variations (e.g. see May & Stevenson 2020;May et al. 2022). This systematic error is likely the source of the discrepancy between our individual and semi-shared fits of the 2014 data, and those of the 2019 data.
The first major discrepancy between our 2014 and 2019 phase curves is the minimum phase curve value. In the individually fit case, we found a zero-consistent value of -77 ± 200 for the 2014 visit. The original analysis of the 2014 data also found a near-zero value of 19 ppm (Stevenson et al. 2017;May & Stevenson 2020), and the reanalysis of Bell et al. (2021) also found a low value of ≈196 ppm. Other reanalyses of the 2014 data by Mendonça et al. (2018), Morello et al. (2019), and May & Stevenson (2020) find larger minimum phase curve values of 841 ppm, 300 ± 150 ppm, and 542 ± 124 ppm, respectively. These values more closely match what we found for our 2019 visits, of 568 ± 153 ppm and 583 ± 330 ppm, and suggest that WASP-43b's minimum phase curve value is indeed likely significantly nonzero as should be true for a nightside with a nonzero temperature.
The other major discrepancy is in the hotspot offset. We find a hotspot offset of -19.6 ± 2.5 for the 2014 visit, similar to the -21.1 • ± 1.8 • value found by Stevenson et al. (2017). However, Mendonça et al. (2018), Morello et al. (2019), May & Stevenson (2020), and Bell et al. (2021) find values of -12 • ± 3 • , -11.3 • ± 2.1 • , -20.6 • ± 2.0 • , and -20.4 • ± 3.6 • respectively. These values display significant scatter, and half are consistent with the hotspot offset value of -11.4 • ± 6.8 • that we found for our 2019b visit. These discrepancies in the hotspot offset are likely not astrophysical. As Bell et al. (2021) discuss, the form of phase curve model chosen seems to drive scatter in the exact value of the hotspot offset. Analyses that use, and find contribution from, the second harmonic term tend to find a larger offset than those which use only a single, fundamental term.

Global cloud properties in hot Jupiter atmospheres
When compared to our observations, neither our thinnest or thickest cloud models could be statistically rejected. As discussed in Section 5.3, though, the thinner model provided a closer prediction for the dayside emission to that of our fully shared phase curve, and same of the thick model for the nightside emission. The inability to reject either model is, in part, due to the limited resolution of our model and degeneracy in the dayside emission due to cloud feedback and evaporation. We believe that the most likely scenario is that WASP-43b has a combination of a cloudy nightside and a relatively cloud-free dayside.
Our findings agree with the general picture of WASP-43b's atmosphere painted by previous observations. All observed phase curves of WASP-43b, with HST /WFC3 (Stevenson et al. 2014) and Spitzer /IRAC (Stevenson et al. 2017) 2 , have seen low relative flux on the nightside that cannot be recreated by cloud-free models (e.g. Kataria et al. 2015). Meanwhile, eclipse observations using HST /UVIS (Fraine et al. 2021), CHEOPS, and TESS (Scandariato et al. 2022) have found very low dayside albedos, suggesting a lack of clouds and hazes on the dayside. Helling et al. (2020)'s models predicted significant formation of mineral clouds and photochemical hazes on WASP-43b's dayside, albeit much thinner than on the nightside. It may be the case that these dayside aerosols are sufficiently deep that they have hid from previous observations (see Fraine et al. (2021) for a good discussion of this idea).
Nightside cloud formation is a common and generally undisputed model prediction. The formation and feedback effects of dayside clouds and hazes, however, are less understood. In fact, the similarity we found between the dayside emission of our thinnest and thickest models highlights the immense challenge of accurately modeling cloud formation in exoplanet atmospheres, and interpreting observations thereof. The microphysics involved in cloud formation are difficult to simulate on their own, and the problem becomes extremely complex if trying to include active formation along with effects like radiative feedback, scattering, and absorption in 3D geometries. It is crucial to keep these complexities in mind when attempting to connect model results to spatially unresolved observations.
Our observations and models, both in this work and in general, are resolution limited. Real atmospheres may, and likely do, vary at small spatial scales. Unless models use a very fine spatial grid, which would be difficult in terms of practicality and expense, it will be difficult to study these small scale variations. Additionally, with current observational capabilities, our data are typically limited to disk-integrated quantities so we are inherently insensitive to any atmospheric variations that happen to integrate out to the same observable value. This will likely remain an issue until more spatially resolved observations, such as the technique of eclipse mapping (see, e.g., Rauscher et al. 2007;de Wit et al. 2012;Majeau et al. 2012;Rauscher et al. 2018;Mansfield et al. 2020;Challener & Rauscher 2022), become more common.

Implications of No Variability
Between our 2019 observations, we place upper limits of 210 ppm or 5.0% variation in the eclipse depth, and 128 ppm or 3.7% variation in the peak-to-peak amplitude. These constraints are in line with the smaller scale atmospheric variability predicted by Showman et al. (2009), Dobbs-Dixon & Agol (2013, Komacek & Showman (2020), and Menou (2020). We do not see the strong cloud-induced variability predicted in the models by Lines et al. (2018).
Our empirical limits on the variability in WASP-43b's phase curve are similar to those for HD 189733b (Agol et al. 2010;Kilpatrick et al. 2020) and HD 209458b (Crossfield et al. 2012;Kilpatrick et al. 2020), which were also made using Spitzer. Considered alone, these infrared observations would indicate that hot Jupiter atmospheres are quite stable in time. The significant variability observed at optical wavelengths by Armstrong et al. (2016) and Jackson et al. (2019) are in tension with this picture, assuming they are indeed the result of weather and not stellar contamination. Theory indeed predicts that the observable effect of weather is stronger at shorter wavelengths (Rauscher et al. 2008). There may indeed be a transition from significant to insignificant variability as wavelength increases. In this case, our observations suggest that this transition lies below 4.5 µm . Future higher-precision observations will enable finer constraints on the atmospheric variability.

SUMMARY & CONCLUSIONS
In this work, we analysed two new Spitzer/IRAC 4.5µm phase curves of the hot Jupiter WASP-43b, and reanalyzed the previous 4.5µm phase curve of this planet observed by Stevenson et al. (2017). Altogether, these three full-orbit observations offered three samples to test for variability in WASP-43b's atmosphere.
We found that WASP-43b's three phase curves are best fit by a single, fully-shared model rather than three individual models. This shows that WASP-43b's phase curve does not change in time. The lack of observed variability in these Spitzer observations confirms theoretical predictions that weather-induced variability may be weak (Showman et al. 2009;Dobbs-Dixon & Agol 2013;Komacek & Showman 2020;Menou 2020), and matches the lack of variability seen in longer wavelength Spitzer observations of other hot Jupiters (Agol et al. 2010;Crossfield et al. 2012;Kilpatrick et al. 2020).
To model our observations, we ran a suite of cloudy GCMs. Our models included active cloud formation and were run forward in time to capture potential variability. Each model had a different maximum allowed cloud thickness, beyond which clouds were truncated, to explore the effect of vertical thickness. Our models showed no time variability in their derived phase curves. We did some preliminary modeling to test how other properties vary within our models, such as the total integrated optical depth. This initial work showed a few percent variability in the total optical depth, but more detailed analysis is necessary before drawing any firm conclusions. Furthermore, by comparison to our Spitzer /IRAC observations, we found that the maximum cloud thickness in WASP-43b's atmosphere cannot be varying by more than three pressure scale heights over time. Together, these results suggest that WASP-43b's clouds are stable in their vertical and spatial extents for periods as long as several years.
Our observations were not able to rule out the thinnest cloud models, or the thickest cloud models. We found that the thick models best resembled WASP-43b's nightside, whereas the thin models best resembled the dayside. This suggests that WASP-43b likely has a cloudy nightside and a dayside which is relatively cloud-free at the pressures probed by the observations. Lastly, we simulated phase curve observations of WASP-43b using JWST/MIRI LRS integrated over a 5-6 µm band to predict how the upcoming ERS observations of WASP-43b (Bean et al. 2018) may improve on our constraints made here with Spitzer. Our simulated MIRI LRS observations were able to distinguish the thin and thick cloud models significantly. This implies that the upcoming ERS observations with MIRI LRS will be able to measure precisely WASP-43b's nightside emission, and better constrain the nightside cloud thickness than our Spitzer data could. On the other hand, we found that our simulated MIRI LRS observations could not make any better constraints on the allowed variation in thickness than the Spitzer data. A spectrally resolved analysis, as opposed to our band-integrated approach here, may prove more effective though.
contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. We also made use of, and appreciate the contributors to NASA's Astrophysics Data System and the Spitzer Heritage Archive.
Facilities: Exoplanet Archive, Spitzer (IRAC) Software: astropy (Astropy Collaboration et al. 2013, BATMAN (Kreidberg 2015), emcee (Foreman-Mackey et al. 2013), Matplotlib (Hunter 2007), Numpy (Harris et al. 2020), and the Python programming language. Figure 8 shows our each of our full light curves, with the normalization value given in the caption, and the corresponding best fit astrophysical and systematic models for that visit. Figure 9 shows "Allan variance" plots for each visit, showing how the residual scatter from our best fit models compares to pure white noise expectations. Figure 8. Plots of our full data sets along with each fitted model. The top three rows show the normalized, unbinned light curve of the 2014, 2019a, and 2019b visits, respectively. The normalization factor used for each is given in their corresponding legend. In these top three panels, the red line represents the full best-fit model including astrophysical and systematic components, the cyan line represents the best-fit linear ramp, and the lime line represents the best-fit astrophysical model. The bottom row plots the normalized, unbinned flux against the x-and y-FWHMs of WASP-43 on our images for each visit. The best-fit x-FWHM model is shown in cyan, and the best fit y-FWHM model is shown in lime, for each visit. These best-fit models are all from our fully shared fit case, which provided the best fit. Figure 9. Plot of the standard deviations of the residuals as a function of bin size for each visit, after applying all components of our fully shared fit model. For each data set, the standard deviations decrease roughly as 1 / √ N as the bin size is increased, suggesting there is not much excess red noise left over in the detrended data, and each visit's plot looks nearly the same.