Abstract Exomoons are the natural satellites of planets orbiting stars outside our solar system, of which there are currently no confirmed examples. We present new observations of a candidate exomoon associated with Kepler-1625b using the Hubble Space Telescope to validate or refute the moon’s presence. We find evidence in favor of the moon hypothesis, based on timing deviations and a flux decrement from the star consistent with a large transiting exomoon. Self-consistent photodynamical modeling suggests that the planet is likely several Jupiter masses, while the exomoon has a mass and radius similar to Neptune. Since our inference is dominated by a single but highly precise Hubble epoch, we advocate for future monitoring of the system to check model predictions and confirm repetition of the moon-like signal.

INTRODUCTION The search for exomoons remains in its infancy. To date, there are no confirmed exomoons in the literature, although an array of techniques has been proposed to detect their existence, such as microlensing (1–3), direct imaging (4, 5), cyclotron radio emission (6), pulsar timing (7), and transits (8–10). The transit method is particularly attractive, however, since many small planets down to lunar radius have already been detected (11), and transits afford repeated observing opportunities to further study candidate signals. Previous searches for transiting moons have established that Galilean-sized moons are uncommon at semimajor axes between 0.1 and 1 astronomical unit (AU) (12). This result is consistent with theoretical work that has shown that the shrinking Hill sphere (13) and potential capture into evection resonances (14) during a planet’s inward migration could efficiently remove primordial moons. Nevertheless, among a sample of 284 transiting planets recently surveyed for moons, one planet did show some evidence for a large satellite, Kepler-1625b (12). The planet is a Jupiter-sized validated world (15) orbiting a solar-mass star (16) close to 1 AU in a likely circular path (12), making it a prime a priori candidate for moons. On this basis, and the hints seen in the three transits observed by Kepler, we requested and were awarded time on the Hubble Space Telescope (HST) to observe a fourth transit expected on 28 to 29 October 2017. In this work, we report on these new observations and their impact on the exomoon hypothesis for Kepler-1625b.

MATERIALS AND METHODS Our original analysis was the product of a multiyear survey and thus utilized an earlier version of the processed photometry released by the Kepler Science Operations Center (SOC). In that study (12), we used the simple aperture photometry (SAP) from SOC pipeline version 9.0 (17), but the most recent and final data release uses version 9.3. In this work, we reanalyzed the Kepler data using the revised photometry, which includes updated aperture contamination factors that also affect our analysis. During this process, we also investigated the effect of varying the model used to remove a long-term trend present in the Kepler data. We detrended the revised Kepler photometry using five independent methods. The first method is the CoFiAM (Cosine Filtering with Autocorrelation Minimization) algorithm (18), which was the approach used in the original study, since it was specifically designed with exomoon detection in mind. In addition, we considered four other popular approaches: a polynomial fit, a local line fit, a median filter, and a Gaussian process (see the Supplementary Materials for a detailed description of each). The detrended photometry is stable across the different methods (see Fig. 1), with a maximum standard deviation (SD) between any two SAP time series of 250 parts per million (ppm), far below the median formal uncertainty of ~ 590 ppm. Although we verified that the Presearch Data Conditioning (PDC) version of the photometry (19, 20) produces similar results (as evident in Fig. 1), we ultimately only used the five SAP reductions in what follows. We produced a “method marginalized” final time series by taking the median of the ith datum across the five methods and propagating the variance between them into a revised uncertainty estimate (see the Supplementary Materials for details). In this way, we produced a robust correction of the Kepler data accounting for differences in model assumptions. Fig. 1 Method marginalized detrending. Comparison of five different detrending methods on two different Kepler data products (SAP and PDC). The top curve shows the Kepler reduction used in (12), and the bottom curve shows the method marginalized product used in this work. The three panels show the three transits observed by Kepler. We fit photodynamical models (21) to the revised Kepler data, using the updated contamination factors from SOC version 9.3, before introducing the new HST data. Bayesian model selection revealed only a modest preference for the moon model, with the Bayes factor (K), going from 2 log K = 20.4 in our original study down to just 1.0 now. Detailed investigation revealed that this is not due to our new detrending approach, as we applied our method marginalized detrending to the original version 9.0 data and recovered a similar result to our original analysis (see the Supplementary Materials for details). Instead, it appeared that the reduced evidence was largely caused by the changes in the SAP photometry going from version 9.0 to 9.3 and, to a lesser degree, by the new contamination factors. This can be seen in Fig. 1, where the third transit in particular experienced a pronounced change between the two versions, and it was this epoch that displayed the greatest evidence for a moon-like signature in the original analysis. With a much larger aperture than Kepler, HST is expected to provide several times more precise photometry. Accordingly, the question as to whether Kepler-1625b hosts a large moon should incorporate this new information, and in what follows, we describe how we processed the HST data and then combined them with the revised Kepler photometry. HST monitored the transit of Kepler-1625b occurring on 28 to 29 October 2017 with Wide Field Camera 3 (WFC3). A total of 26 orbits, amounting to some 40 hours, were devoted to observing the event. The observations consisted of one direct image and 232 exposures using the G141 grism, a slitless spectroscopy instrument that projects the star’s spectrum across the charge-coupled device (CCD). This provides spectral information on the target in the near-infrared from about 1.1 to 1.7 μm. Of these 232 exposures, only 3 were unusable, as they coincided with the spacecraft’s passage through the South Atlantic Anomaly, at which time HST was forced to use its less-accurate gyroscopic guidance system. Each exposure lasted roughly 5 min, resulting in about 45 min on target per orbit. Images were extracted using standard tools made available by the Science Telescope Space Institute (STScI) and are described in the Supplementary Materials. Native HST time stamps, recorded in the Modified Julian Date system, were converted to Barycentric Julian Date (BJDUTC) for consistency with the Kepler time stamps. The BJDUTC system accounts for light travel time based on the position of the target and the observer with respect to the solar system barycenter at the time of observation. As the position of HST is constantly changing, we set the position of the observer to be the center of Earth at the time of observation, for which a small discrepancy of ±23 ms is introduced. This discrepancy can be safely ignored for our purposes. While the telescope performed nominally throughout the observation, three well-documented sources of systematic error were present in our data that required removal. First, thermal fluctuations due to the spacecraft’s orbit led to clear brightness changes across the entire CCD (sometimes referred to as “breathing”), which were corrected for by subtracting image median fluxes (see the Supplementary Materials for details). After computing an optimal aperture for the target, we observed a strong intra-orbit ramping effect (also known as the “hook”) in the white light curve (see Fig. 2), which has been previously attributed to charge trapping in the CCD (22, 23). We initially tried a standard parametric approach for correcting these ramps using an exponential function but found the result to be suboptimal. Instead, we devised a new nonparametric approach described in the Supplementary Materials that substantially outperformed the previous approach. Fig. 2 Hook corrections. (Top) The optimal aperture photometry of our target (left) and the best comparison star (right), where the hooks and visit-long trends are clearly present. Points are colored by their exposure number within each HST orbit (triangles represent outliers). (Middle) A hook correction using the common exponential ramp model on both stars. (Bottom) The result from an alternative and novel hook correction approach introduced in this work. The intra-orbit root mean square (RMS) value is quoted for the hook-corrected light curves. We achieved a final mean intra-orbit precision of 375.5 ppm (versus 440.1 ppm using exponential functions), which was about 3.8 times more precise than Kepler when correcting for exposure time. The transit of Kepler-1625b was clearly observed even before the hook correction. After removal of the hooks, an apparent second decrease in brightness appeared toward the end of the observations, which was evident even in the noisier exponential ramp corrected data (see Fig. 2). Repeating our analysis for the only other bright star fully on the CCD, KIC 4760469, revealed no peculiar behavior at this time, indicating that the dip was not due to an instrumental common mode. Similarly, the centroids of both the target and the comparison star showed no anomalous change around this time (see fig. S6 in the Supplementary Materials). A detailed analysis of the centroid variations of both the target and the comparison star revealed that the 10-millipixel motion observed was highly unlikely to be able to produce the ~500 ppm dip associated with the moon-like signature. Further, we found that the signal was achromatic-appearing in two distinct spectral channels, which was consistent with expectations for a real moon. Finally, a detailed analysis of the photometric residuals revealed that the fits including a moon-like transit were consistent with uncorrelated noise equal to the value derived from our hook correction algorithm. These three tests, detailed in the Supplementary Materials, provide no reason to doubt that the moon-like dip is astrophysical in nature and thus we treat it as such in what follows. Upon inspection of the HST images, we identified a previously uncataloged point source within 2 arcseconds of our target. The star resides at position angle 8.5° east of north, with a derived Kepler magnitude of 22.7. We attribute its new identification to the fact that it is both exceptionally faint and so close to the target that it was always lost in the glare in other images. Using a Gaia-derived distance to the target, we found that, were this point source to be at the same distance, it would be within 4500 AU of Kepler-1625. However, it is not known whether the two sources are physically associated. Its faintness means that it produces negligible contamination to our target spectrum. We estimated that the source has a variability of 0.33% and contributes less than 1 part in 3000 to our final WFC3 white light curve, which means that the net contribution to our target is 1 ppm and can be safely ignored. In addition to the breathing and the hooks, a third well-known source of WFC3 systematic error we see is a visit-long trend (apparent in Fig. 2). These trends have not yet been correlated to any physical parameter related to the WFC3 observations (24), and thus, the conventional approach is a linear slope [for example, (25–27)], although a quadratic model has been used in some instances [for example, (28, 29)]. The time scale of the variations is comparable to the transit itself and thus cannot be removed in isolation; rather, any detrending model is expected to be covariant with the transit model. For this reason, it was necessary to perform the detrending regression simultaneous to the transit model fits. We considered three possible trend models: linear, quadratic, and exponential. All models include an extra parameter describing a flux offset between the 14th and 15th orbits. This is motivated by the fact that the spacecraft performed a full guide star acquisition at the beginning of the 15th orbit (a new “visit”) and ended up placing the spectrum ~0.1 pixels away from where it appeared during the first 14 orbits. Although the white light curve shows no obvious flux change at this time, the reddest channels display substantial shifts motivating this offset term. Finally, we extracted light curves in nine wavelength bins across the spectrum in an attempt to perform transmission spectroscopy. As a planet transits its host star, the atmosphere may absorb different amounts of light depending on the constituent molecules and their abundances (30). This makes the planet’s transit depth wavelength-dependent. An accurate measurement of these transit depths not only provides the potential to characterize the atmosphere’s composition; it is also potentially useful in providing an independent measurement of the planet’s mass (31). While a low–surface gravity planet will show very pronounced molecular features and a steep slope at short wavelengths due to Rayleigh scattering, a high–surface gravity world will yield a substantially flatter transmission spectrum. With the HST WFC3 data prepared, we are ready to combine them with the revised Kepler data to regress candidate models and compare them. We considered four different transit models, which, when combined with three different visit-long trend models, leads to a total of 12 models to evaluate. The four transit models here were designated as P, for the planet-only model; T, for a model that fits the observed transit timing variations (TTVs) in the system agnostically; Z, for the zero-radius moon model, which may produce all the gravitational effects of an exomoon without the flux reductions of a moon transit; and M, which is the full planet plus moon model. Models were generated using the LUNA photodynamical software package (21), and regression was performed via the multimodal nested sampling algorithm MultiNest (32, 33). For each model, we derived not only the joint a posteriori parameter samples but also a Bayesian evidence (also known as the marginal likelihood) enabling direct calculation of the Bayes factor between models.

RESULTS One clear result from our analysis is that the HST transit of Kepler-1625b occurred 77.8 min earlier than expected, indicating TTVs in the system. Bayes factors between models P and T support the presence of significant TTVs for any choice of detrending model (see Table 1), with the T fits returning a χ2 decreased by 17 to 19 (for 1048 data points). Further, if we fit the Kepler data in isolation and make predictions for the HST transit time, the observed time is > 3 σ discrepant (see fig. S12 in the Supplementary Materials). For reference, each Kepler transit midtime has an uncertainty on the order of 10 min, and the SD on linear ephemeris predictions is 25.2 min derived from posterior samples. Identifying TTVs was among the first methods proposed to discover exomoons (8), but certainly, perturbations from an unseen planet could also be responsible. We find that the ≃25-min amplitude TTV can be explained by an external perturbing planet (see the Supplementary Materials), although with only four transits on hand, it is not possible to constrain the mass or location of such a planet, and no other planet has been observed so far in the system. Table 1 Model performance. Bayesian evidences (Z) and maximum likelihoods ( ) from our combined fits using Kepler and new HST data. Kepler plus HST fits. The subscripts are P for the planet model, T for the planetary TTV model, Z for the zero-radius moon model, and M for the moon model. The three columns are for each trend model attempted. View this table: We also found that model Z consistently outperforms model T, though the improvement to the fits is smaller at Δχ2 ≃ 2−5 (see Table 1). This suggests that the evidence for the moon based on timing effects alone goes beyond the TTVs, providing modest evidence in favor of additional dynamical effects such as duration changes (9) and/or impact parameter variation (10), both expected consequences of a moon present in the system. This by itself would not constitute a strong enough case for a moon detection claim, but we consider it to be an important additional check that a real exomoon would be expected to pass. The most compelling piece of evidence for an exomoon would be an exomoon transit, in addition to the observed TTV. If Kepler-1625b’s early transit were indeed due to an exomoon, then we should expect the moon to transit late on the opposite side of the barycenter. The previously mentioned existence of an apparent flux decrease toward the end of our observations is therefore where we would expect it to be under this hypothesis. Although we have established that this dip is most likely astrophysical, we have yet to discuss its significance or its compatibility with a self-consistent moon model. We find that our self-consistent planet plus moon models (M) always outperform all other transit models in terms of maximum likelihood and Bayesian evidences (see Table 1). The moon signal is found to have a signal-to-noise ratio of at least 19. The presence of a TTV and an apparent decrease in flux at the correct phase position together suggest that the exomoon is the best explanation. However, as is apparent from Fig. 3, the amplitude and shape of the putative exomoon transit vary somewhat between the trend models, leading to both distinct model evidences and associated system parameters. Fig. 3 HST detrending. The HST observations with three proposed trends fit to the data (left) and with the trends removed (right). Bottom-right numbers in each row give the Bayes factor between a planet plus moon model (model M) and a planet plus moon model where the moon radius equals zero (model Z), which tracks the significance of the moon-like dip in isolation.

DISCUSSION Although the overall preference of the moon model is arguably best framed by comparison to model P, the significance of the moon-like transit alone is best framed by comparing M and Z alone. Such a comparison reveals a strong dependency of the implied significance on the trend model used. In the worst case, we have the quadratic model with 2 log K ≃ 4, corresponding to “positive evidence” (34), although we note that the absolute evidence is the worst among the three. The linear model is far more optimistic, yielding 2 log K ≃ 18, corresponding to “very strong evidence” (34), whereas the exponential sits between these extremes. The question then arises, which of our trend models is the correct one? Because the linear model is a nested version of the quadratic model, and both models are linear with respect to time, it is more straightforward to compare these two. The quadratic model essentially recovers the linear model, apparent from Fig. 3, with a curvature within 1.5 σ of zero, and yields almost the same best χ2 score to within 1.2. This lack of meaningful improvement causes the log evidence to drop by 2.8, since evidences penalize wasted prior volume. The exponential model appears more competitive with a log evidence of 1.72 lower, but a direct comparison of two different classes of models, such as these, is muddied by the fact that these analyses are sensitive to the choice of priors. The most useful comparison here is simply to state that the maximum likelihoods are within Δχ2 = 0.68 of one another and thus are likely equally justified from a data-driven perspective. Another approach we considered is to weigh the trend models using the posterior samples. Given a planet or moon’s mass, there is a probabilistic range of expected radii based on empirical mass-radius relations (35). Although we exclude extreme densities in our fits, parameters from model M can certainly lead to improbable solutions with regard to the photodynamically inferred (36) masses and radii. To investigate this, we inferred the planetary mass using two methods for each model and evaluated their self-consistency. The first method combines the photodynamically inferred planet-to-star mass ratio (36) with a prediction for the mass based on the well-constrained radius using forecaster, an empirical probabilistic mass-radius relation (35). The second method approaches the problem from the other side, taking the moon’s radius and predicting its mass with forecaster and then calculating the planetary mass via the photodynamically inferred moon-to-planet mass ratio. Our analysis (discussed in more detail in the Supplementary Materials) reveals that all three models have physically plausible solutions and generally converge at ~103 M ⊕ for the planetary mass, with the exception of the quadratic model that had broader support extending down to Saturn mass. We ultimately combined the two mass estimates to provide a final best estimate for each model in Table 2. Table 2 System parameters. Median and ±34.1% quantile range of the a posteriori model parameters from model M, where each column defined a different visit-long trend model. The top panel gives the credible intervals for the actual parameters used in the fit, and the lower panel gives a selection of relevant derived parameters conditioned upon our revised stellar parameters. The quoted inclination of the satellite is the inclination modulo 90°. View this table: As a consistency check, we used our derived transmission spectrum to constrain the allowed range of planetary masses for a cloudless atmosphere (31). Using an MCMC (Markov chain Monte Carlo) with Exo-Transmit (37), we find that masses in the range of > 0.4 Jupiter masses (to 95% confidence) are consistent with the nearly flat spectrum observed, assuming a cloudless atmosphere (see the Supplementary Materials for details). In conclusion, the linear and exponential models appear to be the most justified by the data and also lead to slightly improved physical self-consistency, although we certainly cannot exclude the quadratic model at this time. For this reason, we elected to present the associated system parameters resulting from all three models in Table 2. The maximum a posteriori solutions from each, using model M, are presented in Fig. 4 for reference. Fig. 4 Moon solutions. The three transits in Kepler (top) and the October 2017 transit observed with HST (bottom) for the three trend model solutions. The three colored lines show the corresponding trend model solutions for model M, our favored transit model. The shape of the HST transit differs from that of the Kepler transits owing to limb darkening differences between the bandpasses. We briefly comment on some of the inferred physical parameters for this system. First, we have assumed a circular moon orbit throughout due to the likely rapid effects of tidal circularization. However, we did allow the moon to explore three-dimensional orbits and find some evidence for noncoplanarity. Our solution somewhat favors a moon orbit tilted by about 45° to the planet’s orbital plane, with both pro- and retrograde solutions being compatible. The only comparable known large moon with such an inclined orbit is Triton around Neptune, which is generally thought to be a captured Kuiper Belt object (38). However, we caution that the constraints here are weak, reflected by the posterior’s broad shape, and thus it would be unsurprising if the true answer is coplanar. One jarring aspect of the system is the sheer scale of it. The exomoon has a radius of ≃4 R ⊕ , making it very similar to Neptune or Uranus in size. The measured mass, including the forecaster constraints, comes in at log(M S /M ⊕ ) = (1.2 ± 0.3), which is again compatible with Neptune or Uranus (although note that this solution is in part informed by an empirical mass-radius relation). This Neptune-like moon orbits a planet with a size fully compatible with that of Jupiter at (11.4 ± 1.5) R ⊕ , but most likely a few times more massive. Finally, although the moon’s period is highly degenerate and multimodal, we find that the semimajor axis is relatively wide at ≃40 planetary radii. With a Hill radius of (200 ± 50) planetary radii, this is well within the Hill sphere and expected region of stability (see the Supplementary Materials for further discussion). The blackbody equilibrium temperature of the planet and moon, assuming zero albedo, is ~350 K. Adopting a more realistic albedo can drop this down to ~300 K. Of course, as a likely gaseous pair of objects, there is not much prospect of habitability here, although it appears that the moon can indeed be in the temperature zone for optimistic definitions of the habitable zone. What is particularly interesting about the star is that it appears to be a solar-mass star evolving off the main sequence. This inference is supported by a recent analysis of the Gaia DR2 parallax (39), as well as our own isochrone fits (see the Supplementary Materials). We find that the star is certainly older than the Sun, at ≃9 gigayears in age, and that insolation at the location of the system was thus lower in the past. The luminosity was likely close to solar for most of the star’s life, making the equilibrium temperature drop down to ~250 K for Jovian albedos for most of its existence. The old age of the system also implies plenty of time for tidal evolution, which could explain why we find the moon at a fairly wide orbital separation. The origins of such a system can only be speculated upon at this time. A mass ratio of 1.5% is certainly not unphysical from in situ formation using gas-starved disk models, but it does represent the very upper end of what numerical simulations form (40). In such a scenario, a separate explanation for the tilt would be required. Impacts between gaseous planets leading to captured moons are not well studied but could be worth further investigation. A binary exchange mechanism would be challenged by the requirement for a Neptune to be in an initial binary with an object of comparable mass, such as a super-Earth (38). Formation of an initial binary planet, perhaps through tidal capture, seems improbable due to the tight orbits simulation work tends to produce from such events (41). If confirmed, Kepler-1625b-i will certainly provide an interesting puzzle for theorists to solve.

CONCLUSION Together, a detailed investigation of a suite of models tested in this work suggests that the exomoon hypothesis is the best explanation for the available observations. The two main pieces of information driving this result are (i) a strong case for TTVs, in particular a 77.8-min early transit observed during our HST observations, and (ii) a moon-like transit signature occurring after the planetary transit. We also note that we find a modestly improved evidence when including additional dynamical effects induced by moons aside from TTVs. The exomoon hypothesis is further strengthened by our analysis that demonstrates that (i) the moon-like transit is not due to an instrumental common mode, residual pixel sensitivity variations, or chromatic systematics; (ii) the moon-like transit occurs at the correct phase position to also explain the observed TTV; and (iii) simultaneous detrending and photodynamical modeling retrieves a solution that is not only favored by the data but is also physically self-consistent. Together, these lines of evidence all support the hypothesis of an exomoon orbiting Kepler-1625b. The exomoon is also the simplest hypothesis to explain both the TTV and the post-transit flux decrease, since other solutions would require two separate and unconnected explanations for these two observations. There remain some aspects of our present interpretation of the data that give us pause. First, the moon’s Neptunian size and inclined orbit are peculiar, though it is difficult to assess how likely this is a priori since no previously known exomoons exist. Second, the moon’s transit occurs toward the end of the observations and more out-of-transit data could have more cleanly resolved this signal. Third, the moon’s inferred properties are sensitive to the model used for correcting HST’s visit-long trend, and thus, some uncertainty remains regarding the true system properties. However, the solution we deem most likely, a linear visit-long trend, also represents the most widely agreed upon solution for the visit-long trend in the literature. Finally, it is somewhat ironic that the case for observing Kepler-1625b with HST was contingent on a previous data release of the Kepler photometry that indicated a moon (12), while the most recent data release only modestly favors that hypothesis when treated in isolation. Despite this, we would argue that planets such as Kepler-1625b— Jupiter-sized planets on wide, circular orbits around solar-mass stars— were always ideal targets for exomoon follow-up. There are certainly hints of the moon present even in the revised Kepler data, but it is the HST data—with a precision four times superior to Kepler—that are critical to driving the moon as the favored model. These points suggest that it would be worthwhile to pursue similar Kepler planets for exomoons with HST or other facilities, even if the Kepler data alone do not show large moon-like signatures. Furthermore, our work demonstrates how impactful the changes to Kepler photometry were, at least in this case, as it suggests that other results over the course of the Kepler mission may be similarly affected, particularly for small signals. All in all, it is difficult to assign a precise probability to the reality of Kepler-1625b-i. Formally, the preference for the moon model over the planet-only model is very high, with a Bayes factor exceeding 400,000. On the other hand, this is a complicated and involved analysis where a minor effect unaccounted for, or an anomalous artifact, could potentially change our interpretation. In short, it is the unknown unknowns that we cannot quantify. These reservations exist because this would be a first-of-its-kind detection—the first exomoon. Historically, the first exoplanet claims faced great skepticism because there was simply no precedence for them. If many more exomoons are detected in the coming years with similar properties to Kepler-1625b-i, it would hardly be a controversial claim to add one more. Ultimately, Kepler-1625b-i cannot be considered confirmed until it has survived the long scrutiny of many years, observations and community skepticism, and perhaps the detection of similar such objects. Despite this, it is an exciting reminder of how little we really know about distant planetary systems and the great spirit of discovery that exoplanetary science embodies.

SUPPLEMENTARY MATERIALS Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/4/10/eaav1784/DC1 Supplementary Materials and Methods Fig. S1. The “Phantom” star. Fig. S2. Kepler detrending. Fig. S3. Kepler detrending comparison. Fig. S4. HST image rotation. Fig. S5. HST hook model comparison. Fig. S6. HST centroids. Fig. S7. Spectral analysis. Fig. S8. Transmission spectrum. Fig. S9. Wavelength solution. Fig. S10. Wavelength-dependent pixel sensitivity. Fig. S11. Modeling the uncataloged source contamination. Fig. S12. Transit timing variations. Fig. S13. Residual analysis. Fig. S14. Chromatic test. Fig. S15. Mass constraints. Fig. S16. Model posteriors. Fig. S17. Physical posteriors. Fig. S18. The May 2019 transit. Table S1. Kepler-only fits. Table S2. Transmission spectrum. Table S3. Transit timings. References (42–72)

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.