An alternative interpretation of the exomoon candidate signal in the combined Kepler and Hubble data of Kepler-1625

Kepler and Hubble photometry of a total of four transits by the Jupiter-sized Kepler-1625b have recently been interpreted to show evidence of a Neptune-sized exomoon. The profound implications of this first possible exomoon detection and the physical oddity of the proposed moon, that is, its giant radius prompt us to re-examine the data and the Bayesian Information Criterion (BIC) used for detection. We combine the Kepler data with the previously published Hubble light curve. In an alternative approach, we perform a synchronous polynomial detrending and fitting of the Kepler data combined with our own extraction of the Hubble photometry. We generate five million MCMC realizations of the data with both a planet-only model and a planet-moon model and compute the BIC difference (DeltaBIC) between the most likely models, respectively. DeltaBIC values of -44.5 (using previously published Hubble data) and -31.0 (using our own detrending) yield strongly support the exomoon interpretation. Most of our orbital realizations, however, are very different from the best-fit solutions, suggesting that the likelihood function that best describes the data is non-Gaussian. We measure a 73.7min early arrival of Kepler-1625b for its Hubble transit at the 3 sigma level, possibly caused by a 1 day data gap near the first Kepler transit, stellar activity, or unknown systematics. The radial velocity amplitude of a possible unseen hot Jupiter causing Kepler-1625b's transit timing variation could be some 100m/s. Although we find a similar solution to the planet-moon model as previously proposed, careful consideration of its statistical evidence leads us to believe that this is not a secure exomoon detection. Unknown systematic errors in the Kepler/Hubble data make the DeltaBIC an unreliable metric for an exomoon search around Kepler-1625b, allowing for alternative interpretations of the signal.


Introduction
The recent discovery of an exomoon candidate around the transiting Jupiter-sized object Kepler-1625 b orbiting a slightly evolved solar mass star ) came as a surprise to the exoplanet community. This Neptune-sized exomoon, if confirmed, would be unlike any moon in the solar system, it would have an estimated mass that exceeds the total mass of all moons and rocky planets of the solar system combined. It is currently unclear how such a giant moon could have formed (Heller 2018). Rodenbeck et al. (2018) revisited the three transits obtained with the Kepler space telescope between 2009 and 2013 and found marginal statistical evidence for the proposed exomoon. Their transit injection-retrieval tests into the out-of-transit Kepler data of the host star also suggested that the exomoon could well be a false positive. A solution to the exomoon question was supposed to arrive with the new Hubble data of an October 2017 transit of Kepler-1625 b .
The new evidence for the large exomoon by , however, remains controversial. On the one hand, the Hubble transit light curve indeed shows a significant decrease in stellar brightness that can be attributed to the previously suggested moon. Perhaps more importantly, the transit of Kepler-1625 b occurred 77.8 min earlier than expected from a sequence of strictly periodic transits, which is in very good agreement with the proposed transit of the exomoon candidate, which occurred before the planetary transit. On the other hand, an upgrade of Kepler's Science Operations Center pipeline from version 9.0 to version 9.3 caused the exomoon signal that was presented in the Simple Aperture Photometry (SAP) measurements in the discovery paper  to essentially vanish in the SAP flux used in the new study of . This inconsistency, combined with the findings of Rodenbeck et al. (2018) that demonstrate that the characterization and statistical evidence for this exomoon candidate depend strongly on the methods used for data detrending, led us to revisit the exomoon interpretation in light of the new Hubble data.
Here we address two questions. How unique is the proposed orbital solution of the planet-moon system derived with the Bayesian information criterion (BIC)? What could be the reason for the observed 77.8 min difference in the planetary transit timing other than an exomoon?

Methods
Our first goal was to fit the combined Kepler and Hubble data with our planet-moon transit model (Rodenbeck et al. 2018) and to derive the statistical likelihood for the data to represent the model. In brief, we first model the orbital dynamics of the star-planet-moon system using a nested approach, in which the planet-moon orbit is Keplerian and unperturbed by the stellar gravity. The transit model consists of two black circles, one for the planet and one for the moon, that pass in front of the limb-darkened stellar disk. The resulting variations in the stellar brightness are computed using Ian Crossfield's python code of the Mandel & Agol (2002) analytic transit model 1 . The entire model contains 16 free parameters and it features three major updates compared to Rodenbeck et al. (2018): (1) planetmoon occultations are now correctly simulated, (2) the planet's motion around the local planet-moon barycenter is taken into account, and (3) inclinations between the circumstellar orbit of the planet-moon barycenter and the planet-moon orbit are now included.
We used the emcee code 2 of Foreman-Mackey et al. (2013) to generate Markov chain Monte Carlo (MCMC) realizations of our planet-only model (M 0 ) and planet-moon model (M 1 ) and to derive posterior probability distributions of the set of model parameters (θ). We tested both a standard MCMC sampling with 100 walkers and a parallel-tempering ensemble MCMC (PTMCMC) with five temperatures, each of which has 100 walkers. As we find a better convergence rate for the PTMCMC sampling, we use it in the following. Moreover, PTMCMC can sample both the parameter space at large and in regions with tight peaks of the likelihood function. The PTMCMC sampling is allowed to walk five million steps.
The resulting model light curves are referred to as F i (t, θ), where t are the time stamps of the data points from Kepler and Hubble (N measurements in total), for which time-uncorrelated standard deviations σ j at times t j are assumed, following the suggestion of . This simplifies the joint probability density of the observed (and detrended) flux measurements (F(t)) to the product of the individual probabilities for each data point, We then determined the set of parameters (θ max ) that maximizes the joint probability density function (p(F|θ max , M i )) for a given light curve F(t j ) and model M i and calculated the BIC (Schwarz 1978) The advantage of the BIC in comparison to χ 2 minimization, for example, is in its relation to the number of model parameters (m i ) and data points. The more free parameters in the model, the stronger the weight of the first penalty term in Eq.
(2), thereby mitigating the effects of overfitting. Details of the actual computer code implementation or transit simulations aside, this Bayesian framework is essentially what the Hunt for Exomoons with Kepler survey used to identify and rank exomoon candidates (Kipping et al. 2012), which ultimately led to the detection of the exomoon candidate around Kepler-1625 b after its first detection via the orbital sampling effect (Heller 2014;Heller et al. 2016a).

Data preparation
In a first step, we used Kepler's Pre-search Data Conditioning Simple Aperture Photometry (PDCSAP) and the Hubble Wide Field Camera 3 (WFC3) light curve as published by  based on their quadratic detrending. Then we executed our PTMCMC fitting and derived the ∆BIC values and the posterior parameter distributions. In a second step, we did our own extraction of the Hubble light curve including an exponential ramp correction for each Hubble orbit. Then we performed the systematic trend correction together with the transit fit of a planet-moon model. Our own detrending of the light curves is not a separate step, but it is integral to the fitting procedure. For each calculation of the likelihood, we find the best fitting detrending curve by dividing the observed light curve by the transit model and by fitting a thirdorder polynomial to the resulting light curve. Then we remove the trend from the original light curve by dividing it through the best-fit detrending polynomial and evaluate the likelihood. We also performed a test in which the detrending parameters were free PTMCMC model parameters and found similar results for the parameter distributions but at a much higher computational cost. We note that the resulting maximum likelihood is (and must be) the same by definition if the PTMCMC sampling converges.
Kepler-1625 was observed by Hubble under the GO program 15149 (PI Teachey). The observations were secured from October 28 to 29, 2017, to cover the ∼20 h transit plus several hours of out-of-transit stellar flux . The F130N filter of WFC3 was used to obtain a single direct image of the target, while 232 spectra were acquired with the G141 grism spanning a wavelength range from 1.1 to 1.7 µm. Due to the faintness of the target, it was observed in staring mode (e.g. Berta et al. 2012;Wilkins et al. 2014) unlike the most recent observations of brighter exoplanet host stars, which were monitored in spatial scanning mode (McCullough & MacKenty 2012). Hence, instead of using the IMA files as an intermediate product, we analyzed the FLT files, which are the final output of the calwfc3 pipeline of Hubble and allow a finer manipulation of the exposures during consecutive nondestructive reads. Each FLT file contains measurements between about 100 and 300 electrons per second, with exposure times of about 291 s.
We used the centroid of the stellar image to calculate the wavelength calibration, adopting the relations of Pirzkal et al. (2016). For each spectroscopic frame, we first rejected the pixels flagged by calwfc3 as "bad detector pixels", pixels with unstable response, and those with uncertain flux value (Data Quality condition 4, 32, or 512). Then we corrected each frame with the flat field file available on the Space Telescope Science Institute (STScI) website 3 by following the prescription of the WFC3 online manual. We performed the background subtraction on a column-by-column basis. Due to a number of contaminant stars in the observation field ( Fig. 1, top panel), we carefully selected a region on the detector that was as close as possible to the spectrum of Kepler-1625, close to row 150 in spatial dimension, and far from any contaminant. For each column on the detector, we applied a 5 σ clipping to reject the outliers and then calculated the median background flux value in that column. Following STScI prescriptions, we also removed pixels with an electron-per-second count larger than 5. An example for the background behavior is shown in the bottom panel of Fig. 1. We inspected each frame with the image registration package (Baker & Matthews 2001) to search drifts in both axes of the detector with respect to the very last frame, and then extracted the spectrum of Kepler-1625 by performing optimal extraction (Horne 1986) on the detector rows containing the stellar flux. This procedure automatically removes bad pixels and cosmic rays from the frames by correcting them with a smoothing function. We started the extraction with an aperture of a few pixels centered on the peak of the stellar trace and gradually increased its extension by one pixel per side on the spatial direction until the flux dispersion reached a minimum.
We performed another outlier rejection by stacking all the one-dimensional spectra along the time axis. We computed a median-filtered version of the stellar flux at each wavelength bin and performed a 3 σ clipping between the computed flux and the median filter. Finally, we summed the stellar flux across all wavelength bins from 1.115 to 1.645 µm to obtain the band-integrated stellar flux corresponding to each exposure.
Before performing the PTMCMC optimization, we removed the first Hubble orbit from the data set and the first data point of each Hubble orbit, as they are affected by stronger instrumental effects than the other observations (Deming et al. 2013) and cannot be corrected with the same systematics model. We also removed the last point of the 12th, 13th, and 14th Hubble orbit since they were affected by the passage of the South Atlantic Anomaly (as highlighted in the proposal file, available on the STScI website).

Mass-orbit constraints for a close-in planet
According to , the 2017 Hubble transit of Kepler-1625 b occurred about 77.8 min earlier than predicted, an effect that could be astrophysical in nature and is referred to as a transit timing variation (TTV). As proposed by , this TTV could either be interpreted as evidence for an exomoon or it could indicate the presence of a hitherto unseen additional planet. Various planetary configurations can cause the observed TTV effect such as an inner planet or an outer planet. At this point, no stellar radial velocity measurements of Kepler-1625 exist that could be used to search for additional nontransiting planets in this system.
In the following, we focus on the possibility of an inner planet with a much smaller orbital period than Kepler-1625 b simply because it would have interesting observational consequences. We use the approximation of Agol et al. (2005) for the TTV amplitude (δt) due to a close inner planet, which would impose a periodic variation on the position of the star, and solve their expression for the mass of the inner planet (M p,in ) as a function of its orbital semimajor axis (a p,in ), where a out = 0.87 AU is the semimajor axis of Kepler-1625 b. The validity of this expression is restricted to coplanar systems without significant planet-planet interaction and with a out a in , so that TTVs are only caused by the reflex motion of the star around its barycenter with the inner planet. As we show in Sect. 3.2, the proposed inner planet could be a hot Jupiter. The transits of a Jupiter-sized planet, however, would be visible in the Kepler data. As a consequence, we can estimate the minimum orbital inclination (i) between Kepler-1625 b and the suspected planet to prevent the latter from showing transits. This angle is given as per i = arctan(R /a p,in ) and we use R =  inner planet accordingly, and identify the region in the masssemimajor axis diagram of the inner planet that would lead to an overlap of the Hill spheres and therefore to orbital instability.

PTMCMC sampling and ∆BIC
Regarding the combined data set of the Kepler and Hubble data as detrended by , we find a ∆BIC of −44.5 between the most likely planet-only and the most likely planet-moon solution. A combination of the Kepler and Hubble light curves based on our own extraction of the WFC3 data yields a ∆BIC of −31.0. Formally speaking, both of these two values can be interpreted as strong statistical evidence for an exomoon interpretation. The two values are very different, however, which suggests that the detrending of the Hubble data has a significant effect on the exomoon interpretation. In other words, this illustrates that the systematics are not well-modeled and poorly understood. In Figs. 2a-d, we show our results for the PTMCMC fitting of our planet-moon model to the four transits of Kepler-1625 b including the Hubble data as extracted and detrended by Teachey & Kipping (2018) using a quadratic fit. Although our most likely solution shows some resemblance to the one proposed by , we find that several aspects are different. As an example, the second Kepler transit (Fig. 2b) is fitted best without a significant photometric moon signature, that is to say, the moon does not pass in front of the stellar disk 4 , whereas the corresponding best-fit model of  shows a clear dip prior to the planetary transit (see their Fig. 4). What is more, most of our orbital solutions (blue lines) differ substantially from the most likely solution (red line). In other words, the orbital solutions do not converge and various planet-moon orbital configuration are compatible with the data, though with lower likelihood.
In Figs. 2e-h, we illustrate our results for the PTM-CMC fitting of our planet-moon model to the four transits of Kepler-1625 b including our own extraction and detrending of the Hubble transit. Again, the orbital solutions (blue lines) do not converge. A comparison of panels d and h shows that the different extraction and detrending methods do have a significant effect on the individual flux measurements, in line with the findings of Rodenbeck et al. (2018). Although the time of the proposed exomoon transit is roughly the same in both panels, we find that the best-fit solution for the data detrended with our own reduction procedure does not contain the moon egress (panel h), whereas the best-fit solution of the data detrended by  does contain the moon egress (panel d). A similar fragility of this particular moon egress has been noted by  as they explored different detrending functions (see their Fig. 3).
Our Fig. 3 illustrates the distribution of the differential likelihood for the planet-moon model between the most likely model parameter set (θ max ) and the parameter sets (θ ) found after five million steps of our PTMCMC fitting procedure, p(θ |F, M 1 ) − p(θ max |F, M 1 ). For the combined Kepler and Hubble data detrended by Teachey & Kipping (2018) (left panel) and for our own Hubble data extraction and detrending (right panel), we find that most model solutions cluster around a differential likelihood that is very different from the most likely solution, suggesting that the most likely model is, in some sense, a statistical outlier. We initially detected this feature after approximately the first one hundred thousand PTMCMC fits. Hence, we increased the number of PTMCMC samplings to half a million and finally to five million to make sure that we sample any potentially narrow peaks of the likelihood function near the best-fit model at p(θ |F, M 1 ) − p(θ max |F, M 1 ) = 0 with sufficient accuracy. We find, however, that this behavior of the differential likelihood distribution clustering far from the best-fit solution persists, irrespective of the available computing power devoted to the sampling.  . Right: results from fitting our planet-moon transit model to our own detrending of the Kepler and WFC3 data. In both panels the most likely model is located at 0 along the abscissa by definition. In both cases the models do not converge to the best-fit solution, suggesting that the best-fit solution could in fact be an outlier.
our PTMCMC fitting of the combined Kepler and Hubble data (Hubble data as detrended and published by , and the bottom panel shows our PTMCMC fitting of the Kepler data combined with our own extraction and detrending of the Hubble light curve. The respective median values and standard deviations are noted in the upper right corners of each subpanel and summarized in Table 1. A comparison between the upper and lower corner plots in Fig. 4 reveals that the different detrending and fitting techniques have a significant effect on the resulting posterior distributions, in particular for i s and Ω s , the two angles that parameterize the orientation of the moon orbit. At the same time, however, the most likely values (red dots above the plot diagonal) and median values (blue crosses below the plot diagonal) of the seven parameters shown are well within the 1 σ tolerance.
The following features can be observed in both panels of Fig. 4. The moon-to-star radius ratio (Col. 1, leftmost) shows an approximately normal distribution, whereas the scaled planetmoon orbital semimajor axis (Col. 2) shows a more complicated, skewed distribution. The solutions for the orbital period of the exomoon candidate (Col. 3) show a comb-like structure owing to the discrete number of completed moon orbits that would fit a given value of the moon's initial orbital phase (Col. 4), which is essentially unconstrained. The moon-to-planet mass ratio (Col. 5) then shows a skewed normal distribution with a tail of large moon masses. Our results for the inclination i s between the satellite orbit (around the planet) and the line of sight, and for the longitude of the ascending node of the moon orbit are shown in Cols. 6 and 7. The preference of i s being either near 0 or near ±π (the latter is equivalent to a near-coplanar retrograde moon orbit) illustrates the well-known degeneracy of the prograde/retrograde solutions available from light curve analyses (Lewis & Fujii 2014;Heller & Albrecht 2014). Our planet-only model for the 2017 Hubble transit gives a transit midpoint at 3222.5547 ± 0.0014 d, which is 73.728(±2.016) min earlier than the predicted transit midpoint. This is in agreement with the measurements of , who found that the Hubble transit occurred 77.8 min earlier than predicted. This observed early transit of Kepler-1625 b has a formal ∼3 σ significance. We note, however, that this 3 σ deviation is mostly dictated by the first transit observed with Kepler (see Fig. S12 in . We also note that this transit was preceded by a ∼ 1 d observational gap in the light curve, about 0.5 d prior to the transit, which might affect the local detrending of the data and the determination of the transit mid-point of a planet-only model. Moreover, with most of the TTV effect being due to the large deviation from the linear ephemeris of the first transit, stellar (or any other systematic) variability could have a large (but unknown) effect on the error bars that go into the calculations.

Transit timing variations
In Fig. 5 we show the mass of an unseen inner planet that is required to cause the observed 73.7 min TTV amplitude from our PTMCMC fit as a function of its unknown orbital semimajor axis. The mass drops from 5.8 M Jup at 0.03 AU to 1.8 M Jup at 0.1 AU. Values beyond 0.1 AU cannot be assumed to fulfill the approximations made for Eq. (3) and are therefore shown with a dashed line. The actual TTV amplitude of Kepler-1625 b could even be higher than the ∼73 min that we determined for the Hubble transit, and thus the mass estimates shown for a possible unknown inner planet serve as lower boundaries.
The resulting radial velocity amplitudes of the star of 923 m s −1 (at 0.03 AU) and 151 m s −1 (at 0.1 AU), respectively, are indicated along the curve. Even if the approximations for a coplanar, close-in planet were not entirely fulfilled, our results suggest that RV observations of Kepler-1625 with a high-resolution spectrograph attached to a very large (8 m class) ground-based telescope could potentially reveal an unseen planet causing the observed TTV of Kepler-1625 b. Also shown along the curve in  . Bottom: results for our own detrending of the Kepler and WFC3 data. In both figures, scatter plots are shown with black dots above the diagonal, and projected histograms are shown as colored pixels below the diagonal. The most likely parameters are denoted with an orange point in the scatter plots. Histograms of the moon-to-star radius ratio r s , scaled semimajor axis of the planet-moon system (a ps /R ), satellite orbital period (P s ), satellite orbital phase (ϕ), moonto-planet mass ratio ( f M ), orbital inclination of the satellite with respect to our line of sight (i s ), and the orientation of the ascending node of the satellite orbit (Ω s ) are shown on the diagonal. Median values and standard deviations are indicated with error bars in the histograms.
A95, page 6 of 8 The pale red shaded region is excluded from a dynamical point of view since this is where the planetary Hill spheres would overlap. The extent of this region is a conservative estimate because it assumes a mass of 1 M Jup for Kepler-1625 b and neglects any chaotic effects induced by additional planets in the system or planet-planet cross tides etc. The true range of unstable orbits is probably larger. The black dots show all available exoplanet masses and semimajor axes from the Exoplanet Encyclopaedia, which illustrates that the suspected planet could be more massive than most of the known hot Jupiters. not be able to deliver the signal-to-noise ratios (S/N) required to validate or reject the exomoon hypothesis. With a Gaia G-band magnitude of m G = 15.76 (Gaia Collaboration 2016, 2018) the star is rather faint in the visible regime of the electromagnetic spectrum and the possible moon transits are therefore beyond the sensitivity limits of the TESS, CHEOPS, and PLATO missions. 2MASS observations suggest that Kepler-1625 is somewhat brighter in the near-infrared (Cutri et al. 2003), such that the James Webb Space Telescope (launch currently scheduled for early 2021) should be able to detect the transit of the proposed Neptune-sized moon, for example via photometric time series obtained with the NIRCam imaging instrument.
All things combined, the fragility of the proposed photometric exomoon signature with respect to the detrending methods, the unknown systematics in both the Kepler and the Hubble data, the absence of a proper assessment of the stellar variability of Kepler-1625, the faintness of the star (and the resulting photometric noise floor), the previously stated coincidence of the proposed moon's properties with those of false positives (Rodenbeck et al. 2018), the existence of at least one plausible alternative explanation for the observed TTV effect of Kepler-1625 b, and the serious doubts that we have about the ∆BIC as a reliable metric at least for this particular data set lead us to conclude that the proposed moon around Kepler-1625 b might not be real. We find that the exomoon hypothesis heavily relies on a chain of delicate assumptions, all of which need to be further investigated.
A similar point was raised by , and our analysis is an independent attempt to shed some light on the "unknown unknowns" referred to by the authors. For the time being, we take the position that the first exomoon has yet to be detected as the likelihood of an exomoon around Kepler-1625 b cannot be assessed with the methods used and data currently available.