No Evidence for Lunar Transit in New Analysis of Hubble Space Telescope Observations of the Kepler-1625 System
Observations of the Kepler-1625 system with the Kepler and Hubble Space Telescopes have suggested the presence of a candidate exomoon, Kepler-1625b I, a Neptune-radius satellite orbiting a long-period Jovian planet. Here we present a new analysis of the Hubble observations, using an independent data reduction pipeline. We find that the transit light curve is well fit with a planet-only model, with a best-fit equal to . The addition of a moon does not significantly improve the fit quality. We compare our results directly with the original light curve from [teachey18b], and find that we obtain a better fit to the data using a model with fewer free parameters (no moon). We discuss possible sources for the discrepancy in our results, and conclude that the lunar transit signal found by [teachey18b] was likely an artifact of the data reduction. This finding highlights the need to develop independent pipelines to confirm results that push the limits of measurement precision.
Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138 \move@AU\move@AF\@affiliationHarvard Society of Fellows, 78 Mount Auburn Street, Cambridge, MA 02138 \move@AU\move@AF\@affiliationFlatiron Institute, Simons Foundation, 162 Fifth Ave, New York, NY 10010, USA \move@AU\move@AF\@affiliationFlatiron Institute, Simons Foundation, 162 Fifth Ave, New York, NY 10010, USA
Moons are abundant in the Solar System, and provide clues to the formation history, evolution, and even habitability of the planets they orbit. The great scientific potential of moons has prompted extensive search for lunar companions in exoplanetary systems (exomoons), and creative development of new search techniques (e.g. kipping09a; kipping09b; kipping13; simon10; peters13; heller14; noyola14; hippke15; agol15; sengupta16; vanderburg18).
Recently, a potential exomoon candidate was identified in the Kepler-1625 system (teachey18a). The host planet, Kepler-1625b, has a radius consistent with that of Jupiter and an orbital period of 287 days. The first evidence for the exomoon candidate, Kepler-1625b I, was based on observations from Kepler. The light curve showed drops in stellar flux that were attributed to a transiting exomoon; however, later analysis called this result into question, showing that the moon transit features were highly sensitive to the Kepler reduction pipeline and the algorithm used to detrend the data (teachey18b; rodenbeck18).
Subsequent follow-up observations with HST revived the possibility of an exomoon in the system based on two factors: a small drop in the system flux after the planet’s transit egress and a transit timing variation (teachey18b, hereafter TK18). The best fit moon had a large radius (comparable to that of Neptune), and if real, is unlike any moon in the Solar System.
Since the primary evidence for the exomoon now rests on the HST transit light curve, in this work we perform an independent reduction and fit to the HST data and compare it to the results from TK18.
2 Observations and Data Reduction
The Kepler-1625 system was observed with 26 continuous HST orbits on 28 - 29 October, 2017 (Program GO 15149: PI: A. Teachey). The observations used the Wide Field Camera 3 (WFC3) G141 grism in staring mode, which fixed the spectrum in a constant position on the detector. At the beginning of the visit, there was a single exposure taken with the F130N filter, which is used to determine the position of the spectral trace. The following exposures used the G141 grism. See TK18 for additional description of the observation design.
We reduced the HST data using custom software developed in kreidberg14a. This software has yielded consistent results with multiple independent pipelines (e.g. knutson14b; spake18). We ran our pipeline on the flt data product provided by the Space Telescope Science Institute (STScI). In keeping with previous WFC3 analysis, we discarded the first orbit of data, where the instrument systematics have larger amplitude. We also discarded exposures taken during the South Atlantic Anomaly passage (exposures 107, 116, 125, and 126).
To begin the data reduction, we fit the centroid of the direct image with a two-dimensional Gaussian. The centroid position determines the position of the spectral trace, which we calculated using the coeffients provided in the configuration file from STScI: G141.F130N.V4.32.conf111available at http://www.stsci.edu/hst/wfc3/analysis/grism˙obs/calibrations/wfc3˙g141.html. To process the spectra, we flatfielded the raw data using the spectroscopic flatfield coefficients provided by STScI in WFC3.IR.G141.flat.2.fits, following the instructions in Section 6 of the aXe User Manual222http://ane-info.stsci.edu/. We then created an extraction box centered on the spectral trace. We varied the height and width of the box in 1-pixel increments to find the window that minimized the root-mean-square (rms) deviation from the best fit to the transit light curve. The best was , and , where X and Y are physical pixels in the spectral and spatial direction, respectively.
We reduced the grism exposures with the optimal extraction routine of horne86, which minimizes background noise in the extracted spectrum by preferentially weighting pixels that are dominated by the target spectrum. The inputs for optimal extraction are the background-subtracted data array, and estimate of the error per pixel, an initial guess for the spectrum and its uncertainty, and a mask array for bad pixels. For the initial guess of the spectrum, we did a simple box extraction (sum over all rows in the extraction window), and assumed the variance was equal to the box-extracted spectrum. We subtracted the background from the data array as described in 2.1. For the error array, we used a quadrature sum of the photon noise (the square root of the pixel counts), the read noise (12 photoelectrons for flt files; WFC3 Data Handbook333http://www.stsci.edu/hst/wfc3/documents/handbooks/currentDHB/), and the error due to background subtraction. The initial pixel mask marked all pixels as good.
To optimally extract the spectrum, we first created a smoothed image by median-filtering each row of the data with a 9-pixel-wide window. We then normalized the smoothed image by dividing each column by its sum, and multiplied it by the best guess spectrum. We compared the smoothed image to the real data and masked outliers in the data that are greater than a threshold . We then recomputed the best guess spectrum with the new mask and the optimal weights from horne86. The process is iterated until no outliers greater than the threshold remain. To create the broadband transit light curve, we sum each optimally extracted spectrum over all wavelengths.
The broadband light curve is shown in Figure 2.1, in comparison to the light curve from TK18. We note that there are differences between the two data sets, particularly a kink near the moon-like transit feature identified by TK18.
2.1 Background Subtraction
The star Kepler-1625 is faint (H mag = 14.0) relative to most other exoplanet host stars observed with WFC3, which makes accurate background subtraction especially important for this target. Moreover, the host star is in a crowded field, so the pixels used to estimate the background must be chosen carefully to avoid contamination from other stars. To estimate the background counts, we masked pixels with total counts larger than 800 electrons (2.7 electrons/sec) and took the median count in the unmasked pixels. The per pixel uncertainty due to background subtraction is 1.4826 times the median absolute deviation.
2.2 Pointing Drift Measurement
The position of the spectrum on the detector shifts slightly over time ( pixel/day) due to the spacecraft’s pointing drift. This drift can change the flux measured for the target star: if the spectrum moves onto less sensitive pixels, fewer photoelectrons are recorded. To enable a correction for this effect, we measured the position of the spectrum over time. Figure 2.2 shows the best fit shifts.
To measure shifts in the spatial direction, we summed each flt image over all columns (which we dub the “column sum”). We used the first exposure in the visit as a template, and for each subsequent exposure, we used least-squares minimization to calculate the shift in pixels that minimized the difference between its column sum and the template. The shifts are a fraction of a pixel, so we used the NumPy interp routine to do linear interpolation on a sub-pixel scale. The WFC3 point spread function is undersampled, so we convolved each column sum with a 4-pixel-wide Gaussian before the interpolation (following deming13).
To measure the spectral shifts, we repeated this procedure with two differences: (1) we used the optimally extracted spectrum rather than the column sum; and (2) in addition to calculating the best fit shift, we also calculated a best fit normalization factor (a scalar multiple for the whole spectrum), to ensure that our results are not biased by the varying brightness of the host star during the planet’s transit.
The extracted transit light transit light curve contains both astrophysical signals and instrument systematic noise, which we model simultaneously.
3.1 Astrophysics Model
For the astrophysics, we used the planetplanet package (luger17), a photodynamical code that calculates light curves for multiple occulting bodies orbiting a star. Within planetplanet, the orbits are computed with the N-body integrator REBOUND (rein12), which calculates the three-dimensional motion of the star, planet, and moon over time under the influence of gravity.
In our analysis, we considered two scenarios: a no-moon model and a moon model. The free parameters for the no-moon model were: the stellar radius, the planet radius, the planet’s time of central transit, the planet inclination, and the planet mass. For the moon model, we added a third body with six free parameters: moon radius, moon time of central transit, moon orbital period, moon inclination, moon mass, and longitude of the ascending node relative to that of the planet. We fixed the eccentricity of all bodies to zero. We also fixed the orbital period of the planet-moon barycenter to 287.378949 days (the best fit from TK18). We elected not to vary the period of the planet-moon barycenter because it is poorly constrained from a single transit observation. A longer period would cause a longer transit duration for the planet, but this could equivalently result from a smaller impact parameter, a larger star, or a massive moon that significantly perturbs the planetary orbit.
We used the following priors: the stellar radius was normally distributed, (see next section). The planet radius was uniform from . The planet transit time was uniform over the timespan of the observations, and inclination was uniform from . The planet mass was log-normally distributed, , based on the expectation for Jupiter-radius objects from ning18. For the moon model, we allowed the transit time to vary uniformly over the entire visit. The moon period was uniform between 1.6 to 260 days. These limits span the duration of the HST observations (so there is one possible moon occultation event), to the orbit at 0.5 the Hill radius, based on the stability limit for prograde moon orbits (domingos06). The Hill radius calculation assumed the planet and stellar masses are and . The moon mass varied uniformly from . The longitude of the ascending node was also uniform from . The moon inclination was uniform from , and was defined relative to the line of sight. We assigned zero prior probability to scenarios where the moon did not transit or experienced a grazing transit. We made this choice to put an upper limit on the radius of a fully transiting moon; for grazing transits or non-transits the moon radius could be arbitrarily large.
3.1.1 Stellar Parameters
For both the moon and no-moon scenarios, we used a quadratic stellar limb darkening law and fixed the coefficients to the prediction for a 5700 K, solar metallicity PHOENIX model from espinoza15; .
We estimatea the host star parameters using the Gaia DR2 parallax (Gaia; GaiaDR2) along with UBV photometry from Everett2012 and JHK photometry from 2MASS (2MASS). We employed the isochrone python package (isochrone) with the Dartmouth isochrone grid (Dotter2008) to obtain posterior constraints on the stellar parameters. The resulting parameters indicate that Kepler-1625 has stellar mass M, radius R, and age Gyr. In our analysis, we fixed the stellar mass to the best fit value (), and used a Gaussian prior on the radius, .
3.2 Instrument Systematics Model
There are two systematic trends in the data. One is the orbit-long ramp, attributed to charge traps in the detector filling up over the orbit (zhou17). The other is a visit-long trend over multiple orbits, which could be due to shifts in the target star position onto more/less sensitive pixels.
For our primary fit, we used a linear combination of the spectrum’s X and Y position (shown in Figure 2.2), multiplied by the non-parametric orbital ramp model from TK18, which assigns each of the nine exposures per orbit a normalization constant, . In sum, for exposure number , the systematics model S is:
where and , and are free parameters, and is the exposure number relative to the first exposure in the orbit.
For comparison, we also tested the “exponential” systematics model from TK18, which combines an exponential visit-long trend, an offset after orbit 14 when the guide stars were reacquired, and the non-parametric orbital ramp model.
3.3 Light Curve Fits
We fit the broadband transit light curve using the models described above. We determined the best fit model parameters with least-squares minimization. We also ran a Markov chain Monte Carlo (MCMC) fit to determine the posterior probability of the parameters. For the MCMC, we held the ramp parameters fixed at their best-fit values. The MCMC used the emcee package (foremanmackey13) with 50 walkers and ran for 5000 steps. We discarded the first 20% of the MCMC chain as burn-in. As a quick test for convergence, we divided the remainder of the chain in half and confirmed that the results from the first half were consistent with the second half.
We obtained an excellent fit to the light curve with the no-moon model, as illustrated in Figure 3.1. The residuals to the no-moon model fit have rms equal to 356 parts per million (ppm), which is within 3% of the predicted photon shot noise (367 ppm), and yields a . The binned rms decreases with the square root of the number of points per bin, as expected for photon noise-limited statistics (see rms versus bin size in Figure 4).
We obtained a slightly lower rms with the moon model (351 ppm); however, this is not a large enough improvement in fit quality to merit the addition of six additional free parameters. The moon model has a small increase in to 1.02. According to the Bayesian information criterion (BIC), which penalizes unnecessary model complexity, the moon model is disfavored with . This constitutes strong evidence against the inclusion of a moon (kass95). In addition, as shown in Figure 4.1, the posterior distribution for the moon transit time spans the entire duration of the observations ( confidence). The upper limit on the moon radius is at confidence.
We found that our results were unchanged when we used the exponential systematics model from TK18 (described in § 3.2). For this case, the moon model is also strongly disfavored (). The fit rms is within 5 ppm of the XY decorrelation model.
The posterior probabilites for the planet’s mass and radius are consistent with either a gas giant planet or brown dwarf. The radius is and mass is .
4.1 Comparison with Teachey & Kipping (2018)
TK18 found evidence for the transit of a Neptune-size moon in their analysis of the HST data, in contrast to the findings presented here. To compare our results with theirs, we fit the TK18 light curve directly. We fit the astrophysical signal with both the no-moon and models, and used the exponential systematics model. Figure 3.1 shows the best fit models.
Similar to the findings of TK18, the moon model improved the fit quality by a . Notably, however, the moon model fit to the TK18 data does not perform better than the no-moon model fit to our new data (rms of 362 versus 356 ppm), even with the additional seven free parameters.
The moon model also yields qualitatively different posterior distributions for the two data sets. As shown in Figure 4.1, for the TK18 data the moon radius and transit time are peaked at and . By contrast, the fit to the new data presented here yields an upper limit on the moon radius of at confidence, and the transit time is unconstrained.
Although the two data sets yield different constraints on the moon properties, the planet’s mid-transit time agrees to better than for the two fits. The transit time is earlier than expected based on the Kepler data ( confidence; TK18), suggesting that there are transit timing variations in the system. Such variation could arise from the presence of a moon, as suggested by TK18; however, the variation could also be caused by another planet in the system.
5 Discussion & Conclusions
A natural question arising from our analysis is the source of the discrepancy between TK18 and the new results presented here. We find that with our new data, there is strong evidence against the moon (), even when we use a comparable model to TK18 (a 6-parameter moon model and an exponential fit to the visit-long trend).
If the model is not the source of the difference, it must be the extracted transit light curves themselves. Figure 2.1 shows a direct comparison of the light curves. The count rate we measure is lower than the TK18 light curve, and there is a small bump in the difference between the two data sets near the location of the moon transit identified in the TK18 data (see the bottom panel). This bump may be the source of the moon feature reported in TK18.
We explored several modifications to our pipeline to attempt to reproduce the TK18 data reduction. These included rotating the image by 0.5 degrees, using the same aperture as TK18 to extract the spectrum, and scaling the master sky flat for the background subtraction (rather than just subtracting the median). None of these modifications had a significant effect on our results.
There are a few other steps in the TK18 data reduction that would require substantial modification of our pipeline to recreate, but seem unlikely to be responsible for the difference. One of these is outlier masking. TK18 identify outliers with a Gaussian process fit to the pixel-level light curves, compared to our optimal extraction approach. Despite the difference, both methods flag 0.01% of pixels as bad. We also do not use the STScI software aXeprep to embed the raw image in a larger array; however, this process primarily affects the edge of the image, many pixels distant from the extraction box, so it is unclear how this step would bias the light curve. We conclude that no single choice in the data reduction provides an easy explanation for the difference in our light curves.
During the referee process for this work, we learned of another manuscript that also reanalyzed the HST transit observation (heller19). The best fit favored a moon model similar to that found by TK18; however, an MCMC analysis did not converge on this model, leading the authors to conclude that the highest likelihood solution may be an outlier.
Taken together, these findings illustrate the challenge of pushing measurement precision to the 100 ppm level, and highlight the importance of developing multiple independent pipelines to confirm cutting-edge results.
We thank A. Teachey for helpful discussions and for providing the extracted light curve from TK18. We also thank the anonymous referee for a thoughtful report that improved the manuscript. The HST data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST). STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. Support for MAST for non-HST data is provided by the NASA Office of Space Science via grant NNX13AC07G and by other grants and contracts. We also use data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.