{'fraction_non_alphanumeric': 0.036924872193820846, 'fraction_numerical': 0.019115358968659703, 'mean_word_length': 4.483394271785497, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Modeling the Orbital Sampling Effect of Extrasolar Moons

The orbital sampling effect (OSE) appears in phase-folded transit light curves of extrasolar planets with moons. Analytical OSE models have hitherto neglected stellar limb darkening and non-zero transit impact parameters and assumed that the moon is on a circular, co-planar orbit around the planet. Here, we present an analytical OSE model for eccentric moon orbits, which we implement in a numerical simulator with stellar limb darkening that allows for arbitrary transit impact parameters. We also describe and publicly release a fully numerical OSE simulator (PyOSE) that can model arbitrary inclinations of the transiting moon orbit. Both our analytical solution for the OSE and PyOSE can be used to search for exomoons in long-term stellar light curves such as those by Kepler and the upcoming PLATO mission. Our updated OSE model offers an independent method for the verification of possible future exomoon claims via transit timing variations and transit duration variations. Photometrically quiet K and M dwarf stars are particularly promising targets for an exomoon discovery using the OSE.


CONTEXT AND MOTIVATION
The race toward the first detection of an extrasolar moon picks up pace. While the first attempts to detect exomoons were byproducts of planet-targeted observations (Brown et al. 2001;Pont et al. 2007;Maciejewski et al. 2010), high-accuracy space-based observations of thousands of transiting exoplanets and planet candidates by Kepler (Borucki et al. 1997;Batalha et al. 2013) now allow for dedicated exomoon searches (Kipping et al. 2012;Szabó et al. 2013;Hippke 2015). Upcoming data from the European Space Agency's CHEOPS and PLATO space missions will provide further promising avenues toward an exomoon detection (Hippke & Angerhausen 2015;Simon et al. 2015).
Exomoon detections will be highly valuable for our understanding of the origin and fate of planetary systems because they probe the substructures of planet formation that is not accessible through planet observations alone. As for the solar system, moons provide key insights into the formation of Earth (by a giant collision; Cameron & Ward 1976;Canup 2012), into the temperature distributions within the circumplanetary accretion disks of Jupiter and Saturn (Pollack & Reynolds 1974;Canup & Ward 2002;Sasaki et al. 2010;Heller & Pudritz 2015a), and into the cause of Uranus' tilted spin axis (by gradual collisional tilting; Morbidelli et al. 2012). As of today, no moon has been confirmed around a planet beyond the solar system. Hence, exoplanetary science suffers from a fundamental lack of knowledge about the fine structure of planetary systems. Heller (2014) recently identified a new exomoon signature in planetary transit light curves, which occurs due to the additional darkening of the star by a transiting moon. This phenomenon, which we refer to as the photometric orbital sampling effect (OSE), occurs in phase-folded transit light curves, because an exomoon's sky-projected position with respect to its host planet is variable in subsequent transits but is statistically predictable for a large number of transits (N 10). 1 The photometric OSE is not sensitive to the satellite mass (M s ), but it is very sensitive to its radius (R s ). This makes the photometric OSE particularly sensitive to low-density, water-rich moons like the three most massive moons in the solar system, Ganymede and Callisto around Jupiter and Titan around Saturn. These moons would hardly be detectable in the available Kepler data by combined TTV and TDV measurements, which are sensitive to roughly Earth-mass moons (Szabó et al. 2013). The most advanced exomoon search as of today, the "Hunt for Exomoons with Kepler" (HEK; Kipping et al. 2012), using a photodynamical model, achieves much better detection limits down to a few Ganymede masses, depending on the mass of the host planet, amongst other parameters (Kipping et al. 2015). However, if giant planets were able to give their fully fledged, icy moon systems a piggyback ride to 1 AU, where super-Jovian planets are abun- Note that the photometric OSE is sensitive to an exomoon's radius but not to its mass. Hence, large (potentially low-density) moons are the most promising targets. The open circle denotes a Mars-mass moon with a Ganymede-like composition as predicted by Heller & Pudritz (2015b). The mass-radius relationship for moons of various compositions is according to Fortney et al. (2007).
dant in radial velocity survey data (Dawson & Murray-Clay 2013;Heller & Pudritz 2015a), then the photometric OSE could be a promising method to find these exomoons. The mass-radius diagram in Figure 1 illustrates the detection threshold for moons transiting photometrically quiet M stars. Note that Mars-mass moons with a water-rock composition similar to Ganymede, Callisto, and Titan are significantly larger than Mars (see the open circle). Hence, they could be detectable around photometrically quiet stars. These moons are predicted to form frequently around the most massive super-Jovian planets (Heller & Pudritz 2015b). Hippke (2015) searched for the photometric OSE in the archival data of the Kepler telescope and found indications for an OSE-like signal in the combined transit light curves of hundreds of planets and planet candidates with orbital periods > 35 d. The original description of the OSE by Heller (2014) was purely analytical and thereby fast, but it made several simplifications: stellar limb darkening was neglected; the planet-moon orbit was assumed to be non-eccentric (e = 0) and noninclined (i = 0 • ) with respect to the circumstellar orbit; the planet-moon system was assumed to transit the star along the stellar diameter; that is, the transit impact parameter (b) was set to 0. Consequently, this model was not broadly applicable. Hippke (2015) performed a mostly numerical study that did include non-circular planet-moon orbits and arbitrary inclinations, but he did not explore a wide parameter range of the OSE.
We here first derive a novel equation for the orbital sampling frequency P s (x) of exomoons on eccentric or- The probability density can be derived as Ps(x) = ds(x)/dx. The special case of a circular orbit is given in Figure 1 of (Heller 2014). bits (e ≥ 0; Section 2.1). We then incorporate our formula into a numerical simulator that models the phasefolded transits of exoplanets with moons in front of stars with limb darkening (Section 2.2). Our simulations are compared to real Kepler data. We also perform purely numerical simulations for a wide range of inclined planet-moon orbits (i ≥ 0), planetary impact parameters (b ≥ 0), satellite radii (R s ), and semimajor axes of the satellite's orbit around the planet (a; Section 3).

A DYNAMICAL MODEL FOR THE PHOTOMETRIC OSE 2.1. Sampling Frequency for Eccentric Exomoon Orbits
In Heller (2014), the OSE was described for circular moon orbits only. Given that the eccentricities of the 10 largest moons in the solar system are all smaller than that of the Earth's Moon (≈ 0.055), this approximation seems appropriate. However, the architectures and physical properties of extrasolar planetary systems turned out to be very different from the solar system planets, as we recall the discoveries of terrestrial planets around pulsars (Wolszczan & Frail 1992), Jupiter-mass planets in extremely short-period orbits (Mayor & Queloz 1995), circumbinary planets (Doyle et al. 2011), and a large population of super-Earths with short orbital periods (Batalha et al. 2013). Hence, astronomers should not be surprised if they were to find that at least some exomoons follow significantly eccentric orbits. Eccentricities may, for example, be forced by gravitational interaction with the star (Heller 2012;Spalding et al. 2016Spalding et al. ), other planets (Ćuk 2007Gong et al. 2013;Payne et al. 2013), or even the host planet (Goldreich 1963). Moreover, an analytical solution is always worthwhile as it can give a deeper understanding of the underlying physical processes that generate an observational feature.
We here derive the orbital sampling probability function for exomoons on eccentric circumplanetary orbits. Figure 2(a) shows the moon's circumplanetary orbit in a plane perpendicular to the reader's line of sight (solid line). If an observer were to take pictures of the moon's orbital position with a constant sampling rate, then the individual moon pictures would be distributed in a manner similar to the one depicted in Figure 2(a). Around periapsis, the moon's Keplerian velocity (v) is higher than around apoapsis. Thus, the probability of observing the moon, e.g. during a common stellar transit with its host planet, at a given orbital position is asymmetric with respect to the sky-projected distance to the planet. This fact is visualized in Figure 2(b), where the planetmoon system is sampled from a co-planar perspective. Due to the projection effect, the moon(s) pile(s) up toward the edges of the projected major axis, but the probability distribution to the left of the planet is different from what it looks to the right of the planet. Figure 2(c) shows how we construct the satellite's probability density P s (x) as a function of its projected distance (x) to the planet.
In its most general form, P s (x) describes what fraction of its orbital period (P ps ) the satellite spends in an infinitely small interval (dx) on its sky-projected orbit along the x-axis (see Figure 2(c)). Hence, it is given as where dt is an infinitesimal change in time. With ds as an infinitely small interval on the moon's elliptical orbit, the Keplerian velocity is v(x) = ds(x)/dt, hence and the challenge is then in finding ds(x)/dx. From Figure 2(a) we learn that ds(x) = r(x)dϕ(x). We use the parameterization of a Keplerian orbit and solve it for ϕ(x) via Equation (2) then becomes where we introduced A≡a(1 − e 2 ). In Equation (5), r(x) can be derived by inserting Equation (4) into (3), hence and v(x) in Equation (5) is given by where µ = G(M p + M s ) and G is Newton's gravitational constant. So far, we assumed the line of sight to be parallel to the y-axis in Figure (2), that is, along the orbital semiminor axis (b = a(1 − e 2 ) 1/2 ). Of course, the moon orbit can be rotated in the x-y plane by an angle ω. We can imagine that once ω = 90 • , the observer samples the moon orbit according to P s (y). In this case, the sky-projected orbit appears symmetric with an apparent radius (ã(e, ω)) equal to b. In general, we havẽ a(e, ω) = a cos(ω) We then replace e withẽ = e cos(ω) and A withÃ = a(1 −ẽ 2 ) in Equations (5)-(7) to obtain P s (x) for arbitrary orientations ω. At last, we can swap the proportionality sign in Equation (5) with an equality sign by normalizing where r a and r p are the moon's orbital radii at apoapsis and periapsis, respectively (see Figure 2(c)). We evaluate this integral numerically for arbitrary a, e, and ω and find an approximate solution P s (x) =2(1 + e 4 ) ã a(e, ω) wherer(x) andṽ(x) refer to Equations (6) and (7), respectively, but swapping e forẽ and A forÃ. 2 The error in Equation 10 is ≪ 1 % for e < 0.9. In Figure 3, we plot Equation (10) for a moon orbiting a planet on various eccentric orbits, all of which have a semimajor axis of 10 planetary radii (R p ). The upper panel depicts the orbital geometries and the lower panel shows P s (x) in each case. Four cases with eccentricities between 0 and 0.4 assume ω = 0 as presented in Figure 2, and the e = 0.5 case is rotated by ω = 90 • . Note that in the latter scenario, where the semiminor axis b is substantially smaller than a and the line of sight is parallel to a, P s (x) is significantly higher at a given planetary separation because the moon occupies a smaller circumplanetary region along the x-axis. This means that the OSE becomes particularly prominent in the phase-folded light curve of eccentric moon systems with ω close to 90 • or 270 • . On the other hand, it becomes relatively weak (or "smeared") for ω near 0 • or 180 • .

The Photometric OSE with Limb Darkening
2.2.1. Dynamical Simulations We built a numerical model to simulate the stellar transit of an exoplanet for arbitrary probability functions P s (x) , including multiple functions in the case of multi-moon systems. As an update to Heller (2014), our simulator now considers stellar limb darkening and the transit impact parameter can be varied. In comparison to Heller (2014, Section 2.2.2 therein), where an analytical solution for the actual light curve (the stellar brightness B (n) OSE due to a transiting planet with n exomoons) without stellar limb darkening has been derived, we do not derive an analytical solution for the light curve of the OSE with limb darkening. Instead, we simulate the OSE by generating a computer model of a limb-darkened stellar disk, approximated as a circle touching the the edges of a square sized 1 000 × 1 000 pixels. Then we let the planet and the probability function(s) of the moon(s) transit. In the following, we refer to this approach as our "dynamical OSE model".
Once P s (x) enters the stellar disk, we multiply the stellar intensity (I px ⋆ ) in any pixel that is covered with the probability density P px s (x) in this pixel, so that dx×I px ⋆ ×P px s (x) is the relative amount of stellar brightness that is obscured in this pixel. In this notation, dx corresponds to the pixel width. Planet-moon eclipses are automatically taken into account, because the planet is simulated as a black circle. Consequently, I px ⋆ = 0 inside the planetary radius and P px s (x) cannot contribute to the photometric OSE in the planetary shadow. At any given observation time t, the sum over all occulted pixels within the stellar radius (R ⋆ ) gives us the relative loss in stellar brightness B 1 The dynamic model without limb darkening reproduces the analytical prediction very well. However, only the dynamic limb darkening model will be useful for fitting real observations. This hypothetical star-planet system is similar to the KOI 255.01 system, except that we here assume b = 0. The toy moon is as large as Ganymede (0.41 R ⊕ ) and the planet-moon orbit is 15.47 Rp wide, equivalent to Ganymede's orbit around Jupiter.
As an example, Figure 4 shows the simulated transit of KOI 255.01 together with the OSE of a hypothetical Ganymede-sized moon at a = 15.47 R p , equivalent to Ganymede around Jupiter. The 2.51 R ⊕ super-Earth (R ⊕ being the Earth's radius) is an interesting object as it transits a 0.53 M ⊙ -mass M dwarf with a radius of 0.51 ± 0.06 R ⊙ every 27.52197994 ± 3.295 × 10 −5 d. 3 Hence, the photometric OSE of even a Ganymede-sized moon could be significant. In our simulations, the transit impact parameter is set to b = 0.0 to ease comparison with the analytic model, although it is really given as b = 0.1244 (+0.2132, −0.1243) in the Exoplanet Archive. The LDCs are a 1 = 0.4354, a 2 = 0.2910, and a 3 = 0 = a 4 . In both panels of Figure 4, the solid line refers to our dynamical OSE model with limb darkening, the dashed line shows the dynamical model with the limb darkening option switched off, and the dotted line shows the analytical OSE model by Heller (2014), which also neglects limb darkening.
The left panel shows that the analytical solution is much less accurate than the dynamical OSE model inside the planetary transit trough because it neglects stellar limb darkening. In the wings of the transit curve, however, the analytical solution without stellar limb dark-ening and the dynamic OSE model with limb darkening differ by < 3 × 10 −6 compared to a maximum OSE depth of about 3 × 10 −5 just before the planetary ingress (right panel). The analytic model might thus offer sufficient accuracy for an initial OSE survey of a large data set. With such an approach, a first and preliminary search for OSE candidates within thousands of phase-folded light curves would be a matter of minutes. Even more encouraging, the dynamic OSE model with limb darkening almost resembles a straight line, at least in this configuration where the moon's semimajor axis is roughly as wide as the stellar diameter. A straight line fit would, of course, simplify an initial OSE search even further, as it would mostly depend on the moon's radius squared (in terms of maximum depth) and on the moon's semimajor axis (in terms of duration).
The right panel of Figure 4 also reveals that the dynamical OSE model causes a slightly smaller depth in the light curve for about two-thirds into the OSE ingress. This effect is caused by stellar limb darkening and the stellar brightness in the occulted regions being lower than the average brightness on the disk. In the final third of the OSE ingress, the dynamical OSE model then yields a deeper absorption because of the increasing stellar brightness towards the disk center. This division into two-thirds, in which the OSE model without stellar limb darkening is deeper than the one with limb darkening, and the one-third where things are reversed, is not a universal relation. Toward the stellar limb, e.g., the model without stellar limb darkening produces a deeper OSE signal during the entire transit.

Comparison with Kepler Data
We now apply our model to observations. This investigation is not an in-depth search for moons around a test planet, but it shall serve as an illustration of the OSE by reference to actual observed data. Our survey for confirmed, super-Neptune-sized Kepler planets with orbital periods > 10 d (to ensure Hill stability of an We first retrieved as many quarters as were available among Q0-Q17 of the long-cadence 4 (30 minutes) publicly available Kepler data for Kepler-229 c. 5 We analyzed the PDCSAP FLUX data, from which the Kepler mission has attempted to remove instrumental variability. Nevertheless, these data still exhibit significant variability unrelated to transits, as seen in Figure 5 (left panel). The creation of the PDCSAP fluxes by the Kepler mission involves the removal of common mode variability from the light curves attributable to instrumental effects, which can distort real astrophysical (such as stellar) variability but primarily at medium to long timescales. Since we consider relatively short-period planets, these possible distortions are unlikely to affect our analysis.
To condition each quarter's observations, we subtracted the quarter's mean value from all data points and then normalized by that mean. To these meansubtracted, mean-normalized data, we applied a mean boxcar filter with a width equal to twice the transit duration plus 10 hr. This window is chosen to maximally remove non-transit variations while preserving the shape of the transit and OSE signals. The right panel of Fig-4 The OSE itself is an averaging effect, as it appears in the phasefolded light curve and only after several transits. Hence, for the OSE curve it does not make a difference if the data is taken in short cadence and then binned into 30 minute intervals or if it is taken in 30 minute intervals in the first place. 5 http://archive.stsci.edu/kepler/data search/search.php ure 5 shows the resulting detrended data. Finally, we stitched together all quarters and masked out 10σ outliers. 6 Figure 6 (right panel) shows the final result for Kepler-229 c, after we folded the detrended data on a period of 16.968618 d. Gray dots present the detrended, phase-folded data, and black dots with error bars show the binned data. The solid black line shows the direct output of our dynamical OSE simulator but for a planet without a moon, and the dashed red line shows the transit including the photometric OSE of an injected moon. A red cross on that curve at about 0.5 hr highlights the orbital configuration that is depicted in the left panel.
The left panel of Figure 6 illustrates our dynamical OSE model at work for Kepler-229 c. The planet along with the probability distribution of one hypothetical exomoon can be seen in transit. The injected moon has a radius of 0.7 R ⊕ , the vertical width of the moons P s (x) is to scale to both the planetary and the stellar radius. The moon orbit is set to 8 R p , which is R H /2 and therefore at the boundary of Hill stability for prograde moons (Domingos et al. 2006).
The inset in the right panel zooms into the transit light curve just about an hour prior to planetary ingress. In addition to the black solid (no moon) and red dashed lines (0.7 R ⊕ -sized moon), we also show a blue dotted line indicating the OSE of a hypothetical Earth-sized moon. Note that the width of individual error bars of the binned data is ≈ 10 −4 , which corresponds to the depth of an Earth-sized moon's photometric OSE. With about 5 of such binned data points (or about 500 unbinned data points during the OSE ingress), a search for the exomoon-induced photometric OSE around this planet could yield statistical constraints on the presence of moons the size of Earth and smaller around this particular exoplanet. A more elaborate statistical analysis of exomoon effects in transit light curves is deferred to a future study (R. Heller et al. 2016, in preparation).

NUMERICAL SIMULATIONS OF THE
PHOTOMETRIC OSE 6 We define σ to be the standard deviation estimated as 1.4826 × the median absolute deviation (Bevington & Robinson 2003). To simulate the OSE in more complex configurations with inclined moon orbits, for which an analytical solution is not available, we wrote a numerical OSE simulator in python, which we call "PyOSE". The code and examples are publicly available 7 under the MIT license. 8 All of the following figures were generated with this code.

Parameterization in PyOSE
In PyOSE, the moon's orbital ellipse is defined by its eccentricity (e), circumplanetary semimajor axis (a), its orbital inclination with respect to the circumstellar orbit (i s ), the longitude of the ascending node (Ω), and the argument of the periapsis (ω). Figure 7 (top) shows a hypothetical 0.7 R ⊕ exomoon around Kepler-229 c on a circular, inclined orbit (i s = 83 • , Ω = 30 • , a = 8 R p ) at the time of the planetary mid-transit. In our numerical implementation, the planet-moon ensemble transits the star from left to right, which is an arbitrary choice. The motion of the planet and the moon around their common barycenter during transit can be simulated with PyOSE.
The center panel in Figure 7 shows a river plot (Carter et al. 2012) representation of 250 of these transit light curves, with one shown in each row. Each of these light curves corresponds to a different fixed position of the moon during the stellar transit. The orbital phase of the moon (along the ordinate in Figure 7) corresponds to the mean anomaly. Two horizontal gray regions at phases ≈ 0.1 and ≈ 0.6 show partial planet-moon eclipses. The The bottom panel of Figure 7 shows a sum of the individual moon transits from the river plot, normalized by the number of transits. This procedure corresponds to a phase-folding of a hypothetical observed OSE.
Although PyOSE can simulate transits of planets with moons, we will focus on the moon's OSE signature in the following and study a range of moon transits for various orbital parameters. The transit model is the one presented by Mandel & Agol (2002). For each OSE curve, we chose to sample at least 100 moon transits equally spaced in time (not necessarily in space for e = 0; see Figure 2) to achieve convergence (Heller 2014). The limbdarkened stellar disk is represented by a numerical grid of 1 000 × 1 000 floating point values, and we calculate the total stellar brightness during transit for 1 000 flux data points between first and last contact of the planetary silhouette with the stellar disk. With the moon having a radius of typically 1/100 the stellar radius in the shown simulations, or about 10 px, its transiting silhouette is represented by about 100 black pixels. The error of this approximation compared to a genuine black circle is < 1 %. PyOSE allows for arbitrarily large pixel grids at the cost of CPU time, which is proportional to the total number of pixels in the grid or to the square of R s (in units of pixels). We treat both the planet and the moon as black disks. Thus, it is irrelevant whether the satellite eclipses behind or in front of the planet. 9 For each of the simulated planet-moon transit configurations, PyOSE compares the distance between the planet and moon disk centers (d ps ) to the sum of their radii. If d ps < R p + R s , a mutual eclipse occurs and PyOSE calculates the area (A) of the asymmetric lens defined by the intersection of the two circles as

Mutual Planet-Moon Eclipses
This area does not contribute to the stellar blocking by the moon, as it is covered by the silhouette of the transiting planet.

Parameter Study
In the following, we study variations of the OSE signal due to variations in the parameterization of the starplanet-moon orbital and physical configuration. We use Kepler-229 c as a reference case and modify one parameter at a time, as specified below.

The Moon's Semimajor Axis (a)
For circular moon orbits, changes in a modify the shape of the OSE, while the area under the integral remains unaffected (see Equation 9). For e = 0 or i s = 0, however, both the shape and the integral will change. The left panel of Figure 8 visualizes this effect for two cases; one in which the moon is barely stable from a dynamical point of view (a = 0.5 R H ), and one for a close-in moon near the Roche limit (a ≈ 0.128 R H ≈ 2 R p ). The duration of the OSE signal is longer for moons with larger semimajor axes. Planet-moon eclipses are visualized by bumps in the OSE curve of the moon at a = 0.128 R H , near ±0.07 days. The moon in the wider orbit is not subject to planet-moon eclipses.
The right panel of Figure 8 shows the integral under the OSE curve as a function of a. This is an important quantity as it serves as a measure for the significance of the moon-induced OSE imprint on the light curve. The decline of the OSE integral up to a ≈ 0.6 R H is mostly due to stellar limb darkening: the larger a, the more flux will be blocked closer to the stellar center for this specific orbital configuration. Beyond 0.6 R H , moon transits will occasionally be missed during planetary transits and the OSE signal decreases. The dashed part of the curve,   Figure 9. Variation of the OSE for different satellite radii. The inset zooms into the ingress of the phase-folded satellite transit, showing that the time of first contact between the satellite's silhouette and the stellar disk depends on the satellite radius.
corresponding to moons beyond 0.5 R H , is physically implausible for prograde moons, and valid only for prograde moons (Domingos et al. 2006).

The Moon's Radius (Rs)
Larger moons naturally cause deeper transits. Keeping everything else fixed, variations in R s cause variations in the OSE amplitude roughly proportional to R 2 s , as illustrated in Figure 9 (see also the term (R s /R ⋆ ) 2 in Equation 11). Comparing the upper and lower curves, we see a signal increase by a factor of four (-25 ppm vs. -100 ppm) for a change in R s by a factor of two (from 0.35 to 0.7 R ⊕ ). Changes in transit duration occur due to the different timings of the first and last contact of the planetary silhouette with the stellar disk (see inset in Figure 9).  2015), this is extremely time-consuming and a practicable approach still needs to be demonstrated. Following Kipping et al. (2014, Sect. 2.3 therein), the HEK team indeed restricts their exomoon search in the Kepler data to single-moon systems. Our Equation (12) is an analytical description of the additive OSE in multi-moon systems and describes how the dynamical OSE model handles multi-satellite systems. Our purely numerical simulator PyOSE can handle multi-satellite systems as well.
That being said, with real observational data, distinguishing between an OSE signal of a single large moon and a signal from several smaller moons may prove difficult, if not impossible. This is demonstrated in Figure 11, where we show the individual OSE signals from two moons (solid lines) and their combined OSE signature (dashed line). This dashed curve demonstrates that the photometric OSE is additive as long as moon-moon occultations can be neglected. 10 The dashed curve is also what would be observed in reality, noise effects aside. With noise taking into account, there would be substantial degeneracies among the model parameters and the number of moons.

DISCUSSION
In addition to the orbital and physical parameters treated above, there are several minor effects on the OSE, the most important being stellar limb darkening. We ran a series of simulations using different LDCs for FGK main-sequence stars. For a typical uncertainty in stellar temperature ≈100K, the effect on the OSE signal is of the order of 1 ppm, that is, negligible for the vast majority of observable cases. Nevertheless, if LDCs cannot be constrained otherwise, there could be degenerate solutions between an exomoon-induced OSE and a sole planetary transit with a different LDC parameterization. This is, however, only an issue during planetary transit. Alternatively, one can search for OSE-like flux decreases before and after the planetary transit (Hippke 2015), which cannot possibly be made up for with different LDCs.
10 PyOSE currently neglects moon-moon occultations. They would only occur in < 1 % of the transits, depending on the exact geometry of the circumstellar and circumplanetary orbits.
Some aspects of our OSE model might hardly be accessible with near-future technology. For moon orbits that are co-planar with the circumstellar orbit (i s = 90 • ), moderate moon eccentricities will not cause an OSE signal much different from a circular moon orbit. Nevertheless, there are configurations in which these parameters make all the difference in determining the presence of an exomoon (see Fig. 10). Moreover, the discovery of the formerly unpredicted hot Jupiter population, the unsuspected dominance of super-Earth-sized exoplanets in orbits as short as that of Mercury, and the puzzling abundance of close-in planets with highly misaligned orbits suggests that the solar system does not present a reliable reference for extrasolar planetary systems. Hence, any model used to explore yet undiscovered exomoons will need to be able to search a large parameter space beyond the margins suggested by the solar system.
Both numerical simulations and observations will always be undersampled and only converge to the analytical solution. In terms of observations, this is because of the limited number of observed transits, usually < 100 for an observational campaign over a few years, and because of telescope downtime and observational windows. Beyond that, the orbital periods of both the planet-moon barycenter (around the star) and the moon(s) (around the planet-moon barycenter) determine a non-randomized sampling of the moon orbit in successive transits. If an exoplanet's orbital period were -by whatsoever reason -an integer multiple of its moon's period, then the moon were to appear at the same position relative to the planet in successive transits. The resulting phase-folded light curve would not display the OSE and therefore not converge to our solutions, but it would show two transits: one caused by the planet and one caused by the moon.
PyOSE can simulate large numbers of OSE curves for arbitrary star-planet-multi-moon configurations to test real observations. Beyond the functionalities demonstrated in this paper, PyOSE can add noise (parameterized or injected real noise).
Both our dynamical OSE model (Sect. 2.2.1) and our numerical OSE simulator (PyOSE, Sect. 3) are computationally inexpensive and easy to implement in computer code, which is crucial for the independent verification or rejection of possible exomoon signals. More advanced methods may suffer from a large parameter space to be explored, resulting in a huge number of simulations (10 11 ; Kipping et al. 2013). These could be difficult to verify. For an independent verification of an exomoon search, either the original code needs to be released for review, or an independent implementation is required. To exclude processing, runtime, and hardware errors (bit error rates are typically 10 −14 ), a search based on TBs of data ultimately needs to be repeated on different hardware. The average computational burden for this is large, with an average of roughly 33,000 CPU hours required per candidate using photodynamical modeling. The monetary equivalent, e.g. using Amazon's EC2 on-demand facilities 11 , is about $50,000 U.S. dollars (2015 December prices) for a single candidate check.
The dynamical OSE simulation in Figure 4 contains 48 data points and was computed within 10 to 14 s on two modern computers. 12 Hence, this setup can generate a grid of about 10 4 such OSE light curves per day. Our OSE model involves 11 independent parameters: M ⋆ , R ⋆ , the planet's orbital period around the star (P ⋆p ), R p , b, a 1 , a 2 , a 3 , a 4 , R s , and a. For a well parameterized starplanet system, M ⋆ , R ⋆ , P ⋆p , and b can be observed and fit without considerations of any potential moon, assuming that moon-induced variations of the transit impact parameter (Kipping 2009) are negligible. If one were to carry out a search for the photometric OSE in the Kepler data, a limb darkening law with two LDCs should still do a good job in a first broad survey, leaving us with R p , R s , and a plus the two LDCs to be fit per phase-folded light curve. Testing 10 values per parameter would then imply a grid of 10 5 OSE light curves per planet or planet candidate. With 10 4 LCs simulated per day, one Kepler planet or planet candidate could be checked for an OSE signature within a week, given a standard desktop computer. If dedicated high-speed computational resources could be used, all Kepler planets and candidates (about 4000 as of today) could be checked for a photometric OSE within maybe a month. A simple straight line fit of the OSE, as suggested above, would dramatically reduce this time frame to much less than one day. A detailed statistical analysis, e.g. within a Bayesian framework (Kipping et al. 2012) and using an injection-retrieval method (Hippke 2015), could then be used to infer the significance of the best-fit model for each object.

CONCLUSION
We present a new formula to describe the orbital sampling frequency of a moon on an eccentric orbit around a planet, that is, the probability of a moon residing at a specific sky-projected distance from the planet (Equation 10). This formula assumes co-planar circumstellar and circumplanetary orbits. We implemented it in a dynamical OSE simulator with stellar limb darkening that can be applied to arbitrary transit impact parameters. In contrast to a previously derived framework, in which stellar limb darkening was neglected and b was required to be zero (Heller 2014), our new dynamical OSE simulator can now be applied to observations ( Figure 6).
Using an independent numerical OSE simulator dubbed PyOSE, we examined the moon's part of the OSE parameter space, spanned by its orbital semimajor axis, its physical radius, the inclination of its orbit with respect to the line of sight, the orbit's longitude of the ascending node, and the transit impact parameter of the planet (and therefore of the moon).
The OSE might give evidence of a multi-moon configuration, but the precise characterization of multi-satellite system will be extremely challenging using OSE only. Therefore, the OSE method could be used for preliminary analyses of a large number of systems, while more costly methods (Kipping et al. 2012;Heller & Albrecht 2014;Agol et al. 2015) could be used to focus on the most promising subset of targets. Beyond that, the OSE method can generally be used as an independent means to verify an exomoon claim via planetary TTV and TDV.