© ESO 2020

1 Introduction

Among the many unexpected discoveries of the Kepler mission such as evaporating planets (Rappaport et al. 2012), exocomets (Rappaport et al. 2018), and even evidence of an exo-Trojan population (Hippke & Angerhausen 2015), exomoons have hitherto escaped unambiguous detection. The peculiar case of the exomoon candidate around Kepler-1625 b (Teachey et al. 2018; Teachey & Kipping 2018) is extremely hard to assess based on the very limited amount of data (Heller 2018a; Rodenbeck et al. 2018; Kreidberg et al. 2019).

Many ways have been proposed to look for moons in transit light curves (for a review see Heller 2018b). Sartoretti & Schneider (1999) suggested that a massive moon can distract the transits of its host planet from strict periodicity and induce transit timing variations (TTVs). Kipping (2009) identified an additional effect of the moon on the planet’s sky-projected velocity component tangential to the line of sight that results in transit duration variations (TDVs). These effects relate to the planet’s position with respect the planet–moon barycenter, which is modeled on a Keplerian orbit around the star, and they depend on the planet–satellite orbital semi-major axis (a_ps) and on the mass of the satellite (M_s) but on its radius (R_s). Heller et al. (2016) proposed that the combined planetary TTV and TDV effects can produce distinct ellipsoidal patterns in the TTV-TDV diagram if the moon is sufficiently small to avoid any effects on the photometric transit center.

However, Szabó et al. (2006) noted that the stellar dimming caused by a moon with a sufficiently large radius can affect the photometric center of the transit light curve, defined as τ = ∑_it_iΔm_i∕∑_iΔm_i with t_i being the times when the stellar differential magnitudes (Δm_i) are measured. For small (arbitrarily massive) moons, τ would coincide with the transit midpoint of the planet but for large moons τ can be substantially offset from the position of the midpoint of the planetary transit. Simon et al. (2007) derived an analytical estimate of the photometric TTV effect (dubbed “TTV_p”) and compared it to the magnitude of the barycentric TTV effect (dubbed “TTV_b”). These latter authors define the TTV_p as the difference between TTV_b and τ.

Here we simulate many different transit light curves of a planet with a moon of non-negligible radius and then fit the resulting data with a planet-only transit model to obtain the TTVs and TDVs of a hypothetical sequence of transits. Our aim is to provide the exoplanet community with a tool to inspect their TTV–TDV distributions for possible exomoon candidates without the need for developing a fully consistent, photodynamical transit modeling of a planet with a moon (Kipping 2011; Rodenbeck et al. 2018). We show that sufficiently large moons can distort the previously proposed elliptical pattern in the TTV–TDV distribution into very complicated patterns, which are much harder to discriminate from the TTV–TDV distribution of a planet without a moon. We also identify a new effect that appears in the sequence of transit radius (R_p) measurements and could be used to identify exomoon candidates in large exoplanet surveys.

2 Method

2.1 Setup of the transit model

Figure 1 demonstrates how the transit light curve of an exoplanet with a moon (solid black line) differs ever so slightly from the light curve expected for a single planet (orange dashed line). Firstly, the midpoint of the planetary transit changes between transits with respect to the planet–moon barycenter, thereby causing a TTV_b effect (Sartoretti & Schneider 1999). Moreover, the additional dimming of the star caused by a sufficiently large moon could act to shift the photometric center of the combined transit, an effect referred to as TTV_p (Szabó et al. 2006; Simon et al. 2007). The same principle of the moon’s photometric (rather than the barycentric) effect should be applicable to the TDV, though it has not been explored in the literature so far. Finally, Fig. 1 illustrates how the planet-only fit to the planet–moon model overestimates the radius of the planet. For single planets, the transit depth is roughly proportional to $R_{\text{p}}^2$, although details depend on the stellar limb darkening (Heller 2019). But a moon can affect the radius estimate for the planet and therefore mimic transit radius variations (TRVs) in a sequence of transit measurements.

We investigate both the barycentric and the photometric contributions to the measured TTVs and TDVs of a planet-with-moon transit light curve that is fitted with a planet-only model. To separate the TTV_b and TTV_p effects, we simulate transits of a hypothetical planet–moon system with a zero-radius or zero-mass moon, respectively, and then fit the simulated light curves with a planet-only transit model. Finally, we assume a realistic moon with both mass and radius and redo our fit to the resulting light curve.

We simulate exoplanet–exomoon transits using our implementation of the Mandel & Agol (2002) analytic transit model similar to Rodenbeck et al. (2018). However, as an improvement to Rodenbeck et al. (2018), our model now includes planet–moon occultations and the possibility of inclining the planet–moon orbit with respect to the line of sight. The planet-only transit model is fitted to the simulated light curve using the planet-to-star radius ratio (r_p), the transit duration (t_T), and the orbital phase (φ_b), while the orbital period (P_b) and the limb darkening coefficients (q₁ and q₂ for the quadratic limb darkening model as per Kipping 2013) are kept constant at the respective values shown in Table A.1. Figure 1 shows an example of a fitted planet-only light curve (orange dashed line) to the simulated noiseless transit of a planet–moon system.

We also investigate the possibility of detecting TRVs. Variations of the planetary transit depth have been observed before (Holczer et al. 2016) but to our knowledge they have not been treated in the context of exomoons so far.

2.2 Simulated noise in the light curves

We compare the effect of photon noise between two space missions: Kepler (Borucki et al. 2010) and PLATO (Rauer et al. 2014), the latter of which is currently scheduled for launch in 2026. For both missions we test stars of two different apparent magnitudes, namely 8 and 11. The two stellar magnitudes provide us with two possible amplitudes of the white noise level per data point, which is composed of photon noise and telescope noise. For Kepler long cadence observations with a cadence of 29.4 min¹, we obtain a white noise level of 9 parts per million (ppm) at magnitude 8 and 36 ppm at magnitude 11 (Koch et al. 2010). For PLATO we estimate 113 ppm at magnitude 8 and 448 ppm at magnitude 11, both at a cadence of 25 s. These estimates for PLATO assume an optimal target coverage by all 24 normal cameras.

We also simulate a time-correlated granulation noise component according to Gilliland et al. (2011), where the granulation power spectrum is modeled as a Lorentzian distribution in the frequency domain and then transformedinto the time domain. The granulation amplitude and timescale depend on the surface gravity and temperatureof the star. We use two stellar types, a sun-like star and a red dwarf with similar properties as ɛ Eridani (ɛ Eri).

Fig. 1 Example of a planet-only fit (dashed orange line) to the planet–moon model (solid black line). The star is sun-like, the planet is Jupiter-sized and in a 30 d orbit around the star, and the moon is Neptune-sized and in a 3.55 d orbit around the planet. In this example, the fitted planet-to-star radius ratio (rfit), which goes into the calculation of the transit depth with a power of two, is slightly overestimated due to the photometric signature of the moon. As a comparison between the amplitude of the moon signal and the expected noise, the noise from the Poisson statistics of an mV = 8 star (using a 1 h binning) is 12.6 ppm for PLATO and 6.0 ppm for Kepler. This is two magnitudes smaller than the amplitude of the moon signature but possibly comparable to the effect of stellar noise, from activity for example.

2.3 Parameterization of the test cases

We investigate several hypothetical star–planet–moon systems to explore the observability of the resulting TTV, TDV, and TRV effects. We ensure orbital stability of the moons by arranging them sufficiently deep into the gravitational potential of their host planet, that is, at < 0.5 times the planetary Hill radius. This distance has been shown to guarantee orbital stability of both pro- and retrograde moons based on numerical N-body simulations (Domingos et al. 2006). We choose stellar limb darkening coefficients (q₁, q₂) for both the sun and a K2V star akin to ɛ Eri from Claret & Bloemen (2011).

In Table A.1 we summarize our test cases. For case 1, and for all of its sub-cases, 1a–1p, we choose the planet to be Jupiter-sized and the moon to be Earth-sized. The orbit period of the moon can take two values, either the orbital period of Europa (3.55 d) or that of Io (1.77 d). The barycenter of the planet–moon system is in a P_b = 30 d orbit around its Sun-like star. In test cases 2a–2d, we simulate Neptune-sized moons and vary both orbitalperiods involved to test the effect of the number and duration of the transits. We set P_b to be either 30 d orbit or close to twice that value. However, we shorten the 60 d orbit by 1.69 d to 58.31 d for the Io-like orbit and by 1.59 d to 58.41 d for the Europa-wide orbit in order for the moon to show the same orbital advancement between transits, that is, for the remainder of P_b∕P_s to be the same (Heller et al. 2016). Cases 3a–3d refer to a Saturn-sized planet and a super-Earth moon, cases 4a–4d to a Neptune-sized planet and an Earth-sized moon, and cases 5a–5d to an Earth-Moon analog (though at either a 30 d or a 58.31 d orbital period around their star). In cases 6–9, we redo all these cases except that the star resembles ɛ Eri.

In cases 1a–1p, we study the effect of the moon’s orbital phase on the resulting transit shape and generate transits of planet–moon systems for orbital moon phases ranging from 0 to 1. We allow for different orbital inclinations of the moon orbit (i_s) in these test cases (see Fig. 2), which may prevent occultations if the planet–moon orbit is sufficiently wide for a given inclination. As we show below, the inclination has an effect on the measured transit parameters if the line connecting the planet and the moon is parallel to the line of sight, that is, if the planet and the moon have a conjunction during the transit. For test cases 2a–9d we choose the lowest inclination possible without causing a planet–moon occultation.

Fig. 2 Illustration of the orbital geometry and of our parameterization of the simulated planet–moon systems. φs is the orbital phase relative to the line of sight half-way through the transit.

2.4 Application to actual Kepler data

As a first application of our search for exomoon indicators, we explore the time series of the fitted transit depth from Holczer et al. (2016). We then do a “by-eye” vetting of the transit depth time series, of the autocorrelation functions (ACF), and of the periodograms for each of the 2598 Kepler Objects of Interest (KOIs) listed in Holczer et al. (2016) and identify KOI-1457.01 (Kepler-856 b) as an interesting candidate. Briefly, this is a roughly Saturn-sized validated planet with a false-positive probability of 2.9 × 10⁻⁵ (Morton et al. 2016) and an orbital period of about 8 d around an ɛ Eri-like host star.

3 Results

3.1 Signatures of exomoons in exoplanet-only transit model fits

Figure 3 shows the TDV (panel a), TTV (panel b), TRV (panel c), and combined TTV-TDV (panel d) effects. Each of these four panels is divided into three subplots showing the contributions of the photometric distortion of the light curve due to the moon (top subplots, moon mass set to zero), of the barycentric motion of the planet due to its moon (center subplots, moon radius set to zero), and the combined effect as derived from the planet-only fit to the simulated planet–moon transit light curve (bottom subplots). All panels assume a Jupiter-sized planet at an orbital period of P_b = 30 d, an Earth-like moon with an orbital period of P_s = 3.55 d (= P_Eu), and a sun-like host star. The blue lines refer to an orbital inclination of i_s = 0, that is to say, to co-planar orbital configurations. This results in occultations near moon phases of zero when the moon is in front of the planet as seen from Earth (see Fig. 2), and 0.5 when the planet is behind the moon. Orange lines refer to i_s = 0.2 rad and do not produce planet–moon occultations.

The central panels of Fig. 3, which illustrate the barycentric effects, correspond to the frameworks of Sartoretti & Schneider (1999, for the TTVs), Kipping (2009, for the TDVs), and Heller et al. (2016, for the TTV-TDV figure). However, the contribution of the photometric distortion of the light curves by the moon (upper subplots) results in a very different fit of the TDV, TTV, TRV, and TTV-TDV figures (bottom subplots) than proposed in these latter papers. In particular, the fitted TTV-TDV figure does not at all resemble the ellipse that one would expect based on the barycentric contribution only. As a consequence, the TTV-TDV figure might not be as straightforward an exomoon indicator as proposed by Heller et al. (2016).

The variation of the photometric TDV signal as a function of the moon’s orbital phase (0 ≤ φ_s ≤ 1) shows a relatively complicated behavior with minima when the planet and the moon align along the line of sight (moon phases of 0 or 1) and maxima at moon phases of 0.25, 0.5, and 0.75. For comparison, the variation of the barycentric TDV signal with changing moon phase is strongest when the planet and moon align at either a moon phase of 0 or 0.5. Both contributions combined, this leads to a complicated pattern for the TDVs as derived from the light-curve fits.

The pattern of the photometric contribution to the measured TTV (Fig. 3b) is shifted by half an orbital moon phase relative to the barycentric effect (panel a), though with roughly the same amplitude. The resulting TTV signal in the fitted model light curves then turns out to be extremely small (panel c).

The photometric TRV greatly depends on the projected separation of planet and moon, with φ_s = 0 corresponding to the moon being precisely in front of the planet as seen from Earth, and φ_s = 0.25 corresponding to the moon being separated as much as possible from the planet and as “late” as possible for the transit and so on (see Fig. 2). The larger the sky-projected apparent separation between the planet and the moon, that is, the better the moon’s photometric dip in the light curve is separated from the planet’s photometric signature, the smaller the TRV. On the other hand, if the two bodies are sufficiently close to cause a planet–moon occultation, then the TRV drops significantly (see the dip at moon phase 0.5 in Fig. 3c). In the case of an occultation, there are two extreme scenarios to consider; one in which the occultation persists throughout the entire transit and one in which the occultation starts or ends during the transit. In the first case, the measured transit depth corresponds very closely to the true planet radius. In the second case,the resulting shape of the transit is highly dependent on the exact orbital geometry and dynamics, and the inferred transit depth can differ significantly from the planet-only scenario (see Fig. 3c).

Interestingly, there is a difference in amplitude depending on whether the moon passes behind or in front of the planet. This is due to the slightly different transit shapes of the combined planet–moon system. If the moon’s effective (tangential) speed across the stellar disk is lower than that of the planet, the overall transit shape resembles that of a single planet except for an additional signature in the wings of the transit light curve. We note that this asymmetry of the TRV effect cannot help to discriminate between prograde and retrograde moon orbits, which is both a theoretical and an observational challenge for exomoon characterization (Lewis & Fujii 2014; Heller & Albrecht 2014). Moreover, in practice this effect will be difficult to observe when noise is present. We note that the barycentric TRV shows no variation as a function of the moon phase in the case of i_s = 0 but some variation in the case of i_s = 0.2 rad. This can be explained by the planet crossing the star at a higher or lower impact parameter due to its movement around the planet–moon barycenter. Therefore, at different transit impact parameters, the resulting transit depth will be different due to the stellar limb darkening effect (Heller 2019).

Fig. 3 Exomoon indicators for a hypothetical Jupiter-Earth planet–moon system in a 30 d orbit around a sun-like star as measured from noiseless simulated light curves by fitting a planet-only model. Panel a: transitduration variation (TDV). Panel b: transit timing variation (TTV). Panel c: transit radius variation (TRV). Panel d: TTV vs. TDV diagram. Blue lines refer to a coplanar planet–moon system with is = 0 and orange lines depict the effects for an inclined planet–moon system with is = 0.2 rad. In each panel, the top subplot shows the photometric contribution, the center subplot the barycentric contribution, and the bottom subplot the combined and measured effect.

Fig. 4 TRV, TDV, and TTV effects of an Earth-sized moon in a Europa-like orbit (3.55 d) around a Jupiter-sized planet. The planet–moon barycenter has a 30 d orbit around a sun-like star and the transits (not shown) were sampled with a PLATO-like cadence of 25 s. Left: without noise. Right: with noise contributions from stellar granulation and PLATO-like instrumental noise for an mV = 8 star. The subpanels show the autocorrelation. In all panels, blue symbols refer to a coplanar planet–moon system with is = 0 and orange symbols depict the effects for an inclined planet–moon system with is = 0.2 rad.

3.2 TTV, TDV, and TRV from PLATO-like light curves

The panels on the left of Fig. 4 show the TRV (top), TDV (center), and TTV (bottom) effects, respectively, as measured from a sequence of simulated light curves without noise. Just like in Fig. 1, we modeled the planet–moon transits as such but fitted a planet-only model. The number of the transit in this sequence is shown along the abscissa. Also shown are the respective autocorrelations below each of these three time series. The model systems used for these sequences are our test cases 1b (i_s = 0) and 1j (i_s = 0.2 rad), that is a Jupiter-sized, Jupiter-mass planet with an Earth-sized moon in a 3.55 d (Europa-like) orbit. The center of mass of this planet–moon binary is in a 30 d orbit around a sun-like star (for details, see Table A.1).

The combination of these particular orbital periods around the star (P_b) and around the planet–moon barycenter (P_s) happen to make the moon jump by a phase of 0.448 circumplanetary (or circumbarycentric) orbits between successive transits. This well-known undersampling results in an observed TTV signal period of under 2 P_b (Kipping 2009). The TDV series is also subject to this aliasing effect, but the autocorrelation suggests some periodicity. Finally, the TRV signal, whose period is half of P_s, is not affected by the undersampling. Hence, the periodicity can even be observed by eye in the top panel and indeed also in the corresponding autocorrelation.

Moving on to the addition of noise in our simulated light curves, the panels on the right of Fig. 4 illustrate the same TRV (top), TDV (center), and TTV (bottom) series as in the panels on the left, but now as if the host star were an m_V = 8 solar-type star observed with PLATO. In both the TDV and the TTV panels we see that the magnitude of the error caused by the noise in the light curves is comparable to the TTV and TDV amplitudes. As a consequence, both the two time series and their autocorrelations show only weak hints of periodic variations. For comparison, the TRV signal of the slightly inclined system (i = 0.2 rad, orange dashed line) is still clearly visible, whereas the signal for the i = 0 system is somewhat weaker.

To obtain a more quantitative understanding of the possible periodicities in the TRV, TTV, and TDV sequences and in order to compare the performances of the Kepler and PLATO missions, we investigate the periodograms of the data. In Fig. 5 we present the periodograms of the time series and of the autocorrelations of the TRV (top), TTV (center), and TDV (bottom) sequences. The left column of panels assumes a signal-to-noise ratio (S/N) performance and 30 min cadence akin to the Kepler mission, whereas the right column refers to a 25 s sampling and coverage by all 24 of PLATO’s normal cameras. Each panel is based on 20 randomized noise realizations of 24 transits of the same system as in Fig. 4. In both the Kepler and the PLATO case, the TRVs of each realization show a clear peak at around 9.6 times the orbital period of the barycenter around the star. This peak occurs in the periodograms of both the time series and the autocorrelation of the time series and can be explained as an aliasing effect (see Appendix C for a calculation). However, the periodograms of the TDV series and TDV autocorrelation do not show any significant signals near this period, although the TTVs show a peak in the periodogram close to the period corresponding to the Nyquist frequency.

Of similar interest is our finding that Kepler and PLATO have comparable performances in detection of TRVs and TDVs, but PLATO data produce a more significant signal in the TTV periodogram and autocorrelation. This result is mostly due to the higher cadence of PLATO. Moreover, PLATO will provide many more light curves of bright Sun-like stars than Kepler. For comparison, the brightest exoplanet host star observed with Kepler, HD 179070 (Kepler-21), has an m_V magnitude of 8.3 (Howell et al. 2012) and most Kepler host stars are substantially fainter. As a consequence, any transiting planet with an exomoon in the Kepler data would yield a signal that is significantly weaker than the one shown in Fig. 5.

Fig. 5 Comparison of the simulated performance of Kepler (left) and PLATO (right) in detecting the TRV, TDV, and TTV effects.A sequence of 20 simulated light curves for a Jupiter–Earth planet–moon system around a sun-like star were pre-computed and analyzed, similar to the system used in Fig. 4 (test case 1p in Table A.1). The 20 randomized noise realizations contain both white noise and granulation equivalent to an mV = 8 star. The Nyquist period of the system is equal to 2 Pb. The panels shows the periodograms of the TRV (top), TDV (center), and TTV (bottom) effects and their respective autocorrelations.

3.3 The TRV effect of different star–planet–moon systems

The results presented in Figs. 4 and 5 refer to only two specific test cases, that is, 1b and 1f, respectively. Next, we extend our investigations to all other test cases described in Sect. 2.3 and summarized in Table A.1. Appendix B contains a completed set of plots for each of the subcases, while Fig. 6 shows the results for a subsection of these cases as a cut through the parameter space. All results shown in Fig. 6 assume synchronous observations of an m_V = 8 star with all 24 normal cameras of PLATO. We focus on the TRV because our investigations to this point have indicated that the TRV effect is more pronounced than the TTV and TDV effects.

The four systems in the upper row in Fig. 6 refer to cases 2b (solid orange), 2d (dashed orange), 6b (blue solid), and 6d (blue dashed). The second line of panels shows cases 3b (solid orange), 3d (dashed orange), 7b (solid blue), and 7d (dashed blue). The third line of panels shows cases 4b (solid orange), 4d (dashed orange), 8b (blue solid), and 8d (blue dashed). The bottomrow of panels refers to cases 5b (orange solid), 5d (orange dashed), 9b (blue solid), and 9d (blue dashed).

In cases 2, 3, 6, and 7 involving the large Neptune and super-Earth-sized moons, the TRV signal is always visible in the time series (left column). However, the autocorrelation function of the TRV time series (center column) for systems with P_b = 60 d is truncated at six times the orbital period due to the low number of transits available at that period. This truncation is caused by our restriction of the autocorrelation function to a lag time of at least half the observed number of transits. For comparison, in those cases that assume P_b = 30 d and thus twice the number of transits available for the analysis, the signal in the corresponding periodogram of the autocorrelation is very clear.

For cases 4 and 8 involving an Earth-sized moon, the detectability of the TRV signal depends on the moon’s orbital period and the size of its star. For the case of a Europa-like orbital period of the moon and the ɛ Eridani-like star, a marginally significant peak in the periodogram is present, but no significant peak can be seen if the star is sun-like or the moon’s orbital period is that of Io. We were unable to produce a significant TRV signal in any of the cases with a Moon-sized moon (cases5 and 9).

In the right column of panels in Fig. 6, we observe that our test moons with a Europa-like orbital period all cause a peak in the periodograms at around 9.6 P_b. For all the other cases involving a moon with an Io-like period (all subcases a and c), the position of the peak in the periodgrams of the autocorrelation is near 12 P_b, which is at the edge of the number of orbital periods we deem reasonable for analysis.

The positions of the peak in the periodogram of the TRV autocorrelation function corresponds to the remainder of the fraction between the circumstellar period of the planet and the circumplanetary period of the moon. If the nominal value of this resulting peak position is smaller than two planetary orbital periods, then the actual position of the peak is shifted to a position longer than two times P_b (see Appendix C).

Fig. 6 Transit sequence of the transit radius variation (TRV, left), autocorrelation function of the TRV transit sequence (center), and periodogram of the autocorrelation (right) for various planet–moon configurations around an mV = 8 star asobserved with PLATO (see Sect. 3.3). Blue lines refer to an ɛ Eridani-like star, orange lines to a sun-like star. Subscripts “J”, “S”, “N”, “⊕”, and “☾ ” for the radii and masses refer to Jupiter, Saturn, Neptune, Earth, and the Moon, respectively. In the Jupiter–Neptune and Saturn/Super–Earth cases the TRV signal is clearly visible in the measured TRV time series. Consequently the signal is also visible in the autocorrelation. Due to only calculating the AC up toa time lag of have the observation length, a peak in the periodograms is only visible in the Pb = 30 d cases and not the Pb = 60 d cases. In the Neptune–Earth cases the signal in the TRV time series is barely visible, while the autocorrelation and periodogramfor the Pb = 30 d cases show a clear signal. There is no visible signal in any of the Earth–Moon cases for in the corresponding time series, autocorrelation, and periodogram.

3.4 Kepler-856 b as an example for TRVs

As a first application of our novel TRV method, we searched for TRVs in the published values of the transit depth sequences of 2598 KOIs, or transiting exoplanet candidates from Holczer et al. (2016). As a general result of our manual inspection of those KOIs with the highest periodogram power, we noticed that stellar rotation is the major cause of false positives. However, in the case of Kepler-856 b (KOI-1457.01), the interpretation of a false positive is less evident based purely on the inspection of the light curve. Figure 7 shows our vetting sheet constructed for our TRV survey with Kepler. Kepler-856 b is a Saturn-sized planet with unknown mass in a 8 d orbit around a sun-like star. The transit sequence has been used to statistically validate Kepler-856 b as a planet with less than 1% probability of being a false positive (Morton et al. 2016).

Figure 7a shows the Pre-Search Conditioning Simple Aperture Photometry (PDCSAP) Kepler flux (Jenkins et al. 2010), which has been corrected for the systematic effects of the telescope and normalized by the mean of each Kepler quarter. We note that the apparent periodic variation of the maximum transit depth is probably not the TRV effect that we describe. The apparent variation of the maximum transit depth visible in Fig. 7a is likely an aliasing effect due to the finite exposure time of the telescope.

Figure 7b shows the series of the transit depths of Kepler-856 b, panel c the autocorrelation function, panel d the periodogram power of the transit depth series (blue) and of the autocorrelation function (orange), and panel e shows the periodograms of the 14 Kepler quarters of the light curve (gray) together with their mean values (black). We used panel e as a means to quickly identify stellar variability and compare it to any possible periodicity in the TRV periodogram of panel d.

An interesting observation of Fig. 7d that caught our attention is the striking peak in the periodogram of the autocorrelation function near 200 d and, to a minor extent, a peak in the periodogram of the transit depth variation at the same period. These observed features are reminiscent of the features in the corresponding periodograms of the simulated data shown in Figs. 6, B.1, and B.2, which are caused by large moons. Moreover, the periodogram of the total light curve, which is a combination of all Kepler quarters (red line in Fig. 7e), does not show a peak near 200 d with similar amplitude as the periodograms in panel d, although there are two moderate peaks between 100 d and 200 d. We conclude that there is only marginal evidence for stellar activity as the origin of the TRV signal in panel d.

Could amoon nevertheless be stable around Kepler-856 b from an orbital dynamics point of view? Using the stellar and planetary parameters derived by Morton et al. (2016), we estimate that the planetary Hill radius (R_H) is about 540000 km wide. Given that prograde moons can only be stable out to about 0.5 × R_H (Domingos et al. 2006), the range of prograde and gravitationally stable moon orbits is limited to about 270 000 km, corresponding to about 9 R_p for Kepler-856 b. For comparison, Io, the innermost of the Galilean moons, orbits Jupiter at about 6.1 planetary radii, corresponding to about 8% of Jupiter’s Hill radius. The space for orbital stability is therefore quite narrow around Kepler-856 b compared to the orbits of the Solar System moons, though a close-in sufficiently large moon could still be physically plausible.

Fig. 7 Analysis of the Kepler light curve of Kepler-856 b (KOI-1457.01). Panel a: PDCSAP Kepler light curve. Panel b: fractionaltransit depths variation from Holczer et al. (2016) for a total of 134 transits. Panel c: autocorrelation of the transit depth variation. Panel d: periodogram of the transit depth variation and periodogram of the autocorrelation function. Panel e: Periodograms of the PDCSAP flux of 14 Kepler quarters (gray lines) and their mean (black line) after masking outthe planetary transits.

4 Discussion

Comparing the TTV, TDV, and TRV distributions as a function of the moon’s orbital phase for systems of various orbital inclinations,we noticed that systems without occultations (and significant inclinations) exhibit more sinusoidal patterns. A comparison of the blue (with occultations) and orange (without occultations) lines in any of the upper panels of Figs. 3a–c serves as an illustrative example of this. Configurations without occultations therefore naturally produce higher peaks in the corresponding periodograms. An improvement to our TRV indicator, in particular for systems with occultations, could be achieved if the TRV signal could be modeled and fitted to the TRV time series, or if it could be used as a fitting function in the periodogram. However, this approach could become as mathematically and computationally involved as implementing a planet–moon transit model to the light curve to begin with.

The detectability of moon-induced TRV signals increases with the number of transits observed and it is desirable to at least observe a full period of the TRV signal. Hence, this new exomoon indicator can preferably be detected with longterm transit surveys such as Kepler or PLATO, which deliver (or are planned to deliver) continuous stellar high-precision photometry for several years.

Localized features on a star’s surface such as spots or faculae can also cause variations in the measured exoplanet radius. Hence, any search for moon-induced TRVs needs to include the analysis of stellar variations and verify that any candidate signals are not caused by the star. Analyzing the contributions of the barycentric and photometric components to the observed – or in our case the fitted – TTV signal, we noticed that both components tend to cancel each other; see Fig. 3b. This means that for low-density moons the TTV signal will be dominated by the photometriccomponent, while for high-density moons the barycentric component will be more relevant. For moons with densities similar to the densities of the terrestrial planets and moons of the Solar System, which range between about 3 and 5 g cm⁻³, the overall TTVeffect might actually be very hard to observe as the barycentric and photometric contributions essentially cancel each other out.

The amplitude of the TRV signal increases with increasing orbital separation between the planet and the moon. The TRV amplitude is maximized when the orbital semimajor axis around the planet is larger than two stellar radii, at which point the two transits of the planet and the moon can occur without any overlap in the light curve. Transits without an overlap happen when the moon’s orbital phase is near φ_s = 0.25 or φ_s = 0.75, that is, when the planet and the moon show maximum tangential deflection on the celestial plane. In these cases of a transit geometry, there is no “contamination” of the planetary transit by the moon and the lowest possible transit radius is fitted, thereby maximizing the variation with respect to those transits of a transit sequence in which the planet and the moon move along or near the line of sight, that is, near φ_s = 0 or φ_s = 0.5, at which point the maximum transit radius (with maximum contamination by the moon) is fitted.

5 Conclusions

We use standard tools for the analysis of exoplanet transits to study the usefulness of TTVs, TDVs, and TRVs as exomoon indicators. Our aim is to identify patterns in the data that are now regularly being obtained by exoplanet searches, which could betray the presence of a moon around a given exoplanet without the need for dedicated star–planet–moon simulations.

We first test these indicators on simulated light curves of different stellar systems and find that the TRVs observed over multiple transits show periodic patterns that indicate the presence of a companion around the planet. However, the TTV and TDV signals may be less useful indicators for the presence of an exomoon depending on the density of the moon, because the moon density determines whether the barycentric or photometric component of the respective effect dominates. The TRV signal, for comparison, is always dominated by the photometric contribution and as a consequence has more predictable, understandable variation.

Our search for indications of an exomoon in the 2598 series of transit mid-points, transit durations, and transit depths by Holczer et al. (2016) reveals the hot Saturn Kepler-856 b as an interesting object with strong TRVs that are likely not caused by stellar variability. While an exomoon interpretation would be a natural suggestion to put forward in this study, we caution that the range of dynamically stable moon orbits around Kepler-856 b is rather narrow due to its proximity to the star. The planetary Hill sphere is about 18 R_p in diameter, restricting prograde moons to obits of less than about 9 R_p in width. However, at this proximity to the planet, tides may become important and ultimately present additional obstacles to the longterm survival of a massive moon (Barnes & O’Brien 2002; Heller 2012). Generally speaking, our investigations of the fitted transit depths from Holczer et al. (2016) and of the corresponding Kepler light curves show that rotation of spotted stars or an undersampling effect of the transit light curve can cause false positives.

The previously suggested method of finding exomoons in the TTV-TDV diagram (Heller et al. 2016) might not be as efficient a tool to search for moons as first thought. The TTV-TDV diagram only forms an ellipse in the barycentric regime, that is, when the moon has a very high density. However, for rocky or icy moons, the figure of the TTV-TDV diagram derived from fitting a planet-only model to the combined planet–moon transit light curve is much more complex than an ellipse.

We find that the minimum exomoon radius for a TRV-based exomoon detection around an m_V = 8 solar type star with the Kepler or PLATO missions lies at about the size of the Earth. Our simulations suggest that PLATO and Kepler show a similar performance in the detection of the exomoon-induced TRV effect. That said,Kepler did not actually observe any transiting planets around a star as bright as m_V = 8, whereas PLATO will explicitly target these stars for transiting planets. Moreover, we find that PLATO has a better chance of detecting exomoon-induced TTVs due to its higher cadence.

Acknowledgements

This work was supported in part by the German Aerospace Center (DLR) under PLATO Data Center grant 50OL1701.

Appendix A Full parameter set

Appendix B Plots of cases 2–9

Fig. B.1 As Fig. 6, but for a PLATO-like telescope and a moon period equal to Io’s orbit around Jupiter (top panel) and Kepler-like telescopes and moon periods PEu and PIo (bottom panels).

Fig. B.2 As Figs. 6 and B.1, but for magnitude 11, PLATO and Kepler-like telescopes and moon periods PEu and PIo.

Appendix C Aliasing of an observed TRV

In our model of a star–planet–moon system, the phase of the moon orbit around the planet advances by a constant fraction 0 ≤ r ≤ 1 between successive transits, although the moon orbits the planet many times between transits. This advancement can be calculated from the orbital period of the planet–moon barycenter around the star (P_b) and the orbital period of the satellite around the local planet–moon barycenter (P_s) as $$\frac{P_{\text{b}}}{P_{\text{s}}}=n+r,$$(C.1)

where n is the number of completed orbits between successive transits and r is the remainder (Heller et al. 2016). In a similar way, the remainder of the TRV signal r_TRV can be calculated, but we note that the period of the TRV effect (P_TRV) is half of the orbital period of the moon, P_TRV  =  P_s∕2, as can be seen in Fig. 3 for example². The frequency of the TRV signal then is $$f_{\text{TRV}}=\frac{1}{P_{\text{TRV}}}=\frac{2}{P_{\text{s}}}.$$(C.2)

The TRV signal can only be sampled during transits, a constraint that establishes a sampling frequency f_sam = 1∕P_b and a Nyquist frequency f_Nyq = f_sam∕2.

If $\left|f_{\text{TRV}}\right|$ is larger than f_Nyq then the signal is aliased and a full integer multiple (N) of f_sam is subtracted from f_TRV so that the resulting alias frequency f_{TRV, alias} falls between − f_Nyq and + f_Nyq, $$f_{\text{TRV},\text{alias}}=f_{\text{TRV}}-N{f}_{\text{sam}},\text{\hspace{0.17em}\hspace{0.17em}with} N=\text{round}\left(\frac{f_{\text{TRV}}}{f_{\text{sam}}}\right).$$(C.3)

The resulting period of the alias TRV signal then comes out as $$\begin{array}{ll}{P}_{\text{TRV},\text{alias}}\hfill & =\frac{1}{f_{\text{TRV},\text{alias}}}=\frac{1}{f_{\text{TRV}}-N{f}_{\text{sam}}}=\frac{1}{\frac{2}{P_{\text{s}}}-\frac{N}{P_{\text{b}}}}\hfill \\ \hfill & =\frac{1}{\frac{2}{P_{\text{s}}}-\frac{N}{P_{\text{b}}}}.\hfill \end{array}$$(C.4)

As an example, let us consider the case of a Jupiter-sized planet in a P_b = 30 d orbit around a sun-like star and an Earth-sized moon in an orbit with a Europa-like orbital period of P_s = 3.5512 d as shown in Fig. 5. We find N = 17 and then use Eq. (C.4) to calculate an alias period of the TRV signal as 287.6 d, corresponding to 9.59 P_b. This is where the peak in the TRV periodogram presented in Fig. 5 occurs.

These calculations assume that the TRV signal is composed of a single sine wave. If the TRV signal is not a sine but is composed of multiple sine waves (in the sense of a Fourier transformation), then each frequency of the signal is folded individually, which might lead to alias signals that have a very different form in Fourier space than the original signal.

References

All Tables

All Figures

Fig. 1 Example of a planet-only fit (dashed orange line) to the planet–moon model (solid black line). The star is sun-like, the planet is Jupiter-sized and in a 30 d orbit around the star, and the moon is Neptune-sized and in a 3.55 d orbit around the planet. In this example, the fitted planet-to-star radius ratio (rfit), which goes into the calculation of the transit depth with a power of two, is slightly overestimated due to the photometric signature of the moon. As a comparison between the amplitude of the moon signal and the expected noise, the noise from the Poisson statistics of an mV = 8 star (using a 1 h binning) is 12.6 ppm for PLATO and 6.0 ppm for Kepler. This is two magnitudes smaller than the amplitude of the moon signature but possibly comparable to the effect of stellar noise, from activity for example.

In the text

	Fig. 2 Illustration of the orbital geometry and of our parameterization of the simulated planet–moon systems. φs is the orbital phase relative to the line of sight half-way through the transit.
In the text

Fig. 3 Exomoon indicators for a hypothetical Jupiter-Earth planet–moon system in a 30 d orbit around a sun-like star as measured from noiseless simulated light curves by fitting a planet-only model. Panel a: transitduration variation (TDV). Panel b: transit timing variation (TTV). Panel c: transit radius variation (TRV). Panel d: TTV vs. TDV diagram. Blue lines refer to a coplanar planet–moon system with is = 0 and orange lines depict the effects for an inclined planet–moon system with is = 0.2 rad. In each panel, the top subplot shows the photometric contribution, the center subplot the barycentric contribution, and the bottom subplot the combined and measured effect.

In the text

Fig. 4 TRV, TDV, and TTV effects of an Earth-sized moon in a Europa-like orbit (3.55 d) around a Jupiter-sized planet. The planet–moon barycenter has a 30 d orbit around a sun-like star and the transits (not shown) were sampled with a PLATO-like cadence of 25 s. Left: without noise. Right: with noise contributions from stellar granulation and PLATO-like instrumental noise for an mV = 8 star. The subpanels show the autocorrelation. In all panels, blue symbols refer to a coplanar planet–moon system with is = 0 and orange symbols depict the effects for an inclined planet–moon system with is = 0.2 rad.

In the text

Fig. 5 Comparison of the simulated performance of Kepler (left) and PLATO (right) in detecting the TRV, TDV, and TTV effects.A sequence of 20 simulated light curves for a Jupiter–Earth planet–moon system around a sun-like star were pre-computed and analyzed, similar to the system used in Fig. 4 (test case 1p in Table A.1). The 20 randomized noise realizations contain both white noise and granulation equivalent to an mV = 8 star. The Nyquist period of the system is equal to 2 Pb. The panels shows the periodograms of the TRV (top), TDV (center), and TTV (bottom) effects and their respective autocorrelations.

In the text

Fig. 6 Transit sequence of the transit radius variation (TRV, left), autocorrelation function of the TRV transit sequence (center), and periodogram of the autocorrelation (right) for various planet–moon configurations around an mV = 8 star asobserved with PLATO (see Sect. 3.3). Blue lines refer to an ɛ Eridani-like star, orange lines to a sun-like star. Subscripts “J”, “S”, “N”, “⊕”, and “☾ ” for the radii and masses refer to Jupiter, Saturn, Neptune, Earth, and the Moon, respectively. In the Jupiter–Neptune and Saturn/Super–Earth cases the TRV signal is clearly visible in the measured TRV time series. Consequently the signal is also visible in the autocorrelation. Due to only calculating the AC up toa time lag of have the observation length, a peak in the periodograms is only visible in the Pb = 30 d cases and not the Pb = 60 d cases. In the Neptune–Earth cases the signal in the TRV time series is barely visible, while the autocorrelation and periodogramfor the Pb = 30 d cases show a clear signal. There is no visible signal in any of the Earth–Moon cases for in the corresponding time series, autocorrelation, and periodogram.

In the text

Fig. 7 Analysis of the Kepler light curve of Kepler-856 b (KOI-1457.01). Panel a: PDCSAP Kepler light curve. Panel b: fractionaltransit depths variation from Holczer et al. (2016) for a total of 134 transits. Panel c: autocorrelation of the transit depth variation. Panel d: periodogram of the transit depth variation and periodogram of the autocorrelation function. Panel e: Periodograms of the PDCSAP flux of 14 Kepler quarters (gray lines) and their mean (black line) after masking outthe planetary transits.

In the text

	Fig. B.1 As Fig. 6, but for a PLATO-like telescope and a moon period equal to Io’s orbit around Jupiter (top panel) and Kepler-like telescopes and moon periods PEu and PIo (bottom panels).
In the text

	Fig. B.2 As Figs. 6 and B.1, but for magnitude 11, PLATO and Kepler-like telescopes and moon periods PEu and PIo.
In the text