Issue 
A&A
Volume 617, September 2018



Article Number  A49  
Number of page(s)  15  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201833085  
Published online  12 September 2018 
Revisiting the exomoon candidate signal around Kepler1625 b^{★}
^{1}
Institute for Astrophysics, Georg August University Göttingen,
FriedrichHundPlatz 1,
37077
Göttingen, Germany
^{2}
Max Planck Institute for Solar System Research,
JustusvonLiebigWeg 3,
37077
Göttingen, Germany
email: rodenbeck@mps.mpg.de
^{3}
Sonneberg Observatory,
Sternwartestr. 32,
96515
Sonneberg, Germany
Received:
23
March
2018
Accepted:
26
June
2018
Context. Transit photometry of the Jupitersized exoplanet candidate Kepler1625 b has recently been interpreted as showing hints of a moon. This exomoon, the first of its kind, would be as large as Neptune and unlike any moon we know from the solar system.
Aims. We aim to clarify whether the exomoonlike signal is indeed caused by a large object in orbit around Kepler1625 b, or whether it is caused by stellar or instrumental noise or by the data detrending procedure.
Methods. To prepare the transit data for model fitting, we explore several detrending procedures using second, third, and fourthorder polynomials and an implementation of the Cosine Filtering with Autocorrelation Minimization (CoFiAM). We then supply a light curve simulator with the coplanar orbital dynamics of the system and fit the resulting planet–moon transit light curves to the Kepler data. We employ the Bayesian information criterion (BIC) to assess whether a single planet or a planet–moon system is a more likely interpretation of the light curve variations. We carry out a blind hareandhounds exercise using many noise realizations by injecting simulated transits into different outoftransit parts of the original Kepler1625 light curve: (1) 100 sequences with three synthetic transits of a Kepler1625 blike Jupitersize planet and (2) 100 sequences with three synthetic transits of a Kepler1625 blike planet with a Neptunesized moon.
Results. The statistical significance and characteristics of the exomoonlike signal strongly depend on the detrending method (polynomials versus cosines), the data chosen for detrending, and the treatment of gaps in the light curve. Our injectionretrieval experiment shows evidence of moons in about 10% of those light curves that do not contain an injected moon. Strikingly, many of these falsepositive moons resemble the exomoon candidate, that is, a Neptunesized moon at about 20 Jupiter radii from the planet. We recover between about one third and one half of the injected moons, depending on the detrending method, with radii and orbital distances broadly corresponding to the injected values.
Conclusions. A ΔBIC of − 4.9 for the CoFiAMbased detrending is indicative of an exomoon in the three transits of Kepler1625 b. This solution, however, is only one out of many and we find very different solutions depending on the details of the detrending method. We find it concerning that the detrending is so clearly key to the exomoon interpretation of the available data of Kepler1625 b. Further highaccuracy transit observations may overcome the effects of red noise but the required amount of additional data might be large.
Key words: planets and satellites: detection / techniques: photometric / methods: data analysis
A movie associated to Fig. 4 is available at https://www.aanda.org.
© 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; BenJaffel & Ballester 2014; Agol et al. 2015; Forgan 2017; Vanderburg et al. 2018)^{1}, 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 Jupitersized exoplanet candidate Kepler1625 b, which they provisionally refer to as Kepler1625 b–i. Kepler1625 is a slightly evolved Gtype star with a mass of (M_{⊙} being the solar mass), a radius of (with R_{⊙} as the solar radius), and an effective temperature of K (Mathur et al. 2017). Its Kepler magnitude of 15.756 makes it a relatively dim Kepler target^{2}. 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 phasefolded 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 Kepler1625, we calculate the rootmeansquare of the noise level to be roughly 700 ppm. As a consequence, any moon would have to be at least about (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 , 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 Kepler1625 b. Our aim is to test whether the planetwithmoon hypothesis is favored over the planetonly 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 Kepler1625 light curve, determine the most likely moon parameters, and assess if the planetwithmoon hypothesis is favored over the planetonly hypothesis. Finally, we perform a blind injectionretrieval test. To preserve the noise properties of the Kepler1625 light curve, we inject planetwithmoon and planetonly transits into outoftransit parts of the Kepler1625 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 “prewhitening” 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 prewhitening technique to both the Simple Aperture Photometry (SAP) flux and the Presearch Data Conditioning (PDCSAP) flux of Kepler1625 to determine whether a planetonly 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 Kepler1625 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 closeup 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 prewhitening 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 timewindow around the planetary midtransit, 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 fourthorder 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 moonlike 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 Kepler1625 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 injectionretrieval test, we also use polynomials of second, third, and fourthorder for detrending. While loworder 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énOrnelas 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 , where the factor , 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, (1)
to the light curve (excluding the region t_{c} around the transit) by minimizing the χ^{2} between the light curve and , where and are the free model parameters to be fitted. The parameter k is a number between 1 and k_{max} = round(2D∕t_{p}), where is the timescale below which we want to preserve possible signals. is an input parameter to the detrending algorithm. For each k, we divide the light curve by , giving us the detrended light curves F^{k}. We calculate the firstorder 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 and combine these around each transit into our detrended light curve F. We fit the planetonly 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 and , corresponding to t_{c} = 1.6 d and t_{p} = 3.1 d.
Fig. 1 Kepler light curve of Kepler1625. Left panel: Simple Aperture Photometry (SAP) flux. Right panel: Presearch 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 Kepler1625 b, respectively. These transits were shifted to the panel center and ± 10 d of data are shown around the transit midpoints. Some examples of jumps and gaps in the light curve are shown. Time is given as aBarycentric Kepler Julian Date. 

Open with DEXTER 
Fig. 2 Example of how the detrending procedure alone can produce an exomoonlike transit signal around a planetary transit. We use “transit 5” of Kepler1625 b as an example. Top panel: gray dots indicate the Kepler PDCSAP flux. The lines show a fourthorder polynomial fit for which we exclude 7.5 d (blue) or 4 d (orange) of data around the midpoint (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 planetonly transit model. Bottom panel: dots visualize the detrended light curve using the orange polynomial fit from the top panel. We note the additional moonlike 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. 

Open with DEXTER 
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 planetonly model as (the index referring to the number of moons) and the planet–moon model as . We do not consider models with more than one moon.
2.2.1 Planetonly model
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 semimajor axis (R_{⋆}∕a), the skyprojected apparent distance to the star center relative to the stellar radius can be calculated as (2)
where b is the transit impact parameter and t_{0} is the time of the transit midpoint. We use the python implementation of the Mandel & Agol (2002) analytic transit model by Ian Crossfield^{3} to calculate the transit light curve based on the planettostar radius ratio (r_{p} = R_{p}∕R_{⋆}) and based on a quadratic limbdarkening law parametrized by the limbdarkening parameters q_{1} and q_{2} as given in Kipping (2013). We call this model light curve with zero moons .
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 Neptunemass 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 (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 . We use the limbdarkening 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. 

Open with DEXTER 
2.2.3 Finding the posterior probability distribution
We use the Markov chain Monte Carlo (MCMC) implementation Emcee (ForemanMackey et al. 2013)to estimate the posterior probability distribution of the parameters for model ( or ). 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 , model produces a model light curve . 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 (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 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 planetonly 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 socalled burnin phase, the resulting model fits of which are discarded. We chose a burnin 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 burnin 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 planetonly 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 outoftransit 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. 

Open with DEXTER 
Nominal parameterization of the planet–moon model to reproduce the transit shape suggested by Teachey et al. (2018).
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 with m_{i} parameters is defined by Schwarz (1978) as (5)
where is the set of parameters that maximizes the probability density function for a given the light curve F and model .
The difference of the BICs between two models gives an indication as to which model is more likely. In particular, if model is more likely. We consider ΔBIC < 6 (or Δ BIC > 6) as strong evidence in favor of (or against) model (see, e.g., Kass & Raftery 1995).
The bestfitting set of parameters derived from our MCMC runs () is then used to calculate . 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).
Parameter ranges explored with our planet–moon model.
2.4 Injectionretrieval 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 injectionretrieval experiments. One of us (MH) injected two cases of transits into the outoftransit parts of the original PDCSAP Kepler flux. In one case, a sequence of three planetonly 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 planetwithmoon 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 planetonly versus the planetwithmoon 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 nonocculted 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 limbdarkened star on a twodimensional grid of floatingpoint values. The skyprojected 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 coplanar 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 modelselection algorithm on synthetic light curves with white noise only
As a first validation of our injectionretrieval 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 Kepler1625 light curve but not in the synthetic light curves with noise only could then be attributed to the imperfect detrending of the timecorrelated (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 planetonly 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 planetonly (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 planetonly 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 planetonly systems, the false positive rate is 0% throughout all noise levels. Occasionally, a system is flagged as ambiguous, but overall the algorithm consistently classifies planetonly 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 maximumlikelihood 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 planetonly 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 injectionretrievals 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 nomoon 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 nomoon 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 planetonly system is clearly preferred, is ± 6 with the state of systems with a ΔBIC between those values considered to be ambiguous. 

Open with DEXTER 
2.4.3 Transit injection into real outoftransit data
We inject synthetic transits into the Kepler1625 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 outoftransit parts of the Kepler1625 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 midtransit was communicated with a precision of 0.1 days to avoid the requirement of a prestage transit search. This is justified because the original transits of Kepler1625 b have already been detected, and the transits are visible by eye and do not necessarily need computerbased searches. We provide the 200 datasets to the community for reproducibility^{4} and encourage further blind retrievals.
Fig. 6 Distribution of the median likelihood R_{s} and a_{s} 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 nomoon 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. 

Open with DEXTER 
2.4.4 Detrending of the transitinjected light curves
The detrending procedure for our injectionretrieval experiment differs from the one used to detrend the original light curve around the Kepler1625 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 timeconsuming 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 injectionretrieval 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 Kepler1625 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 Kepler1625 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).
Definition of the detrending identifiers in relation to the respective detrending functions that we explored in our transit injectionretrieval experiment of the Kepler1625 data.
3 Results
3.1 Analysis of the original Kepler1625 b transits
Our first result is a reproduction of a detrended transit light curve of Kepler1625 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, and , and identify such a detrended light curve for and . 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 and . 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 , the stellar radius is , and the density is . 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 we found is –4.954, indicating moderate evidence in favor of an exomoon being in the light curve.
3.2 Injectionretrieval 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 planetonly 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 (nondetection) 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 planetonly 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 secondorder 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 (deeperblue 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 ΔBICbased planetonly versus planet–moon classification is correct but the maximum likelihood values are very different from the injected ones.
The correctly identified planetonly 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 Kepler1625 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. 

Open with DEXTER 
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. 

Open with DEXTER 
4 Discussion
In this article, we compare several detrending methods of the light curve of Kepler1625, some of which were used by Teachey et al. (2018) in their characterization of the exomoon candidate around Kepler1625 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 “prewhitening” method of the data must therefore be used with caution. Our investigations of a polynomialbased fitting and of a trigonometric detrending procedure show that the resulting bestfit 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 twoday 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 Kepler1625 b is unknown a priori, it is unclear how much of the light curve would need to be protected from (or neglected for) the prefit detrending process in order to avoid a detrending of a possible moon signal itself. In turn, we show that in the case of Kepler1625, 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 injectionretrieval 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 Kepler1625 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 planetonly 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 highprecision 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 injectionretrieval experimentsinto the original outoftransit parts of the Kepler light curve. We also extended the previous work by exploring different detrending methods, such as second, third, and fourthorder polynomials as well as trigonometric methods, and find falsepositive 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 Kepler1625 b at about 20 R_{J}.
To summarize, we find tentative statistical evidence for a moon in this particular Kepler light curve of Kepler1625, 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 Kepler1625 b–i. Clearly, stellar and systematic red noise components are the ultimate barrier to an unambiguous exomoon detection around Kepler1625 b and followup 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 nomoon model using different detrending methods for 160 light curves, generated using the PDCSAP flux of Kepler1625, 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. 

Open with DEXTER 
Fig. 10 Distribution of the median likelihood R_{s} and a_{s} for the transits injected into different parts of the Kepler1625 light curve, using different detrending methods. The Δ BIC of the planetonly 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 subpanel. 

Open with DEXTER 
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 (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 t_{c} and base length D, color coded by the resulting ΔBIC using a second and fourthorder 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. 

Open with DEXTER 
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 zerolag 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. 

Open with DEXTER 
References
 Agol, E., Jansen, T., Lacy, B., Robinson, T. D., & Meadows, V. 2015, ApJ, 812, 5 [NASA ADS] [CrossRef] [Google Scholar]
 Aigrain, S., & Irwin, M. 2004, MNRAS, 350, 331 [NASA ADS] [CrossRef] [Google Scholar]
 Aigrain, S., Hodgkin, S. T., Irwin, M. J., Lewis, J. R., & Roberts, S. J. 2015, MNRAS, 447, 2880 [NASA ADS] [CrossRef] [Google Scholar]
 Aigrain, S., Parviainen, H., & Pope, B. J. S. 2016, MNRAS, 459, 2408 [NASA ADS] [Google Scholar]
 Aizawa, M., Uehara, S., Masuda, K., Kawahara, H., & Suto, Y. 2017, AJ, 153, 193 [NASA ADS] [CrossRef] [Google Scholar]
 BenJaffel, L., & Ballester, G. E. 2014, ApJ, 785, L30 [NASA ADS] [CrossRef] [Google Scholar]
 Cabrera, J., & Schneider, J. 2007, A&A, 464, 1133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Domingos, R. C., Winter, O. C., & Yokoyama, T. 2006, MNRAS, 373, 1227 [NASA ADS] [CrossRef] [Google Scholar]
 Dumusque, X., Pepe, F., Lovis, C., et al. 2012, Nature, 491, 207 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Durbin, J., & Watson, G. S. 1950, Biometrika, 37, 409 [Google Scholar]
 ForemanMackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [NASA ADS] [CrossRef] [Google Scholar]
 Forgan, D. H. 2017, MNRAS, 470, 416 [NASA ADS] [CrossRef] [Google Scholar]
 Gibson, N. P., Aigrain, S., Roberts, S., et al. 2012, MNRAS, 419, 2683 [NASA ADS] [CrossRef] [Google Scholar]
 Han, C., & Han, W. 2002, ApJ, 580, 490 [NASA ADS] [CrossRef] [Google Scholar]
 Heller, R. 2014, ApJ, 787, 14 [NASA ADS] [CrossRef] [Google Scholar]
 Heller, R. 2017, Detecting and Characterizing Exomoons and Exorings, eds. H. J. Deeg & J. A. Belmonte (Cham: Springer International Publishing), 1 [Google Scholar]
 Heller, R. 2018, A&A, 610, A39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Heller, R., & Albrecht, S. 2014, ApJ, 796, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Heller, R., Williams, D., Kipping, D., et al. 2014, Astrobiology, 14, 798 [NASA ADS] [CrossRef] [Google Scholar]
 Heller, R., Hippke, M., & Jackson, B. 2016a, ApJ, 820, 88 [NASA ADS] [CrossRef] [Google Scholar]
 Heller, R., Hippke, M., Placek, B., Angerhausen, D., & Agol, E. 2016b, A&A, 591, A67 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hippke, M. 2015, ApJ, 806, 51 [NASA ADS] [CrossRef] [Google Scholar]
 Hippke, M. 2018, Synthetic Dataset for the Exomoon Candidate Around Kepler1625 b, DOI: 10.5281/zenodo.1202034 [Google Scholar]
 Kass, R. E., & Raftery, A. E. 1995, J. Am. Stat. Assoc., 90, 773 [CrossRef] [MathSciNet] [Google Scholar]
 Kipping, D. M. 2009, MNRAS, 392, 181 [NASA ADS] [CrossRef] [Google Scholar]
 Kipping, D. M. 2013, MNRAS, 435, 2152 [NASA ADS] [CrossRef] [Google Scholar]
 Kipping, D. M., Bakos, G. Á., Buchhave, L., Nesvorný, D., & Schmitt, A. 2012, ApJ, 750, 115 [NASA ADS] [CrossRef] [Google Scholar]
 Kipping, D. M., Forgan, D., Hartman, J., et al. 2013a, ApJ, 777, 134 [NASA ADS] [CrossRef] [Google Scholar]
 Kipping, D. M., Hartman, J., Buchhave, L. A., et al. 2013b, ApJ, 770, 101 [NASA ADS] [CrossRef] [Google Scholar]
 Kipping, D. M., Nesvorný, D., Buchhave, L. A., et al. 2014, ApJ, 784, 28 [NASA ADS] [CrossRef] [Google Scholar]
 Kipping, D. M., Schmitt, A. R., Huang, X., et al. 2015, ApJ, 813, 14 [NASA ADS] [CrossRef] [Google Scholar]
 Lecavelier des Etangs, A., Hébrard, G., Blandin, S., et al. 2017, A&A, 603, A115 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lewis, K. M., & Fujii, Y. 2014, ApJ, 791, L26 [NASA ADS] [CrossRef] [Google Scholar]
 Mandel, K., & Agol, E. 2002, ApJ, 580, L171 [NASA ADS] [CrossRef] [Google Scholar]
 Mathur, S., Huber, D., Batalha, N. M., et al. 2017, ApJS, 229, 30 [NASA ADS] [CrossRef] [Google Scholar]
 Moskovitz, N. A., Gaidos, E., & Williams, D. M. 2009, Astrobiology, 9, 269 [NASA ADS] [CrossRef] [Google Scholar]
 Peters, M. A.,& Turner, E. L. 2013, ApJ, 769, 98 [NASA ADS] [CrossRef] [Google Scholar]
 Rajpaul, V., Aigrain, S., & Roberts, S. 2016, MNRAS, 456, L6 [NASA ADS] [CrossRef] [Google Scholar]
 Sartoretti, P., & Schneider, J. 1999, A&AS, 134, 553 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schwarz, G. 1978, Ann. Stat., 6, 461 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Seager, S., & MallénOrnelas, G. 2003, ApJ, 585, 1038 [NASA ADS] [CrossRef] [Google Scholar]
 Simon, A. E., Szabó, G. M., Szatmáry, K., & Kiss, L. L. 2010, MNRAS, 406, 2038 [NASA ADS] [Google Scholar]
 Szabó, R., Szabó, G., Dálya, G., et al. 2013, A&A, 553, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Teachey, A., Kipping, D. M., & Schmitt, A. R. 2018, AJ, 155, 36 [NASA ADS] [CrossRef] [Google Scholar]
 Vanderburg, A., Rappaport, S. A., & Mayo, A. W. 2018, AAS journals, submitted [arXiv:1805.01903] [Google Scholar]
For reviews see Heller et al. (2014) and Heller (2017).
NASA Exoplanet Archive: https://exoplanetarchive.ipac. caltech.edu.
Available at http://www.astro.ucla.edu/~ianc/files as python.py.
Available on Zenodo, [10.5281/zenodo.1202034], Hippke (2018).
Movie
Movie of Fig. 4. (Access here)
All Tables
Nominal parameterization of the planet–moon model to reproduce the transit shape suggested by Teachey et al. (2018).
Definition of the detrending identifiers in relation to the respective detrending functions that we explored in our transit injectionretrieval experiment of the Kepler1625 data.
All Figures
Fig. 1 Kepler light curve of Kepler1625. Left panel: Simple Aperture Photometry (SAP) flux. Right panel: Presearch 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 Kepler1625 b, respectively. These transits were shifted to the panel center and ± 10 d of data are shown around the transit midpoints. Some examples of jumps and gaps in the light curve are shown. Time is given as aBarycentric Kepler Julian Date. 

Open with DEXTER  
In the text 
Fig. 2 Example of how the detrending procedure alone can produce an exomoonlike transit signal around a planetary transit. We use “transit 5” of Kepler1625 b as an example. Top panel: gray dots indicate the Kepler PDCSAP flux. The lines show a fourthorder polynomial fit for which we exclude 7.5 d (blue) or 4 d (orange) of data around the midpoint (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 planetonly transit model. Bottom panel: dots visualize the detrended light curve using the orange polynomial fit from the top panel. We note the additional moonlike 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. 

Open with DEXTER  
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. 

Open with DEXTER  
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 outoftransit 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. 

Open with DEXTER  
In the text 
Fig. 5 Difference between the BIC of the planet–moon model and the nomoon 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 planetonly system is clearly preferred, is ± 6 with the state of systems with a ΔBIC between those values considered to be ambiguous. 

Open with DEXTER  
In the text 
Fig. 6 Distribution of the median likelihood R_{s} and a_{s} 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 nomoon 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. 

Open with DEXTER  
In the text 
Fig. 7 Theobserved second, fourth, and fifth transits of Kepler1625 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. 

Open with DEXTER  
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. 

Open with DEXTER  
In the text 
Fig. 9 Difference between the BIC of the planet–moon model and the nomoon model using different detrending methods for 160 light curves, generated using the PDCSAP flux of Kepler1625, 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. 

Open with DEXTER  
In the text 
Fig. 10 Distribution of the median likelihood R_{s} and a_{s} for the transits injected into different parts of the Kepler1625 light curve, using different detrending methods. The Δ BIC of the planetonly 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 subpanel. 

Open with DEXTER  
In the text 
Fig. A.1 Detrending for different cutout times t_{c} and base length D, color coded by the resulting ΔBIC using a second and fourthorder 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. 

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

Open with DEXTER  
In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.