© ESO 2022

1 Introduction

While the search for planets beyond the Solar System (exoplanets) has greatly benefited from the transit method in search for periodic occultations of a star by its planets (e.g., Batalha et al. 2013; Zeng et al. 2019), no moon around any exoplanet has been securely discovered as of today. With over 3400 exoplanets found by the transit method, we naturally wonder if any exo-moons could be detected around them. The detection of exomoon transits is very difficult due to the complex orbital motion, the occurrence of combined transits of an exoplanet with its moon (or moons), possible planet-moon eclipses, and potentially due to the small physical radii of exomoons, all of which make phase-folding approaches similar to those used for transiting exoplanets (Kovâcs et al. 2002; Hippke & Heller 2019) inefficient. Indeed, methods like the orbital sampling effect (Heller 2014) or transit origami (Kipping 2021) neglect about half of the tiny signal, details depending on the planet-moon orbital separation, and the transit geometry (Heller & Hippke 2022).

Photodynamical modeling, on the other hand, is computationally extremely costly as it consists of up to 19 free parameters for inclined, eccentric orbits. Kipping et al. (2013) report 49.7 yr of CPU time for a single system using a similar computing infrastructure and multimodal nested sampling regression. Kipping et al. (2015) report a mean value of 33 000 h of CPU time per planet-moon system on an unquantified number of AMD Opteron 6272 and 6282 SE processors. These times were recorded in a search for transiting exoplanets with exomoons in a sample of 41 light curves from quarters Q0–Q17 of the Kepler mission using the MultiNest fitting software (Feroz et al. 2009).

Beyond the LUNA Fortran code (Kipping 2011) for transiting exoplanet-exomoon systems, other algorithms for the modeling of mutual transiting and eclipsing bodies exist. planetplanet (Luger et al. 2017), software written in C and wrapped in a python interface, is taylored to model multitransiting exoplanet systems and it includes an analytical framework for the modeling of planet-planet mutual eclipses. Given the recent announcements of transiting exomoon candidates around Kepler-1625 (Teachey et al. 2018) and Kepler-1708 (Kipping et al. 2022) and the lack of available computer code to investigate these and future exomoon claims independently, an open-source exomoon transit detection algorithm is very desirable for the community.

Here, we present Pandora, the first open-source exomoon transit detection algorithm¹ using a full photodynamical model. This state-of-the-art algorithm for exomoon transit modeling and searches requires just a few hours of CPU time for a Kepler-like light curve, which is about four orders of magnitude faster than previous algorithms, of which a factor of a few is due to today’s faster CPUs. Pandora’s speed gain is made possible through efficient computer code implementation in the python programming language (e.g., using a just-in-time compilation of python code using numba, Lam et al. 2015), a number of acceptable physical simplifications (e.g., eccentricity approximations), the symmetry of prograde and retrograde orbits, modern sampling algorithms, and other computational advances.

2 Physical model

2.1 Parameterization

At the heart of Pandora sits a photodynamical model of a planet with a single moon that transits their common host star. Both the planet (with mass M_p) and its moon (with mass M_m) orbit their common center of mass (with mass M_b = M_p + M_m) in Keplerian orbits. The orbit of the planet–moon barycenter around the star is parameterized by the midtime of the first barycenter transit in a given transit sequence (t₀), the orbital period of the barycenter around the star (P_b), the semi-major axis of the barycenter orbit around the star (a_b), and the transit impact parameter of the barycenter (b_b). In general, the transit impact parameter can be parameterized in terms of the orbital inclination with respect to the line of sight as per b_b = a_b tan(π/2 – i_b)/R_s, where i_b = π corresponds to an edge-on perspective. For a summary of all parameters, see Table 1.

Moving on to the parameterization of the planet-moon system, the planet has a radius Rp and the moon has a radius Rm. In Pandora, spatial dimensions are measured in units of the stellar radius and time is measured in units of days. As a consequence, the planet and moon radius, as well as the barycenter’s orbital semi-major axis, are used as fractions of R_s, that is, as (R_p/R_s), (R_m/R_s), and (a_b/R_s). By default, the orbit of M_b is assumed to have zero eccentricity (e_b = 0), but users can choose e > 0 and then they would need to define the orientation of the periastron with respect to the line of sight (ω_b) as well.

The orbits of the planet and the moon around their local center of mass are modeled as Keplerian orbits. By default, the planet–moon orbital eccentricity is zero, but users can parameterize nonzero eccentricities (epm). In the default mode, the planet–moon orbit is modeled with an additional five parameters: orbital inclination i_pm, longitude of the ascending node (Ω_pm), semimajor axis (a_pm), orbital period (P_pm), and time of periap-sis passage (τ_pm). In Pandora, τ_pm is normalized with respect to P_pm, that is τ_pm € (0,1). This choice is motivated by a significant boost in the convergence speed of our Monte Carlo sampling method (Sect. 3). For the calculation of the position of the planet and the moon on their elliptical orbit, however, we used τP_pm. For eccentric orbits, the number of orbital elements increases to six, including e_pm and the argument of periapsis (ω_pm). Finally, the moon mass M_m is required to model the motion of the planet and moon around their joint barycenter.

The star is parameterized by R_s and two limb darkening coefficients (u₁u₂) for the quadratic limb darkening law (Manduca et al. 1977). The quadratic limb darkening law is widely used in the exoplanet community because it reproduces stellar limb darkening sufficiently well for modern applications with space-based high-accuracy stellar photometry and because Mandel & Agol (2002) derived an analytical solution to the resulting light curve for transiting planets. The stellar mass (M_s) follows directly from the barycentric orbital mean motion n_b = 2π/P_b and M_b via Kepler’s third law as $M_s=n_{\text{b}}^2{a}_{\text{b}}^3/G-M_{\text{b}}$ where G is the gravitational constant.

2.2 Circumstellar orbital eccentricity

Our three-body system is a typical nested Keplerian system. For one thing, we assume that the orbits of the planet and the moon are not perturbed by the stellar gravitational force and that they orbit their common local barycenter in eccentric orbits. For another thing, we assume that the planet–moon barycenter follows a Keplerian orbit itself around the star.

As for the orbit of Mb, we are not interested in the full orbital revolution, but instead focus on an orbital section around the periodic transits. In our frame of reference (Fig. 1), the line of the periapsis and the line of sight embraced an angle ω_b, while the actual position of Mb on its orbit relative to the periapsis is given by the true anomaly (f). We plotted the orientation of the periapsis along the x-axis, in which case the orbital velocity of M_b is given as (Murray & Dermott 1999)

$$\upsilon =\left({\upsilon }_x,{\upsilon }_y\right)=\frac{n_{\text{b}}{a}_{\text{b}}}{\sqrt{1-e_{\text{b}}^2}}\left(-\sin \left(f\right),e+\cos \left(f\right)\right).$$(1)

We are interested in the velocity component that is tangential with the celestial plane (perpendicular to the observer’s line of sight), that is to say, the sky-projected velocity. We describe this direction by a unit vector e_t and we introduced a unit vector e_r for the (negative) radial velocity component: $$e_{\text{r}}=\left(\cos \left(f\right),\sin \left(f\right)\right)$$(2) $$e_{\text{t}}=\left(-\sin \left(f\right),\cos \left(f\right)\right).$$(3)

The projection of u on e_r, which we refer to as υ_r, is given as $${\upsilon }_{\text{r}}=\left|\upsilon \right|\cos \sphericalangle \left(\upsilon ,e_{\text{r}}\right)e_{\text{r}}=\left|\upsilon \right|\cos \left(\alpha \right)e_{\text{r}},$$(4)

where a = <(υ, e_r) (Fig. 1). The cosine of this angle between υ and e_r can be calculated as $$\cos \left(\alpha \right)=\frac{\upsilon e_{\text{r}}}{\left|\upsilon \right|\left|e_{\text{r}}\right|}=\frac{e_{\text{b}}\sin \left(f\right)}{\sqrt{\sin {\left(f\right)}^2+{\left(e_{\text{b}}+\cos \left(f\right)\right)}^2}}.$$(5)

Plugging Eq. (5) into Eq. (4) yields $${\upsilon }_{\text{r}}=\frac{n_{\text{b}}{a}_{\text{b}}{e}_{\text{b}}}{\sqrt{1-e_b^2}}\left(\sin \left(f\right)\cos \left(f\right)\sin {\left(f\right)}^2\right).$$(6)

Analogously, projection of u on e_t, which we refer to as υ_t, yields $${\upsilon }_{\text{t}}=\left|\upsilon \right|\cos \sphericalangle \left(\upsilon ,e_{\text{t}}\right)e_{\text{t}}=\left|\upsilon \right|\cos \left(\beta \right)e_{\text{t}},$$(7)

with $$\cos (\beta )=\frac{\upsilon e_{\text{t}}}{\left|\upsilon \right|\left|e_{\text{t}}\right|}=\frac{1+e_{\text{b}}\cos (f)}{\sqrt{\sin {(f)}^2+{(e_{\text{b}}+\cos (f))}^2}}.$$(8)

Plugging Eq. (8) into Eq. (7) yields $${\upsilon }_{\text{t}}=\frac{n_{\text{b}}{a}_{\text{b}}}{\sqrt{1-e_{\text{b}}^2}}(1+e_{\text{b}}\cos (f))(-\sin (f),\cos (f)).$$(9)

For our purpose of using the midtransit tangential velocity, we set f = ω_b and found $$\left|{\upsilon }_{\text{t}}\right|=\frac{n_{\text{b}}{a}_{\text{b}}}{\sqrt{1-e_{\text{b}}^2}}(1+e_{\text{b}}\cos ({\varpi }_{\text{b}})),$$(10)

and using $p=2R_s\sqrt{1-b_{\text{b}}^2}$ as the transit path across the star, we constructed the transit duration of a point-like M_b as $$d=\frac{p}{\left|{\upsilon }_{\text{t}}\right|}=\frac{2R_{\text{s}}\sqrt{1-b_{\text{b}}^2}\sqrt{1-e_{\text{b}}^2}}{n_{\text{b}}{a}_{\text{b}}(1+e_{\text{b}}\cos ({\varpi }_{\text{b}}))}.$$(11)

In Appendix A we extend these considerations to the tangential transit velocity during ingress and egress and demonstrate that our assumption of constant in-transit velocity of the planet-moon barycenter produces errors of much less than 1 part per million (ppm) even under very unfavorable conditions. Most importantly, Eq. (11) is computationally much fasterthan using a Kepler solver for the barycentric orbital motion around the star.

Fig. 1 Orbit geometry of the planet-moon barycenter (Mb) on its eccentric circumstellar orbit as used in Pandora. The situation depicted in this illustration corresponds to the midtransit time, when the true anomaly coincides with the orientation of the periastron (ωb). Nothing is to scale in this illustration. The bottom part of the ellipse has been spared for illustration purposes.

2.3 Planet and moon orbits around the local barycenter

The relative positions of the planet and its moon with respect to their local center of mass are calculated using a 2D Kepler solver. By default, Pandora assumes a circular orbit as the standard parameterization for large-scale exomoon surveys in hundreds to thousands of light curves with known transiting exoplanets. With e_pm = 0 and ω_pm = 0, the 2D Kepler solver requires four of the usual six Keplerian elements plus the time of any given data point as input parameters. It then returns the positions of both bodies using the analytical solution of the circular Kepler orbit and without any costly iterative approximations (Appendix C.3). Hence, in Pandora’s default mode, the CPU run time of the 2D Kepler solver is as costly as any analytic approximation.

For eccentric planet–moon orbits, which might be interesting to study in more detail for potential exomoon candidates or to examine eclipse phenomena, etc., Pandora solves the Kepler equation using the approach by Markley (1995) similarly to the implementation in PyAstronomy (Czesla et al. 2019). The solution is built on a fifth-order refinement of a cubic equation. It can be executed in a single iteration and requires the calculation of four transcendental functions.

The star is at the origin of our (x, y) coordinate system and all distances are measured in stellar radii. The distance of the center of any transiting body (planet or moon) from the center of the stellar disk, which is required for the flux calculations, is thus simply given as $\sqrt{x^2+y^2}$.

Before calculating any data points for the Kepler ellipse, however, Pandora estimates if the planet–moon barycenter is sufficiently close to the stellar disk for any transits to occur in the first place. If the planet-moon barycenter is farther away from the origin of the coordinate system than the sum of the stellar radius and the maximum possible deflection of the moon from the planet, $\left(1+a_{\text{pm}}\right)e_{\text{pm}}$ then occultations are geometrically impossible. In this case, no calculations are made and the flux in the model light curve is set to 1. This saves CPU time out of transit.

Pandora’s planet-moon orbit module has a maximum throughput between 14 million (if e > 0) and 21 million (if e = 0) data points per second per core on an Intel i7-1185G7 processor.

2.4 Stellar flux and limb darkening

To calculate the apparent stellar flux during transits of the planet-moon pair, we used the analytical equations of Mandel & Agol (2002). Our code was adapted from PyTransit (Parviainen 2015), which was released under the GPL open source license. Pandora’s transit module requires the distance between the centers of the star and the transiting body (provided as a series of values by Pandora’s Kepler ellipse module), the radius of the transiting circle in units of the stellar radius, and the quadratic limb darkening parameters. It returns the series of flux values where an unocculted star is set to unity.

To speed up calculations, Pandora uses a hybrid occulta-tion model, where the small-body approximation (Mandel & Agol 2002) where constant stellar limb darkening behind the occulted area is assumed whenever appropriate. The error of the small-body approximation is most pronounced during transit ingress and egress. But even then, it is typically on the order of 1 ppm, as long as R_p/R_s < 0.01, and it only accounts for about 1% of the transit duration, with details depending on the transit impact parameter. For bodies with sizes 0.01 < R_p/R_s < 0.05, Pandora employs a hybrid model. It calculates two exact values and linearly interpolates intermediate values from the small-planet approximation. The resulting model fluxes are accurate to <1 ppm in all cases (see Appendix C.2).

In case of fixed limb-darkening parameters during sampling, Pandora creates a lookup table of occultation values before the sampling commences. This procedure calculates a 2D grid of occultation values covering planet-to-star radius ratios 0 < (R_p/R_s) < 0.2 and distances from the stellar center 0 < z < 1 + (R_p/R_s) with 300 values along each axis, which we verified to result in errors <1 ppm. Afterwards, the occultation values are calculated for each model from this rectilinear 2D grid with a bilinear interpolation from their nearest neighbors. This procedure only requires four floats to be obtained from memory, and at most eight multiplications, one division, and thirteen additions per data point.

The throughput of Pandora’s transit algorithm is 16million data points per second per core on an Intel i7-1185G7 processor for R_p/_s > 0.05, 48million for R_p/R_s = 0.02, and 110million for R_p/R_s < 0.01. For cached values (in the case of fixed limb darkening), the maximum throughput is 200 million points per second, which is more than ten times the number of a regular calculation as per Mandel & Agol (2002). In combination with our approximation for the transit duration in eccentric orbits (Sect. 2.2), fitting planet-only models with Pandora is up to fifty times faster than batman (Kreidberg 2015). We show a performance comparison of all methods as a function of the number of data points in Appendix C.2.

2.5 Planet-moon eclipses

As long as the planet and the moon transit the star without mutual eclipses, Pandora treats these synchronous transits separately, as if the star were occulted by two independent spherical bodies. If, however, the circles of the planet and the moon overlap during ingress or egress of at least one of the two bodies, then things become more complicated. Hence, planet-moon eclipses during ingress or egress require special treatment.

In his derivation of the equation for the area of common overlap of three circles, Fewell (2006) identified nine cases for the intersection of three circles in a plane. His case (1) presents a circular triangle, for which an analytic solution did not exist until then. Kipping (2011) extended this system to 27 cases for various star-planet-moon arrangements plus specific subcases, all of which are distinguished by specific conditions of the relative distances and radii of the three circles. Our approach with Pandora is much more pragmatic, though still computationally affordable and numerically accurate to arbitrary levels. Pandora distinguishes only between two cases of a planet-moon eclipse and still covers all possible cases of previous studies.

2.5.1 Planet-moon eclipses and no contact with the stellar limb

Pandora’s eclipse case (1) is called if the planet and the moon are both fully on the stellar disk during a mutual event, that is, if $$1-\sqrt{x_{\text{m}}^2+y_{\text{m}}^2}<R_{\text{m}}\text{\hspace{1em}\hspace{0.33em}and\hspace{1em}\hspace{0.33em}}1-\sqrt{x_{\text{p}}^2+y_{\text{p}}^2}<R_{\text{p}},$$(12)

while the distance between their centers (d_pm) is smaller than the sum of their radii, $$d_{\text{pm}}<R_{\text{P}}+R_{\text{m}}.$$(13)

Pandora ignores the Rømer effect for exomoons (dubbed “transit timing dichotomy”; Heller & Albrecht 2014), which is on the order of a few seconds for realistic planet–moon systems. As a consequence, it is irrelevant for the computation light curve whether the moon is physically located in front of the planet or behind the planet during eclipse.

The decrease in the stellar apparent brightness due to the planet and the moon were initially calculated separately. Then we analytically calculated the area of the asymmetric lens defined by the intersection of the two circles (see Eq. (13) in Heller et al. 2016a) and, assuming that the stellar limb darkening is constant under this small area, compensated for the twofold consideration of this small area analytically in the initial calculation. The resulting error is <0.1 ppm. The throughput of Pandora’s eclipse case (1) is 2 billion data points per second per core.

2.5.2 Planet-moon eclipses in contact to the stellar limb

Pandora’s eclipse case (2) occurs if at least one of the transiting bodies touches the stellar limb during a planet-moon eclipse. Then the area of common overlap of the three circles is integrated numerically and the correction of the light curve is performed as for Pandora’s eclipse case (1) assuming constant stellar limb darkening for this small area.

This numerical treatment is illustrated in Fig. 2. Panel a shows a sketch of three overlapping circles, where the black line depicts the stellar limb, the blue circle the planet, and the orange circle the moon. The dashed rectangle denotes the region shown in the zoom in panel a, where we show a pixelated representation of the moon. In this particular case, the moon has a diameter of 25 pixels. In Appendix B we demonstrate that this resolution is sufficient to reduce the resulting numerical error below 1 ppm. The different colors for the different areas of overlap in Fig. 2b were chosen solely for illustrative purposes. In particular, the white area in the upper left region is the area of common overlap that Pandora’s numerical eclipse module determines and compensates for in this specific example.

The throughput of Pandora’s eclipse case (2) is 0.25 million data points per second and per core. This is relatively slow compared to Pandora’s out-of-eclipse throughput of between 16 and 200 million data points per second and per core (Sect. 2.4). However, since planet-moon eclipses are only relevant for planet-moon systems with orbital inclinations less than a few degrees to the line of sight and since they are usually very short and only relevant to a few data points, this Pandora module only accounts for a small fraction of Pandora’s total computing time.

2.6 Supersampling

If the exposure time of a light curve is longer than the time scale of the variation due to the signal, then this signal becomes temporally smeared (Kipping 2010). In Pandora, this temporal smearing can be compensated for by using supersampling of the light curve, which naturally comes at a higher computational cost and which is thus optional for users.

The supersampling factor can be set to any integer ≥1, where a factor of one corresponds to no supersampling, and a factor of five to seven is a sensible upper limit beyond which additional accuracy gains become marginal. The supersampling factor has a near-linear influence on the model calculation time.

3 Results

3.1 Pandora software

Pandora is a pure python software implementation of the modules describes in the previous sections. As for its working procedure, it first calculates the position of the planet–moon barycenter for each timestamp t (Sect. 2.2). Then the Kepler ellipse for the planet–moon binary is computed in order to determine the (x, y) positions of the two bodies with respect to the stellar disk (Sect. 2.3). With these positions, the occultation function is called to construct the flux drop resulting by both bodies separately (Sect. 2.4). Finally, a check for mutual eclipses and a correction to the flux is made if necessary (Sect. 2.5).

In total, Pandora consists of 640 source lines of code (SLOCs) in the python programming language. This includes the main module and all functions. In addition, Pandora uses a class-based wrapper function (56 SLOCs) and a video generator (83 SLOCs). Pandora only depends on the numpy library (Harris et al. 2020) and the numba just-in-time compiler. If Pandora’s own occultation code was replaced, for example, with a call to PyTransit to obtain the occultation values, only 318 SLOCs would remain. This makes Pandora a comparably small code base. Most other packages have a lot more lines of code, such as Wōtan (Hippke et al. 2019, Hippke et al. 1406 SLOCs), TLS (Hippke & Heller 2019, 1689 SLOCs), or batman (Kreidberg 2015, 256 SLOCs in python and 976 SLOCs in C), not counting test files. A smaller code base has the advantage of making a package easier to understand, debug, and extend. While being concise, Pandora is also readable with clear variable names, a logical modular structure, and comments which explain each code line if necessary.

For debugging purposes and to test the correct implementation of our model in Pandora, we compared single-body occultations with limb darkening and super-sampling against PyTransit and batman, and the planet-moon orbit against PyAstronomy, including required coordinate transformation to solve for the Kepler ellipse. All test cases agreed to within 0.01 ppm per data point. We also created a series of test light curves and the corresponding video animations, and stepped through these frame by frame to validate each data point.

Pandora is concise and modular. Its individual components can be readily understood, tested, and revised as needed. Pandora is an open source code and we invite the community to run their own experiments and submit feature requests, ideas for improvement, and bug reports directly to the GitHub tracker². We have previously made very positive experience with community contributions to TLS and Wōtan, for which a total of about 30 improved versions have been released over the course of the first two years based on community pull requests (i.e., code contributions), bug reports (usually for small issues and edge cases), and general modifications to improve the workflow.

Pandora can generate light curves from custom orbital geometries. These light curves can have arbitrary time sampling and the data quality can be artificially reduced with white noise, for which the amplitude can be set by the users. Pandora can also create transit animation videos including a star with realistic limb darkening and color (as per Harre & Heller 2021) at a custom resolution and frame rate. The light curves can be used for exomoon searches. The videos can be useful for education and outreach. Pandora can also output the spatial coordinates, which we found useful for testing and debugging. They could also be used to compare to observational radial velocity data.

Fig. 2 Planet-moon eclipse in contact with the stellar limb. (a) The planet (blue circle) and the moon (orange circle) are both in front of the star (black circle segment) and in eclipse. (b) Numerical pixel grid representing the moon. Different areas of overlap are shaded in arbitrary colors. In this specific configuration, the white area in the upper left region of the image is determined by pixel area summation and then compensated for in the light curve. The right half of the moon is not transiting.

3.2 Pandora output

In Fig. 3, we show an example of an exoplanet-exomoon system modeled with Pandora. The star is assumed to have a solar radius and solar-type limb darkening, with u₁ = 0.4089 and u₂ = 0.2556 as per Claret & Bloemen (2011) for a star with a surface gravity of log(g) = 4.5, and an effective temperature of T_eff = 5750 K, a metallicity of [M/H] = 0.0, and zero microturbulent velocity. The planet-moon barycenter has an orbital period of 365.25 days around the star and a transit impact parameter b = 0.3. The planet is as large and as heavy as Jupiter. The moon is as large and as heavy as Neptune, thus (M_m/M_p = 0.05395), with an orbital period of 1.28 days around the planet-moon barycenter.

This period, which implies an orbital semi-major axis of 5 Rp, ensures that the moon is beyond the Roche radius of the planet, which we determined to be at about 2.3 R_p (Weidner & Horne 2010) under the assumption of fluid-like objects. For comparison, Io, the innermost of the Galilean moons around those of Jupiter, is at roughly 6.1 R_p around Jupiter. The planet-moon orbit is slightly inclined with respect to the line of sight (i_pm = 80°). We do not claim that our model system is particularly representative of real exomoons. However, it is physically possible and well-suited to illustrate the features of Pandora.

In addition to the model light curve, we generated a simulated observation that is roughly representative of an idealized, photometrically quite, Sun-like star as if it were observed by the PLATO mission. PLATO is expected to launch in 2026 in search for transiting exoplanets around bright stars (Rauer et al. 2014). We chose a cadence of 10min as will be used for PLATO’s P5 stellar sample of ≥ 245 000 dwarf and subdwarf stars of spectral types F5-K7 with apparent visual magnitudes 11 ≤ m_V ≤ 13. We added a white noise component to each data point that is randomly drawn from a normal distribution with a standard deviation of 100 ppm. This noise component is roughly representative of a photometrically quiet Sun-like star from PLATO’s P5 sample³, but it ignores any kind of astrophysical variability. Specifically, the photon noise from this test star would correspond to an apparent visual magnitude m_V ~ 11, though the details of the photon noise floor also depend on the actual number of cameras (6, 12, 18, or 24) with which the target is observed. These simulations assume ideal conditions for the purpose of illustrating Pandora’s functionality. For details on the actual detectability of exoplanet transits with PLATO in the presence of stellar noise, see Heller et al. (2021).

Figure 3a shows a snapshot of the video animation generated with Pandora that is taken during midingress of the moon. We configured the system in such a way that the planet is already in transit when the moon begins its transit. Figure 3b refers to the midtransit time of the planet-moon barycenter, at which point the planet and the moon exhibit near-maximum tangential deflection with respect to the line of sight. In Fig. 3c, the planet is in egress during an eclipse, in which case Pandora switches from the analytical treatment of the light curve to numerical simulations. Figure 3d presents the resulting transit light curve model (solid black line) and simulated PLATO observation (black dots) as well as a breakdown of the contributions from the moon (dotted-dashed orange line) and the planet (dashed blue line). Table 2 summarizes the computing time that Pandora spent for the various processes.

Fig. 3 Output demo of Pandora for a system of a Sun-like star, a Jupiter-sized planet in a one-year orbit around that star, and a Neptune-sized moon in a 1.28 days orbit around the giant planet (for details see Sect. 3.2). (a)–(c) Video renderings of the transit. (d) Light curve. The dotted-dashed orange line shows the transit light curve of the exomoon and the dashed blue line shows the transit light curve of the exoplanet using the analytical solution in both cases. The yellow parts of the moon and planet light curves illustrate numerical simulations. The solid black line shows the combined model. Black dots show a simulated observation roughly representative of an mV ~ 11 Sun-like star from the PLATO mission. Digital star colors are from Harre & Heller (2021).

3.3 Simulated exomoon search

To demonstrate Pandora’s performance as an exomoon search tool, we first simulated a single light curve based on the same star-planet-moon system as described in Sect. 3.2, but we extended the data set to four transit epochs. The resulting mock observations with PLATO are shown with black dots in Fig. 4.

Then we fit the simulated PLATO data using UltraNest (Buchner 2021), a general-purpose Bayesian inference package for parameter estimation and model comparison. It allows one to fit arbitrary models for multimodal or non-Gaussian parameter spaces as can be expected for exoplanet-exomoon systems due to alias effects in the orbital period of the planet-moon system. UltraNest also features computer parallelization and the resumption of incomplete runs.

We chose uninformed flat priors for all parameters, but fixed limb darkening to the injected values. All eccentricities were set to zero. As for the sampling method, we ran the StepSampler of UltraNest with 4000 accepted steps until the sample was considered independent, and we let the sampler decide the number of currently active walkers, which typically resulted in hundreds to thousands of active walkers.

Upon convergence after about 5 h of run time on a single core of an Intel i7-1185G7 processor, we found the solution shown in the corner plot of Fig. 5. The sampler performed roughly 250 million model evaluations in total, which equates to an average of about 13 900 models and log-likelihood calculations per second. We expect linear scaling with more cores and more light curves due to the independence of a search in available data, for example when searching for moons in Kepler or TESS data.

In Fig. 4, we show 100 randomly drawn models from the posterior distribution. The transparent black lines illustrate the total model, the transparent orange lines refer to the transit contribution by the exomoon, and the transparent blue lines show the stellar dimming caused by the exoplanet.

As for the recovery of the injected planet-moon transit, the model parameters are indicated with blue dots in Fig. 5 and the posteriors are shown as gray scale density distributions. As for the confidence intervals of the probability density function, the integral under the 2D normal distribution is given as 1 - e^-(x/σ)2/2, where r is the distance from the mean value. The innermost contour, which refers to one standard deviation (σ), thus contains $1-e^{-2^2/2}=39.3\%$. And, in analogy, the 2σ and 3σ intervals contain 86.4% and 98.9%, respectively. With that being said, the posterior distribution that we find in Fig. 5 is multimodal for most parameter pairs. Hence, the integrated probability contained inside these contours and the formal confidence intervals of our best solution do not refer to normally distributed errors. Some ground truths are formally several σ away from the posterior peaks; however, the maximum a posteriori values are not far off from the ground truth in an absolute sense. For example, the posterior contours of Ω_pm are clearly outside of the truth value, but the error is just ~1% and this is inconsequential. The resulting light curve from the maximum likelihood solution (Fig. 4) is essentially a perfect match to the injected data.

Fig. 4 Simulated planet-moon observations (black points) with a noise amplitude of 100 ppm per data point and a 10 min cadence, roughly corresponding to an idealized observation of an mv ~ 11 photometrically quiet Sun-like star in PLATO’s P5 sample. The system is the same as in Fig. 3 (see Sect. 3.2 for details), (a) Same transit as in Fig. 3d, but with a new realization of the random noise. From our Ul traNest recovery, we selected 100 model parameters from the posteriors and we show their light curves for the moon (orange), planet (blue), and total (gray) light curves.

Fig. 5 Corner plot from the UltraNest retrieval of the system shown in Fig. 3, but with four transit epochs. We note that Rs is given in units of R⊙, ab in units of AU, and Mp in units of Jupiter masses. Due to the small number of transits, the posterior is multimodal for combinations of Ppm and τpm. Additional transits would constrain the system much better.

4 Discussion

To our knowledge, Pandora is the first photodynamical open-source exomoon transit detection algorithm. The LUNA photodynamical modeling code by Kipping (2011) is a proprietary code written in the FORTRAN programming language and it has been used by the Hunt for Exomoons with Kepler survey (Kipping et al. 2012). Our Pandora code is open source, written in the python programming language, and it appears to be about four to five orders of magnitude faster than LUNA. We hope that all of this will make Pandora accessible to a wide range of users and start a community approach in the exomoon search.

Our Pandora code can be used with any nested sampler of choice, such as UltraNest, MultiNest (Feroz et al. 2009), or Dynesty (Speagle 2020). Performance differences between nested samplers, when using the same parameters, are negligible considering a complete run. At the beginning of a sampling run, when the sampling efficiency is high ($\left(\gtrsim 1\%\right)$), a relevant fraction of computing time is spent by the sampler to identify new proposal regions. This applies to the first few percent of a search. Here, MultiNest, written in Fortran, is two times faster than the Python implementations of UltraNest and Dynesty. Later in the run, when the sampling efficiency goes down, this fraction becomes negligible. There may exist novel sampling methods in the newer packages that could converge with fewer steps, but we have not explored them yet.

Other open-source packages such as planetplanet (Luger et al. 2017) and starry (Luger et al. 2019) cover some combination of Kepler solvers and multibody occultations, but they have not been adapted for exomoon purposes. Alternative methods that analyze only the light curve that is not affected by the planetary transit (Simon et al. 2012; Heller 2014; Heller et al. 2016a; Kipping 2021) or TTV-TDV effects (Heller et al. 2016b) are simpler to implement, but also less sensitive (Heller & Hippke 2022).

We are currently improving the fitting and sampling procedure with Pandora, so that priors can be chosen in a more sophisticated way. To fit the quadratic limb darkening coefficients more efficiently, for example, Pandora will offer a conversion routine to calculate $u_1=2\sqrt{q_1}{q}_2$ and $u_2=\sqrt{q_1}\left(1-2q_2\right)$ based on q₁ and q₂ from the unit hypercube. This procedure has been shown to reduce the prior volume (Kipping 2013).

The posterior distribution of our simulated exomoon search with Pandora in Fig. 5 is multimodal and it shows an interesting aliasing effect for the orbital period of the planet-moon system. This effect was first theoretically described by Kipping (2009). It has also been observed before in analyses of the four transits of the giant planet Kepler-1625 b and its Neptune-sized exomoon candidate, namely in Fig. S16 of Teachey & Kipping (2018) and in Fig. 4 of Heller et al. (2019). We expect this to be a general effect of Monte Carlo approaches for exoplanet-exomoon data analysis that includes a small number of transits. In other words, if only a few transits are available, there are multiple combinations of planetary masses, moon periods, and moon orbital positions that result in equally likely solutions. This degeneracy can be increasingly resolved with additional transit epochs. Some of our tests also suggest that eclipses can have an important effect on the resolution of this degeneracy. In extremely close exoplanet-exomoon systems, in which P_pm/2 < d, multiple eclipses can occur and possibly constrain the planet-moon orbital period substantially.

Beyond stellar photometry, stellar radial velocity measurements can help to constrain the total mass M_p + M_m. The simultaneous fitting of transit photometry and stellar radial velocities thus has the potential to aid further in the resolution of the parameter degeneracy.

5 Conclusion

Pandora is the first open source astrophysical model software for stellar light curves with a transiting exoplanet-exomoon system. Pandora has a short code base to increase human readability and facilitate code development, and it is optimized for computational speed. It uses well-established analytical solutions to the transit light curve with stellar limb darkening in most scenarios. In rare cases, if a planet-moon eclipse occurs at the same time when one of the two transiting bodies touches the stellar disk, Pandora calculates the intersection of the star, the planet, and the moon numerically assuming constant stellar brightness behind the small intersecting area. Pandora also offers the option to use the small-body approximation for the moon, which accelerates computations even further. We derived an analytical solution to the transit duration of the planet-moon barycenter in eccentric orbits around the star to avoid computationally expensive solving of Kepler’s equation.

Pandora also supports supersampling for realistic modeling of simulated observations. The output formats include light curves of the moon and the planet, as well as the combined transit model and simulated observations with arbitrary white noise amplitudes. Pandora also comes with a video animation generator that is suitable for testing and debugging, but also for education and outreach.

Finally, we demonstrate Pandora’s potential for exomoon searches by first simulating a physically plausible - though not necessarily representative – system of a Sun-like star, a Jupiter-sized transiting planet at 1AU from the star, and a Neptune-sized moon in a relatively tight orbit of five planetary radii in a 1.28 days period around the planet. We attribute a cadence and noise properties to the simulated model that corresponds to a photometrically inactive, Sun-like star as it could be observed with the 2026 PLATO mission. Our Monte Carlo search for the exomoon with UltraNest takes about 5 h on a standard laptop.

Our nested sampling of simulated exoplanet-exomoon light curves reveals an aliasing effect of the orbital period in planet-moon systems. This is a fundamental aspect of transiting exoplanet-exomoon systems, which Pandora is exquisitely suited for to explore in the future. Our next step is to search for exomoons around all of the thousands of transiting exoplanets and exoplanet candidates from the Kepler mission and to study the most interesting candidates in detail.

Acknowledgements

The authors thank the anonymous referee for a thorough report. RH acknowledges support from the German Aerospace Agency (Deutsches Zentrum für Luft- und Raumfahrt) under PLATO Data Center grant 50OO1501.

Appendix A Error estimates for constant in-transit tangential velocity approximation

The assumption that the in-transit velocity component of the planet-moon barycenter that is tangential to the celestial plane is constant in Eq. (11) provides a substantial computational acceleration of Pandora compared to a dedicated solving algorithm that approximates Kepler’s equation for the eccentric anomaly. Here we estimate the error that is introduced by our approximation.

As shown in Fig. 1, the true anomaly at ingress of the barycenter is (c_b - y), with y = arctan(R_s/r). During egress, the true anomaly is (c_b + y). As a consequence, the velocity component that is tangential to our line of sight during ingress and egress is $$\left|{\upsilon }_{\text{t},\text{in}}\right|=\frac{n_{\text{b}}{a}_{\text{b}}}{\sqrt{1-e_{\text{b}}^2}}(1+e_{\text{b}}\cos ({\varpi }_{\text{b}}-\gamma ))$$(A.1) $$\left|{\upsilon }_{\text{t},\text{eg}}\right|=\frac{n_{\text{b}}{a}_{\text{b}}}{\sqrt{1-e_{\text{b}}^2}}(1+e_{\text{b}}\cos ({\varpi }_{\text{b}}-\gamma )),$$(A.2)

respectively. In Fig. A.1 we compare the constant velocity assumption |v_t| from Eq. (10) with |u_tin| and |ø_t,_eg | from Eqs. (A.1) and (A.2), respectively. We also computed the maximum deviation between the ingress and midtransit tangential velocities, Δi¾ = ¾_n/¾_n - 1, as a function of c_b. For all curves, the stellar radius was set to R_s = R_ø and the stellar mass to M_s = M_ø.

In Fig. A.1(a) we consider a planet-moon barycenter with an orbital period of 30 days, which means that a_b = 0.19 AU. The three curves assume e_b = 0.01, e_b = 0.1, and e_b = 0.2, respectively (see labels). The black dashed lines refer to |v_t|, the orange solid line to |«y_n|, and the light blue solid line to |o_t,_eg|. The lower panel shows the resulting Δv_t values and we find that the maximum value, or error in our assumption of constant in-transit tangential velocity, reaches maximum values of about 5 × 10^-3 for e_b = 0.2, roughly 2.5 × 10^-3 for e_b = 0.1, and approximately 2.5 × 10^-4 for e_b = 0.01. These maximum errors occur for c_b = 90° and c_b = 270°, when the change in the tangential orbital velocity is largest.

In Fig. A.1(b) we consider a planet-moon barycenter with an orbital period of 365.25 days, which means that ab = 1 AU, again for e_b = 0.01, e_b = 0.1, and e_b = 0.2. The resulting maximum error in our assumption of constant in-transit tangential velocity is about 9.7 × 10^-4 for eb = 0.2, roughly 4.7 × 10^-4 for eb = 0.1, and approximately 4.6 × 10^-5 for eb = 0.01.

Fig. A.1 Transit tangential velocity of the planet-moon barycenter during ingress (vt,in), midtransit (vt), and egress (vt,eg) as a function of the orientation of the periastron with respect to the line of sight (cb). (a) The planet-moon barycenter has an orbital period of 30 days around a Sunlike star. The upper panel shows vt,in (red solid line), vt (black dashed line), and vt,eg (blue dashed line) for orbital eccentricities of 0.2, 0.1, and 0.01. The lower panel shows the maximum variation of the transit velocity, which we illustrate as Δvt = ot,in/ot,in - 1 and in units of permille. (b) Same as (a), but with a circumstellar orbital period of 365.25 days. The resulting variation of the in-transit tangential velocity is reduced by almost one order of magnitude, as shown in the lower panel.

We estimated the resulting error in the flux measurements for the transit light curve of a Jupiter-sized planet (R_p = 0.1 R_s) in the worst-case scenario for these above cases, in which P_b = 30 d, e_b = 0.2, and c_b = 90°. We used the batman package and its internal Kepler orbit solver to generate a highly accurate model light curve. This approach takes the variable orbital velocity during transit due to eccentricity into account. We calculated 10,000 in-transit data points, each of which is based on the accurately calculated orbital position of the planet as a function of time. For comparison, we calculated a similar light curve with batman for e = 0 and then stretched this light curve in time via multiplication with a factor s = (1 - eb)^1/2/(1 + e_b cos(c_b)) as per Eq. (11) to produce the transit duration that corresponds to our assumption of constant in-transit velocity. In this particular case, we obtained s = 0.9797959.

Fig. A.2 shows the resulting flux error as a function of time, with the origin of the time axis set to the middle of the transit. The maximum flux difference is about 4 × 10^-7 during ingress. We also tested the error in our approximation for different planetary radii and noticed that the error scales roughly linearly with R_p. For Earth-sized transiting objects, the error is a few times 10^-8.

Fig. A.2 Example of the flux error resulting from our approximation for the transit duration in an eccentric (e = 0.2) circumstellar orbit. We assume a worst-case orientation of periastron $\left({\varpi }_{\text{b}}={90}^{\circ }\right)$ and a Jupiter-sized planet to exaggerate the error. The error was computed as the difference between two light curves computed with the batman package. One light curve was initially simulated as a circular orbit and then stretched with a factor of $s={\left(1-e_{\text{b}}^2\right)}^{1/2}/\left(1+e_{\text{b}}\cos \left({\varpi }_{\text{b}}\right)\right)$ as per Eq. (11). The other light curve was computed with e = 0.2 and using a computationally much more demanding Kepler solver for the orbital motion of the planet.

Appendix B Error estimates for numerical modeling of planet-moon eclipses

If a planet-moon eclipse occurs in combination with an ingress or egress of either the moon or the planet, then Pandora switches into a numerical modeling mode, that is, to eclipse case (2) (Sect. 2.5.2). The moon is modeled as a pixelated area with a diameter of n pixels. The circumference of the moon then is πn and the error of the approximated moon area is on the order of one pixel along its circumference, or (πn)^-1. For n = 25, as in the default Pandora mode, the resulting error for the moon area is ~1.3 %. As an example, an Earth-sized moon transiting a Sun-like star causes a flux drop of ~84 ppm. The numerical error in the light curve is then ~1 ppm and therefore more than an order of magnitude below the astrophysical variability of a photometrically quiet star like the Sun on a time scale of minutes to hours.

In addition to these analytical estimates, we undertook a detailed suite of numerical simulations with Pandora and compared the area of a hypothetical Earth-sized moon from the pixelated area in Pandora’s eclipse case (2) to the actual moon area. The resulting error and run time per data point are shown in Fig. B.1 The blue line illustrates that the error in the light curve is 1 ppm for a moon diameter of 25 pixels, which is in perfect agreement with our theoretical estimates. The orange line denotes that the resulting computer run time is 7.4µs per data point.

Fig. B.1 Error estimates (blue line, left ordinate) and CPU run time per data point (orange line, right ordinate) as a function of the moon diameter in Pandora’s numerical eclipse case (2) (see Sect. 2.5.2). The error refers to the absolute error in the flux of the transit light curve of an Earth-sized body in front of a Sun-like star. The error falls below 1 ppm for grid sizes > 25 px, which is Pandora’s default mode. Users of Pandora can choose suitable grad sizes depending on accuracy requirements.

As for the frequency and duration of planet-moon eclipses, they are rare in the Solar System and they can be expected to be rare in most exomoon scenarios. The highest frequency and longest duration occurs in an edge-on geometry, when the planet–moon orbit is in alignment with the line of sight. In this case, the maximum fraction of the planet-moon orbit spent in eclipse is roughly 2R_p/(πa_pm). For the Earth–Moon system, this corresponds to about 1% of the time; also, for the Galilean moons around Jupiter, this corresponds to about 10.4% (Io), 6.6% (Europa), 4.1% (Ganymede), and 2.3% (Callisto), respectively.

The fraction of planet-moon eclipses which happens during ingress or egress, that is, Pandora’s eclipse case (2) that we solved numerically, is even smaller. For an Earth transit across the diameter of the Sun, ingress and egress account for 1.8% of the total transit duration. For Jupiter, it is 20%. Consequently, the average time for which a planet-moon occultation occurs during ingress or egress is only 0.018% of the transit data points for an Earth-Moon system, or 2.1% (Io), 1.3% (Europa), 0.82% (Ganymede), and 0.46% (Callisto) for the Galilean moons, on average.

Appendix C Algorithmic insights for higher speed

For computational speed optimization of Pandora, it has proven crucial to break the code down into modules that can be bench-marked separately. It is often hard to estimate whether a change to the code makes it slower or faster, but speed can – and in our opinion should – be measured.

C.1 Technical implementation

All modules were implemented in pure python, so Pandora does not require any compilation of C or Fortran code. Pandora also uses just-in-time (JlT) compilation with numba, which reduces Pandora’s run time by a factor of about 1100. All numbers with respect to the performance of Pandora given in this paper are numba-based.

Fig. C.1 Flux errors resulting from the small-body approximation with linear interpolation. Left: For Rp/Rs = 0.05, the maximum distance between the center of the stellar disk and the occulting body that provides errors < 1 ppm is ɀ = 0.65. Right: For Rp/Rs = 0.02, flux errors are < 1 ppm up to z = 0.95. Black curves are realizations for 1400 optical limb darkening coefficients of all stellar types (Claret 2017).

To further assess the performance of our implementation, we rewrote our algorithm for the Kepler ellipse in the circular case in the C language. The pure-C version, compiled with GCC-11.1, is equally fast within our measurement errors of a few percent. Similarly, we benchmarked our occultation code against batman (C language), ALFM (Agol et al. 2020, C++ and Julia), and PyTransit (also in python-numba). Depending on the parameter space, we find Pandora to be 5–20% faster. In summary, the attractiveness of Pandora is in our combination of JIT compilation for optimal speed and python as a programming language to ensure readability by a wide range of users.

Pandora first determines the spatial coordinates for all points, and afterwards all occultation fluxes, etc. This allows each module to iterate over all values sequentially, making full use of a CPU core’s pipeline and cache, resulting in a ~ 10× performance gain. We have also optimized Pandora to make use of instruction-level parallelism in modern superscalar processors. In each clock cycle, multiple instructions can be executed, increasing the throughput by a factor of a few. A key requirement to achieve this is the reduction of conditional (if) statements.

C.2 Conditional small-body approximation with linear interpolation

The small-body approximation assumes constant limb darkening under the occulted area when calculating fluxes. It is ~ 16× faster to calculate than the exact Mandel-Agol equations. Errors from this approximation are most pronounced during ingress and egress, and they become very small (< ppm) for bodies R_p/R_s < 0.01. Larger bodies suffer from relevant errors during most of the transit, and the size of this error is dependent on the limb darkening coefficients. To reduce this error, we performed a linear interpolation between the exact and the approximated functions. This requires only the calculation of two exact values, which takes a negligible amount of time, and allows one to determine most of the occultation light curves with the small-body approximation. For R_p/R_s = 0.05, the approximation can be used for 65% of the data points in a light curve, growing to 95% for R_p/R_s = 0.02, keeping small < 1 ppm errors Fig. C.1. For bodies R_p/R_s > 0.05, the method is not useful because errors become significant.

Fig. C.2 Performance comparison of the four methods provided by Pandora to calculate occultation fluxes. The Mandel-Agol model (black) is the most precise, but also slowest, reference version. The hybrid model (red) uses interpolation between the small-body and the Mandel-Agol methods, and it is applicable for RP/RS < 0.05 to keep errors below 1 ppm. The cached algorithm (blue) is still faster, but only usable for sampling with fixed limb darkening. Finally, the small-body approximation is the fastest method, but only applicable (errors below 1 ppm) for RP/RS < 0.01. All methods show better performance for longer light curves due to constant factors. Numbers are given for a single core on an Intel Core-i7 1185G.

The interpolation method requires one to calculate the exact flux F_x(ɀ, k, u₁, u₂) for only two radial distances ɀ = 0 and ɀ = 0.65, the latter being determined empirically (k is the occulter radius in units of R_p/R_s, and u₁ and u₂ are the quadratic limb darkening coefficients). We then calculated all other values ɀ < ɀ_cutoff with the small-body approximation and subtracted the difference between these values and the linear interpolation through F₁ and F₂ to reduce the errors. We chose the cutoff so that errors remained < 1 ppm, allowing for ɀ < 0.65 (k = 0.05), ɀ < 0.70 (k = 0.04), ɀ < 0.80 (k = 0.03), ɀ < 0.95 (k = 0.02), and ɀ < 0.98 (k = 0.01). We hard-coded these values in Pandora. The computational performance of these methods, as a function of the number of datapoints, is shown in Fig. C.2.

C.3 Substitution of trigonometric functions

Through extensive profiling, it became clear that the major computational cost (90%) for the Kepler ellipse calculation was originally in the calculation of trigonometric functions. Standard implementations such as the one in PyAstronomy calculate a range of sin, cos, tan, and arctan terms for each resulting data point, where the 3D (x, y, z) positions are determined via $$k=\pi \left(\text{time}-\tau \right)/\text{per}$$(C.1) $$Q=2\text{arctan}\left(\tan /\left(k\right)\right)$$(C.2) $$V=\sin \left(Q\right)\cos \left(i\right)$$(C.3) $$x=\left(\cos \left(\text{Ω}\right)\cos \left(\text{Ω}\right)-\sin \left(\text{-c3Ω}\right)V\right)a$$(C.4) $$y=\left(\sin \left(\text{Ω}\right)\cos \left(\text{Ω}\right)-\cos \left(\text{Ω}\right)V\right)a$$(C.5) $$z=\sin \left(Q\right)\sin \left(i\right)a$$(C.6)

First, we can substitute Q = arctan(tan(k)) = k. Then, profiling shows that almost all time is spent on the expressions sin(2k) and cos(2k) (where $k=\pi \left(\text{time-}\tau \right)/\text{per})$, which is trivially fast to calculate). As trigonometric functions are very expensive, we can substitute the sine and cosine calculations for one tan calculation through the following identity: $$\cos \left(2k\right)=\frac{1-{\tan }^2\left(k\right)}{1+{\tan }^2\left(k\right)}$$(C.7)

and $$\sin \left(2k\right)=\frac{2\tan \left(k\right)}{1+{\tan }^2\left(k\right).}$$(C.8)

On standard CPU architectures (Intel Core i7 and Apple M1), the substitution is ~ 10% faster. Additional trigonometry is still formally required to determine the Kepler ellipse, but only once, and not for each of the ~ 1000 data points on the ellipse orbit. That is to say, each data point can be calculated with only one trigonometric calculation of tan(k). Finally, we can remove the calculation of the depth component ɀ because the orientation of the planet-moon system is irrelevant as both are black disks.

C.4 Numerical accuracy requirements

In computer hardware, floating-point numbers are represented as base-2 (or “binary”) fractions. With a finite number of bits available for a given number, any number is an approximation. High-level programming languages such as python use double-precision floating point data types with 64-bit accuracy, resulting in almost 16 significant decimal digits of precision. The standard math library follows the IEEE-754 arithmetic, which aims to provide accuracy to the last significant digit, that is, an error < 0.5 in the unit in the last place (ULP). As errors can add up over the course of multiple calculations, it is useful to test the requirements for a given problem. With that being said, it appears to us that calculating every value in Pandora to 16 digits would be excessive..

In the case of the Kepler ellipse, we calculated a tangent for each point on the orbit and then multiplied this value with various other parameters. We required an accuracy well below the noise present in the data, which is typically 10^-5. The signal of Earth’s Moon transiting a Sun-like star would be about 6 ppm (6 × 10^-6), which implies a precision requirement of at least seven decimal digits to achieve a model that is substantially more accurate than the signal. This requirement is about nine orders of magnitude weaker than the 16 digits typically available in our python implementation of Pandora.

There are two ways to reduce computational effort at the cost of accuracy: Single precision floats and/or a different trigonometric algorithm. The numba package is built upon the LLVM compiler stack, which offers a set of fastmath approximations, resulting in up to four ULPs of errors. The speed gain with this option is 20–50%. Single precision floats provide 7.225 decimal digits of precision on average, which is insufficient given our requirement of seven digits and multiple calculations on intermediate values. The speed difference would be a factor of two on most architectures. Working with double precision (16 digits), we can accept many ULPs of fast approximation errors. Consequently, Pandora uses 64-bit arithmetic, but low accuracy (four ULPs) approximations.

C.5 Potential for future algorithmic optimization

Various alternative algorithms exist which could further improve Pandora’s performance. We have deferred the implementation of these for future releases, given sufficient interest by the community.

As an example, our Kepler solver for the eccentric case uses the method by Markley (1995). Faster algorithms exist, for example one is based on a CORDIC-like solver (Zechmeister 2018, 2021). As eccentric moons are rarely expected, we have deferred to these optimizations.

The occultation algorithm ALFM (Agol et al. 2020) uses a different method to calculate the circle-circle intersect than the one used by Pandora. The standard method of computing the overlap of two circles requires the calculation of a square-root and of two arccos terms. The new “kite expression” replaces the a ccos with inverse a ctan, keeping the same constant factors, while achieving higher numerical accuracy. Our benchmarks in numba and C-language yield identical performance within a measurement uncertainty of 1%. Due to the lower complexity of the standard approach, we decided to keep the classical circle-circle intersect formula in Pandora for now. The improved accuracy of the kite expression is irrelevant for our purpose, as the errors from the standard method are < 10^-8 in all cases (Fig. 2 in Agol et al. 2020).

The occultation part of Pandora also contains code to calculate the complete elliptical integral of the third kind. A well-known standard method for this calculation, which is also used in EXOFAST (Eastman et al. 2013), is the iterative approximation by Bulirsch (1965a,b). As an alternative, Fukushima (2013) describes a method based on half and double argument transformations, which is benchmarked by the author to be 20–50% faster than Bulirsch’s method. We are grateful to Toshio Fukushima for sharing his Fortran code (~ 700 lines of code, including precalculated constants) to perform a comparison. In our use case, however, our numba version of Bulirsch’s method is slightly faster.

C.6 GPU-based calculations

We have implemented a CUDA version of Pandora to run on Nvidia graphics cards (GPUs). A performance test on a Nvidia GTX 1080 card, with 3584 Cuda-cores, yields a sustained throughput of 225,000 model calculations per second, including the corresponding log-likelihood evaluations. This is a speed-up by a factor of 20 compared to a single core of an Intel Core i7-7700k processor. Unfortunately, it appears that current versions of nested samplers such as UltraNest cannot process this load sufficiently fast and that they are limited to ~ 50,000 point proposals per second due to internal overhead. We plan to address these issues in a subsequent version of Pandora because a 20-fold performance gain would reduce typical convergence times from five hours to fifteen minutes.

Appendix D Usage examples

D.0.1 Defining model parameters

A full list of model parameters is available in the online documentation.

D.0.2 Creating a light curve

D.0.3 Obtaining coordinates

The return values px, py determine the x and y components of the planet, while mx and my are for the moon; they were sampled at each time stamp.

D.0.4 Generating a transit animation video

References

All Tables

All Figures

	Fig. 1 Orbit geometry of the planet-moon barycenter (Mb) on its eccentric circumstellar orbit as used in Pandora. The situation depicted in this illustration corresponds to the midtransit time, when the true anomaly coincides with the orientation of the periastron (ωb). Nothing is to scale in this illustration. The bottom part of the ellipse has been spared for illustration purposes.
In the text

Fig. 2 Planet-moon eclipse in contact with the stellar limb. (a) The planet (blue circle) and the moon (orange circle) are both in front of the star (black circle segment) and in eclipse. (b) Numerical pixel grid representing the moon. Different areas of overlap are shaded in arbitrary colors. In this specific configuration, the white area in the upper left region of the image is determined by pixel area summation and then compensated for in the light curve. The right half of the moon is not transiting.

In the text

Fig. 3 Output demo of Pandora for a system of a Sun-like star, a Jupiter-sized planet in a one-year orbit around that star, and a Neptune-sized moon in a 1.28 days orbit around the giant planet (for details see Sect. 3.2). (a)–(c) Video renderings of the transit. (d) Light curve. The dotted-dashed orange line shows the transit light curve of the exomoon and the dashed blue line shows the transit light curve of the exoplanet using the analytical solution in both cases. The yellow parts of the moon and planet light curves illustrate numerical simulations. The solid black line shows the combined model. Black dots show a simulated observation roughly representative of an mV ~ 11 Sun-like star from the PLATO mission. Digital star colors are from Harre & Heller (2021).

In the text

Fig. 4 Simulated planet-moon observations (black points) with a noise amplitude of 100 ppm per data point and a 10 min cadence, roughly corresponding to an idealized observation of an mv ~ 11 photometrically quiet Sun-like star in PLATO’s P5 sample. The system is the same as in Fig. 3 (see Sect. 3.2 for details), (a) Same transit as in Fig. 3d, but with a new realization of the random noise. From our Ul traNest recovery, we selected 100 model parameters from the posteriors and we show their light curves for the moon (orange), planet (blue), and total (gray) light curves.

In the text

	Fig. 5 Corner plot from the UltraNest retrieval of the system shown in Fig. 3, but with four transit epochs. We note that Rs is given in units of R⊙, ab in units of AU, and Mp in units of Jupiter masses. Due to the small number of transits, the posterior is multimodal for combinations of Ppm and τpm. Additional transits would constrain the system much better.
In the text

Fig. A.1 Transit tangential velocity of the planet-moon barycenter during ingress (vt,in), midtransit (vt), and egress (vt,eg) as a function of the orientation of the periastron with respect to the line of sight (cb). (a) The planet-moon barycenter has an orbital period of 30 days around a Sunlike star. The upper panel shows vt,in (red solid line), vt (black dashed line), and vt,eg (blue dashed line) for orbital eccentricities of 0.2, 0.1, and 0.01. The lower panel shows the maximum variation of the transit velocity, which we illustrate as Δvt = ot,in/ot,in - 1 and in units of permille. (b) Same as (a), but with a circumstellar orbital period of 365.25 days. The resulting variation of the in-transit tangential velocity is reduced by almost one order of magnitude, as shown in the lower panel.

In the text

Fig. A.2 Example of the flux error resulting from our approximation for the transit duration in an eccentric (e = 0.2) circumstellar orbit. We assume a worst-case orientation of periastron $\left({\varpi }_{\text{b}}={90}^{\circ }\right)$ and a Jupiter-sized planet to exaggerate the error. The error was computed as the difference between two light curves computed with the batman package. One light curve was initially simulated as a circular orbit and then stretched with a factor of $s={\left(1-e_{\text{b}}^2\right)}^{1/2}/\left(1+e_{\text{b}}\cos \left({\varpi }_{\text{b}}\right)\right)$ as per Eq. (11). The other light curve was computed with e = 0.2 and using a computationally much more demanding Kepler solver for the orbital motion of the planet.

In the text

Fig. B.1 Error estimates (blue line, left ordinate) and CPU run time per data point (orange line, right ordinate) as a function of the moon diameter in Pandora’s numerical eclipse case (2) (see Sect. 2.5.2). The error refers to the absolute error in the flux of the transit light curve of an Earth-sized body in front of a Sun-like star. The error falls below 1 ppm for grid sizes > 25 px, which is Pandora’s default mode. Users of Pandora can choose suitable grad sizes depending on accuracy requirements.

In the text

Fig. C.1 Flux errors resulting from the small-body approximation with linear interpolation. Left: For Rp/Rs = 0.05, the maximum distance between the center of the stellar disk and the occulting body that provides errors < 1 ppm is ɀ = 0.65. Right: For Rp/Rs = 0.02, flux errors are < 1 ppm up to z = 0.95. Black curves are realizations for 1400 optical limb darkening coefficients of all stellar types (Claret 2017).

In the text

Fig. C.2 Performance comparison of the four methods provided by Pandora to calculate occultation fluxes. The Mandel-Agol model (black) is the most precise, but also slowest, reference version. The hybrid model (red) uses interpolation between the small-body and the Mandel-Agol methods, and it is applicable for RP/RS < 0.05 to keep errors below 1 ppm. The cached algorithm (blue) is still faster, but only usable for sampling with fixed limb darkening. Finally, the small-body approximation is the fastest method, but only applicable (errors below 1 ppm) for RP/RS < 0.01. All methods show better performance for longer light curves due to constant factors. Numbers are given for a single core on an Intel Core-i7 1185G.

In the text