© ESO 2018

1 Introduction

Where are they? – With more than 180 moons discovered around the eight solar system planets and over 3500 planets confirmed beyond the solar system, an exomoon detection could be imminent. While many methods have indeed been proposed to search for moons around extrasolar planets (Sartoretti & Schneider 1999; Han & Han 2002; Cabrera & Schneider 2007; Moskovitz et al. 2009; Kipping 2009; Simon et al. 2010; Peters & Turner 2013; Heller 2014; Ben-Jaffel & Ballester 2014; Agol et al. 2015; Forgan 2017; Vanderburg et al. 2018)¹, only a few dedicated surveys have actually been carried out (Szabó et al. 2013; Kipping et al. 2013a,b, 2014; Hippke 2015; Kipping et al. 2015; Lecavelier des Etangs et al. 2017; Teachey et al. 2018), one of which is the Hunt for Exomoons with Kepler (HEK for short; Kipping et al. 2012).

In the latest report of the HEK team, Teachey et al. (2018) find evidence for an exomoon candidate around the roughly Jupiter-sized exoplanet candidate Kepler-1625 b, which they provisionally refer to as Kepler-1625 b–i. Kepler-1625 is a slightly evolved G-type star with a mass of $M_{\star }={1.079}_{-0.138}^{+0.100}{M}_{\odot }$ (M_⊙ being the solar mass), a radius of $R_{\star }={1.793}_{-0.488}^{+0.263}{R}_{\odot }$ (with R_⊙ as the solar radius), and an effective temperature of $T_{\text{eff},\star }={5548}_{-72}^{+83}$ K (Mathur et al. 2017). Its Kepler magnitude of 15.756 makes it a relatively dim Kepler target². The challenge of this tentative detection is in the noise properties of the data, which are affected by the systematic noise of the Kepler space telescope and by the astrophysical variability of the star. Although the exomoon signal did show up both around the ingress/egress regions of the phase-folded transits (known as the orbital sampling effect; Heller 2014; Heller et al. 2016a) generated by Teachey et al. (2018) and in the sequence of the three individual transits, it could easily have been produced by systematic errors or stellar variability, as pointed out in the discovery paper.

The noise properties also dictate a minimum size for an exomoon detected around a given star and with a given instrument. In the case of Kepler-1625, we calculate the root-mean-square of the noise level to be roughly 700 ppm. As a consequence, any moon would have to be at least about $\sqrt{700}\text{ppm} \times 1.793R_{\odot }\approx 5.2R_{\oplus }$ (R_⊕ being the Earth’s radius) in size, about 30% larger than Neptune, in order to significantly overcome the noise floor in a single transit. The three observed transits lower this threshold by a factor of $\sqrt{3}$, suggesting a minimum moon radiusof ≈3 R_⊕. In fact, the proposed moon candidate is as large as Neptune, making this system distinct from any planet–moon system known in the solar system (Heller 2018).

Here wepresent a detailed study of the three publicly available transits of Kepler-1625 b. Our aim is to test whether the planet-with-moon hypothesis is favored over the planet-only hypothesis.

First wedevelop a model to simulate photometric transits of a planet with a moon (see Sect. 2.2.2). Then we implement adetrending method following Teachey et al. (2018) and explore alternative detrending functions. Subsequently, we detrend the original Kepler-1625 light curve, determine the most likely moon parameters, and assess if the planet-with-moon hypothesis is favored over the planet-only hypothesis. Finally, we perform a blind injection-retrieval test. To preserve the noise properties of the Kepler-1625 light curve, we inject planet-with-moon and planet-only transits into out-of-transit parts of the Kepler-1625 light curve.

2 Methods

The main challenge in fitting a parameterized, noiseless model to observed data is the removal of noise on timescales similar or larger than the timescales of the effect to be searched for; at the same time, the structure of the effect must be untouched, an approach sometimes referred to as “pre-whitening” of the data (Aigrain & Irwin 2004). The aim of this approach is to remove unwanted variations in the data caused by, for example, stellar activity, systematic errors, or instrumental effects. This approach bears the risk of both removing actual signal from the data and of introducing new systematic variability. The discovery and refutal of the exoplanet interpretation of variability in the stellar radial velocities of α Centauri B serves as a warning example (Dumusque et al. 2012; Rajpaul et al. 2016). Recently developed Gaussian process frameworks, in which the systematics are modeled simultaneously with stellar variability, could be an alternative method (Gibson et al. 2012). This has become particularly important for the extended Kepler mission (K2) that is now working with degraded pointing accuracy (Aigrain et al. 2015).

That being said, Teachey et al. (2018) applied a pre-whitening technique to both the Simple Aperture Photometry (SAP) flux and the Pre-search Data Conditioning (PDCSAP) flux of Kepler-1625 to determine whether a planet-only or a planet–moon model is more likely to have caused the observed Kepler data. In the following, we develop a detrending and model fitting procedure that is based on the method applied by Teachey et al. (2018), and then we test alternative detrending methods.

During Kepler’s primary mission, the star Kepler-1625 has been monitored for 3.5 yr in total, and five transits could have been observed. This sequence of transits can be labeled as transits 1, 2, 3, 4, and 5. Due to gaps in the data, however, only three transits have been covered, which correspond to transits 2, 4, and 5 in this sequence. Figure 1 shows the actual data discussed here. The entire SAP (left) and PDCSAP (right) light curves are shown in the top panels, and close-up inspections of the observed transit 2, 4, and 5 are shown in the remaining panels. The time system used throughout the article is the Barycentric Kepler Julian Date (BKJD), unless marked as relative to a transit midpoint.

A key pitfall of any pre-whitening or detrending method is the unwanted removal of signal or injection of systematic noise, the latter potentially mimicking signal. In our case of an exomoon search, we know that the putative signal would be restricted to a time-window around the planetary mid-transit, which is compatible with the orbital Hill stability of the moon. This criterion defines a possible window length that we should exclude from our detrending procedures. For a planet of ten Jupiter masses in a 287 d orbit around a 1.1 M_⊙ star (as per Teachey et al. 2018), this window is about 3.25 days either side of the transit midpoint (see Appendix A).

Although this window length is astrophysically plausible to protect possible exomoon signals, many other choices are similarly plausible, but they result in significantly different detrendings. Figure 2 illustrates the effect on the detrended light curve if two different windows around the midpoint of the planetary transit (here transit 5) are excluded from the fitting. We chose a fourth-order polynomial detrending function and a 7.5 d (blue symbols) or a 4 d (orange symbols) region around the midpoint to be excluded from the detrending, mainly for illustrative purposes. In particular, with the latter choice, we produce a moon-like signal around the planetary transit similar to the moon signal that appears in transit 5 in Teachey et al. (2018). For the former choice, however, this signal does not appear in the detrended light curve.

Teachey et al. (2018) use the Cosine Filtering with Autocorrelation Minimization (CoFiAM) detrending algorithm to detrend both the SAP and PDCSAP flux around the three transits of Kepler-1625 b. CoFiAM fits a series of cosines to the light curve, excluding a specific region around the transit. CoFiAM preserves the signal of interest by using only cosines with a period longer than a given threshold and therefore avoids the injection of artificial signals with periods shorter than this threshold. Teachey et al. (2018) also test polynomial detrending functions but report that this removes the possible exomoon signal. We choose to reimplement the CoFiAM algorithm as our primary detrending algorithm so as to remain as close as possible in our analysis to the work in Teachey et al. (2018). In our injection-retrieval test, we also use polynomials of second, third, and fourth-order for detrending. While low-order polynomials cannot generally fit the light curve as well as the series of cosines, the risk of injecting artificial signals may be reduced.

2.1 Trigonometric detrending

We implement the CoFiAM detrending algorithm as per the descriptions given by Kipping et al. (2013b) and Teachey et al. (2018). In the following, we refer to this reimplementation as trigonometric detrending as opposed to polynomial functions that we test as well (see Sect. 2.4.4).

The light curves around each transit are detrended independently. For each transit, we start by using the entire SAP flux of the corresponding quarter. We use the SAP flux instead of the PDCSAP flux to reproduce the methodology of Teachey et al. (2018) as closely as possible. The authors argue that the use of SAP flux avoids the injection of additional signals into the light curve that might have the shape of a moon signal. First, we remove outliers using a running median with a window length of 20 h and a threshold of three times the local standard deviation with the same window length. In order to achieve a fast convergence of our detrending and transit fitting procedures, we initially estimate the transit midpoints and durations by eye and identify data anomalies: for example, gaps and jumps (e.g. the jump 2 d prior to transit 2 and the gap 4 d after transit 4, see Fig. 1).

Jumps in the light curve can have different origins. The jumps highlighted around transit 2 in Fig. 1 are caused by a reaction wheel zero crossing event; the jump 5 d after transit 4 is caused by a change in temperature after a break in the data collection. Following Teachey et al. (2018), who ignore data points beyond gaps and other anomalous events for detrending, we cut the light curve around any of the transits as soon as it encounters the first anomaly, leaving us with a light curve of a total duration D around each transit (see top left panel in Fig. 3). In Sect. 2.4.4, we investigate the effect of including data beyond gaps. The detrending is then applied in two passes, using the first pass to get accurate transit parameters. In particular, we determine the duration (t_T) between the start of the planetary transit ingress and the end of the transit egress (Seager & Mallén-Ornelas 2003) and the second pass to generate the detrended light curve.

First pass. Using the estimated transit midpoints and durations, we calculate the time window (t_c, see top left panel in Fig. 3) around a given transit midpoint to be cut from the detrending fit as $t_{\text{c}}=f_{t_{\text{c}}}{t}_{\text{T}}$, where the factor $f_{t_{\text{c}}}$, relating the time cut around the transit to the transit duration, is an input parameter for the detrending algorithm. Specifically, t_c denotes the total length of time around the transit excluded from the detrending. We fit the detrending function, $${\mathcal{G}}^k(t,\stackrel{\rightharpoonup }{a},\stackrel{\rightharpoonup }{b})=a_0+{\sum }_{l=1}^k{a}_l\cos \left(l\frac{2\pi }{2D}t\right)+b_l\sin \left(l\frac{2\pi }{2D}t\right),$$(1)

to the light curve (excluding the region t_c around the transit) by minimizing the χ² between the light curve and ${\mathcal{G}}^k(t,\stackrel{\rightharpoonup }{a},\stackrel{\rightharpoonup }{b})$, where $\stackrel{\rightharpoonup }{a} = (a_0,a_1,\mathrm{...},a_k)$ and $\stackrel{\rightharpoonup }{b} = (b_1,b_2,\mathrm{...},b_k)$ are the free model parameters to be fitted. The parameter k is a number between 1 and k_max = round(2D∕t_p), where $t_{\text{p}}=f_{t_{\text{p}}}{t}_{\text{T}}$ is the timescale below which we want to preserve possible signals. $f_{t_{\text{p}}}$ is an input parameter to the detrending algorithm. For each k, we divide the light curve by ${\mathcal{G}}^k$, giving us the detrended light curves F^k. We calculate the first-order autocorrelation according to the Durbin & Watson (1950) test statistic for each F^k (excluding again the region around the transit). For each transit, we select the F^k with the lowest autocorrelation $F_{\text{min}}^k$ and combine these $F_{\text{min}}^k$ around each transit into our detrended light curve F. We fit the planet-only transit model to the detrended light curve F and compute the updated transit midpoints and duration t_T.

Second pass.The second pass repeats the steps of the first pass, but using the updated transit midpoints and durations as input. The resulting detrended light curve F is then usedfor our model fits with the ultimate goal being to assess whether or not an exomoon is a likely interpretation of the light curve signatures. We estimate the noise around each transit by taking the variance of F, excluding the transit region.

Figure 3 shows the detrending function as well as the detrended light curve for $f_{t_{\text{c}}}=2.2$ and $f_{t_{\text{p}}}=4.4$, corresponding to t_c = 1.6 d and t_p = 3.1 d.

Fig. 1 Kepler light curve of Kepler-1625. Left panel: Simple Aperture Photometry (SAP) flux. Right panel: Pre-search Data Conditioning Simple Aperture Photometry (PDCSAP) flux. The top panels show theentire light curves, respectively. The second, third, and fourth rows illustrate zooms into transits 2, 4, and 5 of Kepler-1625 b, respectively. These transits were shifted to the panel center and ± 10 d of data are shown around the transit mid-points. Some examples of jumps and gaps in the light curve are shown. Time is given as aBarycentric Kepler Julian Date.

Fig. 2 Example of how the detrending procedure alone can produce an exomoon-like transit signal around a planetary transit. We use “transit 5” of Kepler-1625 b as an example. Top panel: gray dots indicate the Kepler PDCSAP flux. The lines show a fourth-order polynomial fit for which we exclude 7.5 d (blue) or 4 d (orange) of data around the mid-point (dashed parts), respectively. Middle panel: dots show the detrended light curve derived from the blue polynomial fit in the top panel. The blue line illustrates a planet-only transit model. Bottom panel: dots visualize the detrended light curve using the orange polynomial fit from the top panel. We note the additional moon-like transit feature caused by the overshooting of the orange polynomial in the top panel. The orange line shows a planet–moon transit model with moon parameters as in Table 1 (see Fig. 4 for transit dynamics). As an alternative interpretation, the blue detrending function filters out an actually existing moon signature while the orange detrending fit preserves the moon signal.

2.2 Transit model

We construct two transit models, one of which contains a planet only and one of which contains a planet with one moon. We denote the planet-only model as ${\mathcal{M}}_0$ (the index referring to the number of moons) and the planet–moon model as ${\mathcal{M}}_1$. We do not consider models with more than one moon.

2.2.1 Planet-only model

${\mathcal{M}}_0$ assumes a circular orbit of the planet around its star. Given the period of that orbit (P) and the ratio between stellar radius and the orbital semi-major axis (R_⋆∕a), the sky-projected apparent distance to the star center relative to the stellar radius can be calculated as $$z=\sqrt{{\left[\frac{a}{R_{\star }}\sin \left(\frac{2\pi (t-t_0)}{P}\right)\right]}^2+{\left[b\cos \left(\frac{2\pi (t-t_0)}{P}\right)\right]}^2},$$(2)

where b is the transit impact parameter and t₀ is the time of the transit midpoint. We use the python implementation of the Mandel & Agol (2002) analytic transit model by Ian Crossfield³ to calculate the transit light curve based on the planet-to-star radius ratio (r_p = R_p∕R_⋆) and based on a quadratic limb-darkening law parametrized by the limb-darkening parameters q₁ and q₂ as given in Kipping (2013). We call this model light curve with zero moons ${\mathcal{F}}_0$.

2.2.2 Planet–moon model

In our planet–moon model, we assume a circular orbit of the local planet–moon barycenter around the star with an orbital period P_B, a semimajor axis a_B, and a barycentric transit midpoint time t_0,B. The projected distance of the barycenter to the star center relative to the stellar radius is calculated in the same way as in Eq. (2). The planet and moon are assumed to be on circular orbits around their common center of mass with their relative distances to the barycenter determined by the ratio of their masses M_p and M_s to the total mass M_p + M_s. The individual orbits of both the planet and the moon are defined by the total distance between the two objects a_ps, the planet mass M_p, the moon mass M_s and by the time of the planet–moon conjunction t_0,s, that is, the time at which the moon is directly in front of the planet as seen from an observer on Earth.

This model is degenerate in terms of the sense of orbital motion of the moon. In other words, a given planet–moon transit light curve can be produced by both a prograde and a retrograde moon (Lewis & Fujii 2014; Heller & Albrecht 2014). We restrict ourselves to prograde moons. The planet mass is set to a nominal ten Jupiter masses, as suggested by Teachey et al. (2018) and in agreement with the estimates of Heller (2018). This constraint simplifies the interpretation of the results substantially since the moon parameters are then unaffected by the planetary parameters. The moon mass is assumed to be much smaller than that of the planet. In fact, for a roughly Neptune-mass moon around a planet of ten Jupiter masses, we expect a TTV amplitude of 3–4 min and a TDV amplitude of 6–7 min, roughly speaking. Therefore, we simplify our model and set M_s = 0, which means that a_ps is equal to the distance between the moon and the planet–moon barycenter, a_s. The moon is assumed to have a coplanar orbit around the planet and, thus, to have the same transit impact parameter as the planet.

With these assumptions, the projected distance of the planet center to the star center relative to the stellar radius z_p is equal to that of the barycenter z_B. The projected distance of the moon center to the star center relative to the stellar radius z_s is given by $$\begin{array}{ll}{z}_{\text{s}}^2=\hfill & {\left[\frac{a_{\text{B}}}{R_{\star }}\sin \left(\frac{2\pi (t-t_{0,\text{B}})}{P_{\text{B}}}\right)+\frac{a_{\text{ps}}}{R_{\star }}\sin \left(\frac{2\pi (t-t_{0,\text{s}})}{P_{\text{s}}}\right)\right]}^2\hfill \\ \hfill & +{\left[b\cos \left(\frac{2\pi (t-t_{0,\text{B}})}{P_{\text{B}}}\right)\right]}^2,\hfill \end{array}$$(3)

where P_s is the orbital period of the moon calculated from the fixed masses and a_ps.

We calculate the transit light curves of both bodies and combine them into the total model light curve, which we call ${\mathcal{F}}_1$. We use the limb-darkening parameter transformation from Kipping (2013). For computational efficiency, we do not consider planet–moon occultations. For the planet–moon system of interest, occultations would only occur during about half of the transits (assuming a random moon phase) even if the moon orbital plane were perfectly parallel to the line of sight. Such an occultation would take about 1.5 h and would only affect 5–10% of the total moon signal duration.

In Table1 we give an overview of our nominal parameterization of the planet–moon model. In Fig. 4 we show the orbital dynamics of the planetand moon during transits 2, 4, and 5 using the nominal parameters in Table 1. This nominal parameterization was chosen to generate a model light curve that is reasonably close to the preferred model light curve found in Teachey et al. (2018), but it does not represent our most likely model fit to the data.

Fig. 3 Left panel: Kepler SAP flux around the transits used for the trigonometric detrending, our reimplementation of the CoFiAM algorithm. The data points denoted by open circles around the transits are excluded from the detrending fit. The black line shows the resulting light curve trend without the transit. Right panel: detrended transit light curves as calculated by the trigonometric detrending.

2.2.3 Finding the posterior probability distribution

We use the Markov chain Monte Carlo (MCMC) implementation Emcee (Foreman-Mackey et al. 2013)to estimate the posterior probability distribution of the parameters for model ${\mathcal{M}}_{(i)}$ (${\mathcal{M}}_0$ or ${\mathcal{M}}_1$). For this purpose, we need to formulate the probability density of a light curve as well as the prior of the parameters.

All threetransits taken together, we have a total of N detrended flux measurements (see Sect. 2.1). Given a set of parameters $\stackrel{\rightharpoonup }{\theta }$, model ${\mathcal{M}}_i$ produces a model light curve ${\mathcal{F}}_i(t,\stackrel{\rightharpoonup }{\theta })$. We assume that the noise is uncorrelated (see Appendix B) and Gaussian with a standard deviation σ_j at time t_j. This simplifies the joint probability density to a product of the individual probabilities. The joint probability density function of the detrended flux F(t) is given by $$p(F|\stackrel{\rightharpoonup }{\theta },{\mathcal{M}}_i)={\displaystyle \prod _{j=1}^N\frac{1}{\sqrt{2\pi {\sigma }_j^2}}}\exp \left(-\frac{{\left(F(t_j)-{\mathcal{F}}_i(t_j,\stackrel{\rightharpoonup }{\theta })\right)}^2}{2{\sigma }_j^2}\right).$$(4)

The noise dispersion σ_j has a fixed value for each transit.

Table 2 shows the parameter ranges that we explore. A prior is placed on the stellar mass according to the mass of ${1.079}_{-0.138}^{+0.100}{M}_{\odot }$ determined by Mathur et al. (2017). The stellar mass for a given parameter set is determined from the system’s total mass using P_B and a_B and subtracting the fixed planet mass of ten Jupiter masses.

A total of 100 walkers are initiated with randomly chosen parameters close to the estimated transit parameters. For the sake of fast computational performance, the walkers are initially separated into groups of 16 for the planet-only model and 24 for the planet–moon model (twice the number of parameters plus 2, respectively), temporarily adding walkers to fill the last group. To transform the initially flat distribution of walkers into a distribution according to the likelihood function, the walkers have to go through a so-called burn-in phase, the resulting model fits of which are discarded. We chose a burn-in phase for the walkers of 500 steps in both groups. Afterwards, we discard the temporarily added walkers, merge the walkers back together, and perform a second burn-in phase of 2200 steps with a length determined by visual inspection. Finally, we initiate the main MCMC run with a total of 8000 steps. We run the MCMC code on the detrended light curve using both the planet-only and the planet–moon models.

Fig. 4 Left panel: example of a simulated planet–moon transit light curve for transits 2, 4, and 5 using the nominal parameterization given in Table 1. The relative flux is the difference to the out-of-transit model flux and is givenin parts per thousand (ppt). Right panel: visualization of the orbital configurations during transits 2 (left column), 4 (center column), and 5 (right column). Labels 1–5 in the light curves refer to configurations 1–5 (see labels along the vertical axis). An animated version of this figure is available online.

2.3 Model selection

We use the Bayesian information criterion (BIC) to evaluate how well a model describes the observationsin relation to the number of model parameters and data points. The BIC of a given model ${\mathcal{M}}_i$ with m_i parameters is defined by Schwarz (1978) as $$\text{BIC}({\mathcal{M}}_i|F)=m_i\ln N-2\ln p(F|{\stackrel{\rightharpoonup }{\theta }}_{\text{max}},{\mathcal{M}}_i),$$(5)

where ${\stackrel{\rightharpoonup }{\theta }}_{\text{max}}$ is the set of parameters that maximizes the probability density function $p(F|\stackrel{\rightharpoonup }{\theta },{\mathcal{M}}_i)$ for a given the light curve F and model ${\mathcal{M}}_i$.

The difference of the BICs between two models gives an indication as to which model is more likely. In particular, $\Delta \text{BIC}({\mathcal{M}}_1,{\mathcal{M}}_0) \equiv \text{BIC}({\mathcal{M}}_1) - \text{BIC}({\mathcal{M}}_0) < 0$ if model ${\mathcal{M}}_1$ is more likely. We consider ΔBIC  <  6 (or Δ BIC  >  6) as strong evidence in favor of (or against) model ${\mathcal{M}}_1$ (see, e.g., Kass & Raftery 1995).

The best-fitting set of parameters derived from our MCMC runs (${\stackrel{\rightharpoonup }{\theta }}_{\text{max}}$) is then used to calculate $\Delta \text{BIC}({\mathcal{M}}_1,{\mathcal{M}}_0)$. For our calculations, we only use those parts of the light curve around the transits that could potentially be affected by a moon (3.25 d on each side of the transits, determined by the Hill radius R_Hill and the orbital velocity of the planet–moon barycenter; see Appendix A).

2.4 Injection-retrieval test

In order to estimate the likelihood of an exomoon feature to be due to either a real moon or due to noise, we perform several injection-retrieval experiments. One of us (MH) injected two cases of transits into the out-of-transit parts of the original PDCSAP Kepler flux. In one case, a sequence of three planet-only transits (similar to the sequence of real transits 2, 4, and 5) was injected, where the planet was chosen to have a radius of 11 R_⊕. In another case, a sequence of three transits of a planet-with-moon system, with properties similar to the proposed Jupiter–Neptune system, was injected. Author KR then applied the Baysian framework described above in order to evaluate the planet-only versus the planet-with-moon hypotheses, and in order to characterize the planet and (if present) its moon. As an important trait of our experiment, the author KR did not know which of the light curves contained only a planet and which contained also a moon.

2.4.1 Transit injections into light curves

For the injection part, we use PyOSE (Heller et al. 2016a,b) to create synthetic planet and moon ensemble transits. This code numerically integrates the non-occulted areas of the stellar disk to calculate the instantaneous flux of the star, which makes it a computationally slow procedure. We therefore use the analytical model described Sect. 2.2 for the retrieval part. In our model, the moon’s orbit is defined by its eccentricity (e_s, fixed at 0), a_s, its orbital inclination with respect to the circumstellar orbit (i_s, fixed at 90°), the longitude of the ascending node, the argument of the periapsis, and the planetary impact parameter (b, fixed at zero). Due to the small TTV and TDV amplitudes compared to the 29.4 min exposure of the Kepler long cadence data, we neglect the planet’s motion around the planet/moon barycenter (although PyOSE can model this dynamical effect as well) and assume that the moon orbits the center of the planet.

Our numerical code creates a spherical limb-darkened star on a two-dimensional grid of floating-point values. The sky-projected shapes of both the planet and the moon are modeled as black circles. The spatial resolution of the simulation is chosen to be a few million pixels so that the resulting light curve has a numerical error of < 1 ppm, which is negligible compared tothe ≈700 ppm noise level of the Kepler light curve. The initial temporal resolution of our model is equivalent to 1000 steps per planetary transit duration, which we then downsample to the observed 29.4 min cadence. The creation of one such light curve of a planet with a moon takes about one minute on a modern desktop computer.

We create a set of 100 such transit simulations of the planet–moon ensemble, where the two bodies move consistently during and between transits. All orbits are modeled to be co-planar and mutual planet–moon occultations are also included. For each transit sequence, the initial orbital phase of the planet–moon system is chosen randomly.

With P_B = 287.378949 d and P_s = 2.20833 d, the moon advances by ≈ 0.13 in phase between each subsequent planetary transit (P_B∕P_s ≈ 130.13). During a planetary transit, the moon advances by ≈0.36 rad in phase (the planetary transit duration is 0.7869 ± 0.0084 d).

We also create a set of 100 such transits that only have a transiting planet without a moon. In these cases, the planetary radius was increased slightly to match the average transit depth of planet and moon.

2.4.2 Testing the model-selection algorithm on synthetic light curves with white noise only

As a first validation of our injection-retrieval experiment and our implementation of the Bayesian statistical framework, we generate a new set of white noise light curves to test only the model comparison part of our pipeline without any effects that could possibly arise from imperfect detrending. Any effects that we would see in our experiments with the real Kepler-1625 light curve but not in the synthetic light curves with noise only could then be attributed to the imperfect detrending of the time-correlated (red) noise.

The author MH generated 200 synthetic light curves with ten different levels of white noise, respectively, ranging from root mean squares of 250 ppm to 700 ppm in steps of 50 ppm. This results in a total of 2000 synthetic light curves.The method described in Sect. 2.4.1 was used to inject three transits of a planet only into 100 light curves per noise level and three transits of a planet with a moon into the remaining 100 light curves per noise level. The initial orbital phases were randomly chosen and are different from the ones used to generate thelight curves in Sect. 2.4.3. The author MH delivered these light curves to the author KR withoutrevealing their specific contents. The author KR then ran our model selection algorithm to find the Δ BIC for each of the 2000 systems. After the ΔBICs were found, MH revealed the planet-only or the planet–moon nature of each light curve.

Figure 5 shows the resulting ΔBICs for each of the 2000 light curves, separated into the planet-only (left panel) and planet–moon injected systems (right panel) and sorted by the respective white noise level (along the abscissa). Each vertical column contains 100 light curves, respectively. For a noise level of 250 ppm, for example, our algorithm finds no false positive moons in the planet-only data, that is, no system with a Δ BIC < −6, while 1 case remains ambiguous (− 6  <  ΔBIC  <  6) and the other 99 cases are correctly identified as containing no moons. In the case of an injected planet–moon system instead, the algorithm correctly retrieves the moon in 100% of the synthetic light curves, that is, Δ BIC < −6 for all systems.

More generally, for the simulated planet-only systems, the false positive rate is 0% throughout all noise levels. Occasionally, a system is flagged as ambiguous, but overall the algorithm consistently classifies planet-only systems correctly as having no moon. Referring to the injected planet–moon system (right panel), our false negative rate rises steadily with increasing noise level. In fact, it reaches parity with the true positive rate between about 650 and 700 ppm.

In Fig. 6, we present a_s and R_s for each of the maximum-likelihood fits shown in Fig. 5. Each panel in Fig. 6 refers to one white noise level, that is, to one column in Fig. 5 of either the planet-only or the planet–moon injected system. In the case of an injected planet only (upper panels), the most likely values of a_s are distributed almost randomly over the range of values that we explore. On the other hand, R_s is constrained to a small range from about 1.5 R_⊕ at 250 ppm to roughly 3 R_⊕ at 700 ppm with the standard variation naturally increasing with the noise level.

The lower part of Fig. 6 shows the outcome of our planet–moon injection-retrievals from the synthetic light curves with white noise only. The correct parameters are generally recovered at all noise levels. In fact, we either recover the moon with a similar radius and orbital separation as the injection values (symbolized by blue points) or we find the moon to have very different radius and orbit while also rejecting the hypothesis of its presence in the first place (symbolized by red points). The distribution of these false negatives in the a_s – R_s plane resembles the distribution of the true negatives in the corresponding no-moon cases. The ambiguous runs with a Δ BIC around zero still mostly recover the injected moon parameters. This is especially clear for the 700 ppm level, with 50% more ambiguous runs than true positives, where most of the runs still recover the injected parameters.

Fig. 5 Difference between the BIC of the planet–moon model and the no-moon model for 2000 artificial white noise light curves at different noise levels, injected with simulated transits. In the left panel (100 × 10 light curves), a planet and moon transit was injected; in the right panel (100 × 10 light curves), only the planet. Each light curve consists of three consecutive transits. Each column is sorted by the Δ BIC. The Δ BIC threshold, over which a planet–moon or planet-only system is clearly preferred, is ± 6 with the state of systems with a ΔBIC between those values considered to be ambiguous.

2.4.3 Transit injection into real out-of-transit data

We inject synthetic transits into the Kepler-1625 PDCSAP data prior to our own detrending (see Sect. 2.4.1). We use the PDCSAP flux instead of SAP flux because (1) it was easier for us to automate the anomaly detection and (2) PDCSAP flux has been cleaned of common systematics. Since the PDC pipeline removes many of the jumps in the data, we can focus on a single type of anomaly, that is, gaps. Gaps are relatively easy to detect in an automated way, removing the requirement of visual inspection of each light curve. For the injection, we select out-of-transit parts of the Kepler-1625 light curve that have at least 50 d of mostly uninterrupted data (25 d to both sides of the designated time of transit injection), but accept the presence of occasional gaps with durations of up to several days during the injection process.

The set of 200 synthetic light curves was provided by MH to KR for blind retrieval without any disclosure as to which of the sequences have a moon. The time of mid-transit was communicated with a precision of 0.1 days to avoid the requirement of a pre-stage transit search. This is justified because the original transits of Kepler-1625 b have already been detected, and the transits are visible by eye and do not necessarily need computer-based searches. We provide the 200 datasets to the community for reproducibility⁴ and encourage further blind retrievals.

Fig. 6 Distribution of the median likelihood Rs and as for all the runs for the different noise levels, with the runs injecting planet and moon in the top panel and runs injecting only a planet in the bottom panel. The Δ BIC of the planet–moon model compared to the no-moon model for all runs is indicated by the color. Generally, runs with a low Δ BIC (indicating the presence of a moon) also are in the vicinity of the injected parameter.

2.4.4 Detrending of the transit-injected light curves

The detrending procedure for our injection-retrieval experiment differs from the one used to detrend the original light curve around the Kepler-1625 b transits (see Sect. 2.1) in two respects.

First, we test the effect of the detrending function. In addition to the trigonometric function, we detrend the light curve by polynomials of second, third, and fourth order.

In addition, we test if the inclusion or exclusion of data beyond any gaps in the light curve affects the detrending. In one variation of our detrending procedure, we use the entire ± 25 d of data (excluding any data within t_c) around a transit midpoint. In another variation, we restrict the detrending to the data up to the nearest gap (if present) on both sides of the transit.

To avoid the requirement of time-consuming visual inspections of each light curve, we construct an automatic rule to determine the presence of gaps, which are the most disruptive kind of artifact to our detrending procedure. We define a gapas an interruption of the data of more than half a day. Whenever we do detect a gap, we cut another 12 h at boththe beginning and the end of the gap, since our visual inspection of the data showed that many gaps are preceded or followed by anomalous trends (see, e.g., the gap 4 d after transit 4 in Fig. 1).

We ignore any data points within t_c around the transit midpoint (see Sect. 2.1). If a gap starts within an interval [t_c ∕2, t_c∕2  +  12 h] around the transit midpoint, then we lift our constraint of dismissing a 12 h interval around gaps and use all the data within [t_c ∕2, t_c∕2  +  12 h] plus any data up to 12 h around the next gap.

If all these cuts result in no data points for the detrending procedure to one side of one of the three transits in a sequence, then we ignore the entire sequence for our injection-retrieval experiment. This is the case for 40 out of the 200 artificially injected light curves. This high loss rate of our experimental data is a natural outcome of the gap distribution in the original Kepler-1625 light curve. We exclude these 40 light curves for all variations of the detrending procedure that we investigate. All things combined, these constraints produce synthetic light curves with gap characteristics similar to the original Kepler-1625 b transits (see Fig. 3), that is, we allow the simulation of light curves with gaps close to but not ranging into the transits. The four detrending functions and our two ways of treating gaps yield a total of eight different detrending methods that we investigate (see Table 3).

3 Results

3.1 Analysis of the original Kepler-1625 b transits

Our first result is a reproduction of a detrended transit light curve of Kepler-1625 b that has the same morphology and moon characterization as the one proposed by Teachey et al. (2018) and that has a negative Δ BIC. We explore the variation of the free parameters of our trigonometric detrending procedure, $f_{t_{\text{p}}}$ and $f_{t_{\text{c}}}$, and identify such a detrended light curve for $f_{t_{\text{p}}}=4.4$ and $f_{t_{\text{c}}}=2.2$. Figure 7 shows the resulting light curve.

In Fig. 8, we show the results of our MCMC analysis of this particular light curve, which yields a moon with $a_{\text{s}}={16.3}_{-1.9}^{+5.0} R_{\text{J}}$ and $R_{\text{s}}={2.87}_{-0.94}^{+0.61} R_{\oplus }$. While both the moon radius and semimajor axis are well constrained, the distribution of the initial planet–moon orbital conjunction (t_0,s) fills out almost the entire allowed parameter range from − 1∕2 P_s to + 1∕2 P_s. The planetary radius is ${0.863}_{-0.051}^{+0.072} R_{\text{J}}$, the stellar radius is $R_{\star }={1.57}_{-0.09}^{+0.11} R_{\odot }$, and the density is ${\rho }_{\star }={0.26}_{-0.05}^{+0.04} {\rho }_{\odot }$. The point of maximum likelihood in the resulting MCMC distribution is at a_s = 14.7 R_J, R_s = 3.4 R_⊕, R_⋆ = 1.57 R_⊙, ρ_⋆ = 0.23 ρ_⊙, and R_p = 8.63 R_J. The $\Delta \text{BIC}({\mathcal{M}}_1,{\mathcal{M}}_0)$ we found is –4.954, indicating moderate evidence in favor of an exomoon being in the light curve.

3.2 Injection-retrieval experiment

In Fig. 9, we show the ΔBIC for the 160 simulated Kepler light curves that were not rejected by our detrending method due to gaps very close to a transit. The left panel shows our results for the analysis of planet-only injections and the right panel refers to planet–moon injections. The tables in the panel headers list the true negative, false positive, true positive, and false positive rates as well as the rates of ambiguous cases. With “positive” (“negative”), we here refer to the detection (non-detection) of a moon.

In particular, we find the true negative rate (left panel, Δ BIC ≥ 6) to be between 65% and 87.5% and the true positive rate (right panel, ΔBIC ≤ − 6) to be between 31.25% and 46.25% depending on the detrending method, respectively.

The rates of false classifications is between 8.75% and 17.5% for the injected planet-only systems with a falsely detected moon (false positives) and between 30% and 41.25% for the injected planet–moon systems with a failed moon recovery (false negatives).

The rates of classification as a planet–moon system depend significantly on the treatment of gaps during the detrending procedure. Whenever the light curve is cut at a gap, the detection rates for a moon increase – both for the false positives and for the true positives. Among all the detrending methods, this effect is especially strong for the trigonometric detrending. The false positive rate increases by almost a factor of two from 8.75% (T/N) to 16.25% (T/G) and the true positive rate increases by 15%–46.25%. The effect on the true negative rate is strongest for the trigonometric detrending, decreasing from 87.5% when the light curve is not cut at gaps (T/N) to 72.5% if the light curve is cut (T/G). The false negative rate for the second-order polynomial detrending decreases from 41.25% (P2/N) to 30% (P2/G) when gaps are cut, while the false negative rates of the other detrending methods remain almost unaffected.

Of all the light curves with an injected planet only, 21.25% have an ambiguous classification with at least one of the detrending methods showing a negative ΔBIC and a different method showing a positive ΔBIC above the threshold. For the light curves with an injected planet–moon system, there are 18.75% with ambiguous classification and another 18.75% of the injected planet–moon systems are classified unanimously as true positives by all detrending methods.

Figure 10 shows the distribution of the retrieved moon parameters a_s and R_s as well as the corresponding ΔBIC (see color scale) for each of the detrending methods.

For the light curves with an injected planet–moon system (lower set of panels), the maximum likelihood values of a_s and R_s of the true positives (blue) generally cluster around the injected parameters. In particular, we find that the moon turns out to be more likely (deeper-blue dots) when it is fitted to have a larger radius. The parameters of the false positives (blue dots in the upper set of panels) are more widely spread out, with moon radii ranging between 2 and 5 R_⊕ and the moon semimajor axes spread out through essentially the entire parameter range that we explored. The clustering of median a_s at around 100 R_J is an artifact of taking the median over a very unlocalized distribution along a_s. For the polynomial detrending methods, there are a certain number of what one could refer to as mischaracterized true positives. In these cases, the ΔBIC-based planet-only versus planet–moon classification is correct but the maximum likelihood values are very different from the injected ones.

The correctly identified planet-only systems show a similar distribution of a_s and R_s as in our experiment with white noise only and a 700 ppm amplitude (Fig. 5).

Most surprisingly, and potentially most worryingly, the false positives (blue dots in the upper set of panels in Fig. 10) cluster around the values of the moon parameters found by Teachey et al. (2018), in particular if the light curve is cut at the first gap.

Fig. 7 Theobserved second, fourth, and fifth transits of Kepler-1625 b. Black dots refer to our detrendedlight curve from the trigonometric detrending procedure, and orange curves are the model light curves generated using the 100 best fitting parameter sets of the MCMC run. The Δ BIC, calculated from the most likely parameters, is − 4.954.

Fig. 8 Posterior probability distribution of the moon parameters generated by the MCMC algorithmfor the light curve detrended by the trigonometric detrending. The black vertical lines show the median of the posterior distribution, the black horizontal lines indicate the 1σ range around the median. The red vertical lines show the point of maximum likelihood. The locations of the Galilean moons are included in the lower left panel for comparison.

4 Discussion

In this article, we compare several detrending methods of the light curve of Kepler-1625, some of which were used by Teachey et al. (2018) in their characterization of the exomoon candidate around Kepler-1625 b. However, we do not perform an exhaustive survey of all available detrending methods; for example, we leave out Gaussian processes (Aigrain et al. 2016).

We show that the sequential detrending and fitting procedure of transit light curves is prone to introducing features that can be misinterpreted as signal, in our case as an exomoon. This “pre-whitening” method of the data must therefore be used with caution. Our investigations of a polynomial-based fitting and of a trigonometric detrending procedure show that the resulting best-fit model depends strongly on the specific detrending function; for example, on the order of the polynomial or on the minimum timescale (or wavelength) of a cosine. This is crucial for any search of secondary effects in the transit light curves – caused by moons, rings, evaporating atmospheres and so on – and is in stark contrast to a claim by Aizawa et al. (2017), who stated that neither the choice of the detrending function nor the choice of the detrending window of the light curve would have a significant effect on the result. We find that this might be true at the level of visual inspection by eye but not at the level of 100 ppm or below. Part of the difference between our findings and those of Aizawa et al. (2017) could be in the different timescales we investigate. While they considered the effect of stellar flairs on timescales of less than a day, much less than the two-day transit duration of their specific target, our procedure operates on various timescales of up to several weeks. Moreover, we develop a dynamical moon model to fit multiple transits, whereas Aizawa et al. (2017) study only a single transit.

Since the actual presence and the putative orbital position of a hypothetical exomoon around Kepler-1625 b is unknown a priori, it is unclear how much of the light curve would need to be protected from (or neglected for) the pre-fit detrending process in order to avoid a detrending of a possible moon signal itself. In turn, we show that in the case of Kepler-1625, different choices for this protected timescale around the transit yield different confidences and different solutions for a planet–moon system. We find that the previously announced solution by Teachey et al. (2018) is only one of many possibilities with similar likelihoods (specifically, BIC). This suggests, but by no means proves, that all of these solutions could, in fact, be due to red noise artifacts (e.g. stellar or instrumental) rather than indicative of a moon signal.

Our finding of higher true positive rate compared to a false positive rate from injection-retrieval experiments could be interpreted as moderate evidence in favor of a genuine exomoon. This interpretation, however, depends on the number of transiting planets and planet candidates around stars with similar noise characteristics that were included in the Teachey et al. (2018) search. Broadly speaking, if more than a handful of similar targets are studied, the probability of at least one false positive detection becomes quite likely.

5 Conclusions

We investigated the detrending of the transit light curve of Kepler-1625 b with a method very similar to the one used by Teachey et al. (2018) and then applied a Bayesian framework with MCMC modeling to search for a moon. Our finding of a Δ BIC of − 4.954 favors the planet–moon over the planet-only hypothesis. Although significant, this tentative detection fails to cross the threshold of − 6, which we would consider strong evidence of a moon. Our ΔBIC value would certainly change if we were to include the additional data from the high-precision transit observations executed in October 2017 with the Hubble Space Telescope (Teachey et al. 2018) in our analysis. Moreover, by varying the free parameters of our detrending procedure, we also find completely different solutions for a planet–moon system, that is, different planet–moon orbital configurations during transits and different moon radii or planet–moon orbital semimajor axes.

As an extension to this validation of the previously published work, we performed 200 injection-retrieval experimentsinto the original out-of-transit parts of the Kepler light curve. We also extended the previous work by exploring different detrending methods, such as second-, third-, and fourth-order polynomials as well as trigonometric methods, and find false-positive rates between 8.75% and 16.25%, depending on the method. Surprisingly, we find that the moon radius and planet–moon distances of these false positives are very similar to the ones measured by Teachey et al. (2018). In other words, in 8.75%–16.25% of the light curves that contained an artificially injected planet only, we find a moon that is about as large as Neptune and orbits Kepler-1625 b at about 20 R_J.

To summarize, we find tentative statistical evidence for a moon in this particular Kepler light curve of Kepler-1625, but we also show that a significant fraction of similar light curves, which contained a planet only, would nevertheless indicate a moon with properties similar to the candidate Kepler-1625 b–i. Clearly, stellar and systematic red noise components are the ultimate barrier to an unambiguous exomoon detection around Kepler-1625 b and follow-up observations have the potential of solving this riddle based on the framework that we present.

Of all the detrending methods we investigated, the trigonometric method, which is very similar to the CoFiAM method of Teachey et al. (2018), can produce the highest true positive rate. At the same time, however, this method also ranks among those producing the highest false positive rates as well. To conclude, we recommend that any future exomoon candidate be detrended with as many different detrending methods as possible to evaluate the robustness of the classification.

Fig. 9 Difference between the BIC of the planet–moon model and the no-moon model using different detrending methods for 160 light curves, generated using the PDCSAP flux of Kepler-1625, injected with three simulated transits. In the left panel (80 light curves), a planet and moon transit was injected; in the right panel (80 light curves), only the planet. Each light curve consists of three consecutive transits. Each row of eight detrending methods uses the same light curve. The rows are sorted by their mean Δ BIC, with black lines indicating the ΔBIC = {−6, 0, 6} positions for the mean ΔBIC per row.

Fig. 10 Distribution of the median likelihood Rs and as for the transits injected into different parts of the Kepler-1625 light curve, using different detrending methods. The Δ BIC of the planet-only model compared to the planet–moon model is indicated by the symbol color. The values of the moon semimajor axis (abscissa) and radius (ordinate) suggested by Teachey et al. (2018) are indicated with thin, gray lines in each sub-panel.

Movie

Movie of Fig. 4. Access here

Acknowledgements

We thank James Kuszlewicz and Jesper Schou for useful discussions. This work was supported in part by the German Aerospace Center (DLR) under PLATO Data Center grant 50OL1701. This paper includes data collected by the Kepler mission. Funding for the Kepler mission is provided by the NASA Science Mission directorate. This work has made use of data provided by NASA and the Space Telescope Science Institute. K.R. is a member of the International Max Planck Research School for Solar System Science at the University of Göttingen. K.R. contributed to the analysis of the simulated light curves, to the interpretation of the results, and to the writing of the article.

Appendix A Effect of the window length on the Bayesian information criterion

Given the constraint of orbital stability, a moon can only possibly orbit its planet within the planet’s Hill sphere. Therefore, transits may only occur within a certain time interval around the midpoint of the planetary transit. This time t_Hill can be calculated as $$t_{\text{Hill}}=\frac{\eta R_{\text{Hill}}}{v_{\text{orbit}}}=\eta \frac{P}{2\pi }\sqrt{\frac{M_{\text{p}}}{3M_{\star }}},$$(A.1)

where v_orbit is the orbital velocity of the planet–moon system around the star, M_p and M_⋆ are the planet and star mass,

respectively, P is the orbital period of the planet–moon system, and R_Hill is the Hill radius of the planet. η is a factor between 0 and 1, which has been numerically determined for prograde moons (η ≈   0.5) and for retrograde moons (η ≈ 1), details depending on the orbital eccentricities (Domingos et al. 2006). We focus on prograde moon and choose η = 0.5. For a 10 M_J planet in a 287 d orbit around a 1.1 M_⊙ star, the Hill time is t_Hill = 3.25 d.

As shown in Fig. 2, the length of the light curve, which is neglected for the polynomial fit, has a strong effect on the resulting detrended light curve. Figure A.1 shows the effect that different cutout times t_c and detrending base lines D can have on whether a moon is detected or not.

Fig. A.1 Detrending for different cutout times tc and base length D, color coded by the resulting ΔBIC using a second- and fourth-order polynomial function. While some of the detrending models corresponding to a large negative ΔBIC are clearly the result of incorrect detrending, it is much less clear for many other detrending models.

Appendix B Autocorrelation of detrended light curves

The autocorrelations of the detrended light curves are shown in Fig. B.1. For all three transits, the autocorrelation is close to zero, except for the zero-lag component. This suggests that it is reasonable to model the noise covariance matrix as a diagonal matrix (see Sect. 2.2.3).

Fig. B.1 The autocorrelation of the difference between the detrended light curve and the best fitting model.

References

All Tables

All Figures

Fig. 1 Kepler light curve of Kepler-1625. Left panel: Simple Aperture Photometry (SAP) flux. Right panel: Pre-search Data Conditioning Simple Aperture Photometry (PDCSAP) flux. The top panels show theentire light curves, respectively. The second, third, and fourth rows illustrate zooms into transits 2, 4, and 5 of Kepler-1625 b, respectively. These transits were shifted to the panel center and ± 10 d of data are shown around the transit mid-points. Some examples of jumps and gaps in the light curve are shown. Time is given as aBarycentric Kepler Julian Date.

In the text

Fig. 2 Example of how the detrending procedure alone can produce an exomoon-like transit signal around a planetary transit. We use “transit 5” of Kepler-1625 b as an example. Top panel: gray dots indicate the Kepler PDCSAP flux. The lines show a fourth-order polynomial fit for which we exclude 7.5 d (blue) or 4 d (orange) of data around the mid-point (dashed parts), respectively. Middle panel: dots show the detrended light curve derived from the blue polynomial fit in the top panel. The blue line illustrates a planet-only transit model. Bottom panel: dots visualize the detrended light curve using the orange polynomial fit from the top panel. We note the additional moon-like transit feature caused by the overshooting of the orange polynomial in the top panel. The orange line shows a planet–moon transit model with moon parameters as in Table 1 (see Fig. 4 for transit dynamics). As an alternative interpretation, the blue detrending function filters out an actually existing moon signature while the orange detrending fit preserves the moon signal.

In the text

Fig. 3 Left panel: Kepler SAP flux around the transits used for the trigonometric detrending, our reimplementation of the CoFiAM algorithm. The data points denoted by open circles around the transits are excluded from the detrending fit. The black line shows the resulting light curve trend without the transit. Right panel: detrended transit light curves as calculated by the trigonometric detrending.

In the text

Fig. 4 Left panel: example of a simulated planet–moon transit light curve for transits 2, 4, and 5 using the nominal parameterization given in Table 1. The relative flux is the difference to the out-of-transit model flux and is givenin parts per thousand (ppt). Right panel: visualization of the orbital configurations during transits 2 (left column), 4 (center column), and 5 (right column). Labels 1–5 in the light curves refer to configurations 1–5 (see labels along the vertical axis). An animated version of this figure is available online.

In the text

Fig. 5 Difference between the BIC of the planet–moon model and the no-moon model for 2000 artificial white noise light curves at different noise levels, injected with simulated transits. In the left panel (100 × 10 light curves), a planet and moon transit was injected; in the right panel (100 × 10 light curves), only the planet. Each light curve consists of three consecutive transits. Each column is sorted by the Δ BIC. The Δ BIC threshold, over which a planet–moon or planet-only system is clearly preferred, is ± 6 with the state of systems with a ΔBIC between those values considered to be ambiguous.

In the text

Fig. 6 Distribution of the median likelihood Rs and as for all the runs for the different noise levels, with the runs injecting planet and moon in the top panel and runs injecting only a planet in the bottom panel. The Δ BIC of the planet–moon model compared to the no-moon model for all runs is indicated by the color. Generally, runs with a low Δ BIC (indicating the presence of a moon) also are in the vicinity of the injected parameter.

In the text

	Fig. 7 Theobserved second, fourth, and fifth transits of Kepler-1625 b. Black dots refer to our detrendedlight curve from the trigonometric detrending procedure, and orange curves are the model light curves generated using the 100 best fitting parameter sets of the MCMC run. The Δ BIC, calculated from the most likely parameters, is − 4.954.
In the text

Fig. 8 Posterior probability distribution of the moon parameters generated by the MCMC algorithmfor the light curve detrended by the trigonometric detrending. The black vertical lines show the median of the posterior distribution, the black horizontal lines indicate the 1σ range around the median. The red vertical lines show the point of maximum likelihood. The locations of the Galilean moons are included in the lower left panel for comparison.

In the text

Fig. 9 Difference between the BIC of the planet–moon model and the no-moon model using different detrending methods for 160 light curves, generated using the PDCSAP flux of Kepler-1625, injected with three simulated transits. In the left panel (80 light curves), a planet and moon transit was injected; in the right panel (80 light curves), only the planet. Each light curve consists of three consecutive transits. Each row of eight detrending methods uses the same light curve. The rows are sorted by their mean Δ BIC, with black lines indicating the ΔBIC = {−6, 0, 6} positions for the mean ΔBIC per row.

In the text

Fig. 10 Distribution of the median likelihood Rs and as for the transits injected into different parts of the Kepler-1625 light curve, using different detrending methods. The Δ BIC of the planet-only model compared to the planet–moon model is indicated by the symbol color. The values of the moon semimajor axis (abscissa) and radius (ordinate) suggested by Teachey et al. (2018) are indicated with thin, gray lines in each sub-panel.

In the text

	Fig. A.1 Detrending for different cutout times tc and base length D, color coded by the resulting ΔBIC using a second- and fourth-order polynomial function. While some of the detrending models corresponding to a large negative ΔBIC are clearly the result of incorrect detrending, it is much less clear for many other detrending models.
In the text

	Fig. B.1 The autocorrelation of the difference between the detrended light curve and the best fitting model.
In the text