Issue 
A&A
Volume 662, June 2022



Article Number  A37  
Number of page(s)  16  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/202243129  
Published online  14 June 2022 
Pandora: A fast opensource exomoon transit detection algorithm
^{1}
Sonneberg Observatory,
Sternwartestraße 32,
96515
Sonneberg,
Germany
email: michael@hippke.org
^{2}
Visiting Scholar, Breakthrough Listen Group, Berkeley SETI Research Center, Astronomy Department,
UC Berkeley
USA
^{3}
MaxPlanckInstitut für Sonnensystemforschung,
JustusvonLiebigWeg 3,
37077
Göttingen,
Germany
email: heller@mps.mpg.de
^{4}
Institut für Astrophysik, GeorgAugustUniversität Göttingen,
FriedrichHundPlatz 1,
37077
Göttingen,
Germany
Received:
17
January
2022
Accepted:
22
February
2022
We present Pandora, a new software to model, detect, and characterize transits of extrasolar planets with moons in stellar photometric time series. Pandora uses an analytical description of the transit light curve for both the planet and the moon in front of a star with atmospheric limb darkening and it covers all cases of mutual planet–moon eclipses during transit. The orbital motion of the starplanetmoon system is computed with a high accuracy as a nested Keplerian problem. We have optimized Pandora for computational speed to make it suitable for largescale exomoon searches in the new era of spacebased highaccuracy surveys. We demonstrate the usability of Pandora for exomoon searches by first simulating a light curve with four transits of a hypothetical Jupiter with a giant Neptunesized exomoon in a oneyear orbit around a Sunlike star. The 10 min cadence of the data matches that of the upcoming PLATO mission and the noise of 100 parts per million is dominated by photon noise, assuming a photometrically quiet, m_{V} = 11 Sunlike star for practicality. We recovered the simulated system parameters with the UltraNest Bayesian inference package. The runtime of this search is about five hours on a standard computer. Pandora is the first photodynamical opensource exomoon transit detection algorithm, implemented fully in the python programming language and available for the community to join the search for exomoons.
Key words: methods: data analysis / planets and satellites: detection / techniques: photometric
© 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 exomoons 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 planetmoon eclipses, and potentially due to the small physical radii of exomoons, all of which make phasefolding 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 planetmoon 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 planetmoon 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 exoplanetexomoon 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 planetplanet mutual eclipses. Given the recent announcements of transiting exomoon candidates around Kepler1625 (Teachey et al. 2018) and Kepler1708 (Kipping et al. 2022) and the lack of available computer code to investigate these and future exomoon claims independently, an opensource exomoon transit detection algorithm is very desirable for the community.
Here, we present Pandora, the first opensource exomoon transit detection algorithm^{1} using a full photodynamical model. This stateoftheart algorithm for exomoon transit modeling and searches requires just a few hours of CPU time for a Keplerlike 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 justintime 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_{0}), the orbital period of the barycenter around the star (P_{b}), the semimajor 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 edgeon perspective. For a summary of all parameters, see Table 1.
Moving on to the parameterization of the planetmoon 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 semimajor 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 periapsis 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_{1}u_{2}) 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 spacebased highaccuracy 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 where G is the gravitational constant.
2.2 Circumstellar orbital eccentricity
Our threebody 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 xaxis, in which case the orbital velocity of M_{b} is given as (Murray & Dermott 1999)
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 skyprojected 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: (2) (3)
The projection of u on e_{r}, which we refer to as υ_{r}, is given as (4)
where a = <(υ, e_{r}) (Fig. 1). The cosine of this angle between υ and e_{r} can be calculated as (5)
Plugging Eq. (5) into Eq. (4) yields (6)
Analogously, projection of u on e_{t}, which we refer to as υ_{t}, yields (7)
Plugging Eq. (8) into Eq. (7) yields (9)
For our purpose of using the midtransit tangential velocity, we set f = ω_{b} and found (10)
and using as the transit path across the star, we constructed the transit duration of a pointlike M_{b} as (11)
In Appendix A we extend these considerations to the tangential transit velocity during ingress and egress and demonstrate that our assumption of constant intransit velocity of the planetmoon 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 planetmoon barycenter (M_{b}) 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 largescale 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 fifthorder 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 .
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 planetmoon 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, 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 planetmoon 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 i71185G7 processor.
Model parameters.
2.4 Stellar flux and limb darkening
To calculate the apparent stellar flux during transits of the planetmoon 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 occultation model, where the smallbody approximation (Mandel & Agol 2002) where constant stellar limb darkening behind the occulted area is assumed whenever appropriate. The error of the smallbody 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 smallplanet approximation. The resulting model fluxes are accurate to <1 ppm in all cases (see Appendix C.2).
In case of fixed limbdarkening 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 planettostar 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 i71185G7 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 planetonly 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 Planetmoon 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, planetmoon 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 starplanetmoon 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 planetmoon eclipse and still covers all possible cases of previous studies.
2.5.1 Planetmoon 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 (12)
while the distance between their centers (d_{pm}) is smaller than the sum of their radii, (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 Planetmoon 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 planetmoon 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 outofeclipse throughput of between 16 and 200 million data points per second and per core (Sect. 2.4). However, since planetmoon eclipses are only relevant for planetmoon 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 nearlinear 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 classbased wrapper function (56 SLOCs) and a video generator (83 SLOCs). Pandora only depends on the numpy library (Harris et al. 2020) and the numba justintime 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 singlebody occultations with limb darkening and supersampling against PyTransit and batman, and the planetmoon 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^{2}. 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 Planetmoon 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 exoplanetexomoon system modeled with Pandora. The star is assumed to have a solar radius and solartype limb darkening, with u_{1} = 0.4089 and u_{2} = 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 planetmoon 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 planetmoon barycenter.
This period, which implies an orbital semimajor 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 fluidlike objects. For comparison, Io, the innermost of the Galilean moons around those of Jupiter, is at roughly 6.1 R_{p} around Jupiter. The planetmoon 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 wellsuited 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, Sunlike 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 F5K7 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 Sunlike star from PLATO’s P5 sample^{3}, 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 planetmoon barycenter, at which point the planet and the moon exhibit nearmaximum 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 (dotteddashed 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 Sunlike star, a Jupitersized planet in a oneyear orbit around that star, and a Neptunesized 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 dotteddashed 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 m_{V} ~ 11 Sunlike 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 starplanetmoon 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 generalpurpose Bayesian inference package for parameter estimation and model comparison. It allows one to fit arbitrary models for multimodal or nonGaussian parameter spaces as can be expected for exoplanetexomoon systems due to alias effects in the orbital period of the planetmoon 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 i71185G7 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 loglikelihood 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 planetmoon 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 . 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 planetmoon 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 m_{v} ~ 11 photometrically quiet Sunlike 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 R_{s} is given in units of R_{⊙}, a_{b} in units of AU, and M_{p} in units of Jupiter masses. Due to the small number of transits, the posterior is multimodal for combinations of P_{pm} and τ_{pm}. Additional transits would constrain the system much better. 
4 Discussion
To our knowledge, Pandora is the first photodynamical opensource 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 (), 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 opensource 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 TTVTDV 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 and based on q_{1} and q_{2} 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 planetmoon 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 Kepler1625 b and its Neptunesized 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 exoplanetexomoon 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 exoplanetexomoon systems, in which P_{pm}/2 < d, multiple eclipses can occur and possibly constrain the planetmoon 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 exoplanetexomoon system. Pandora has a short code base to increase human readability and facilitate code development, and it is optimized for computational speed. It uses wellestablished analytical solutions to the transit light curve with stellar limb darkening in most scenarios. In rare cases, if a planetmoon 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 smallbody approximation for the moon, which accelerates computations even further. We derived an analytical solution to the transit duration of the planetmoon 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 Sunlike star, a Jupitersized transiting planet at 1AU from the star, and a Neptunesized 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, Sunlike 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 exoplanetexomoon light curves reveals an aliasing effect of the orbital period in planetmoon systems. This is a fundamental aspect of transiting exoplanetexomoon 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 intransit tangential velocity approximation
The assumption that the intransit velocity component of the planetmoon 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 (A.1) (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 planetmoon 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 intransit 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 planetmoon 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 intransit 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 planetmoon barycenter during ingress (v_{t},_{in}), midtransit (v_{t}), and egress (v_{t},_{eg}) as a function of the orientation of the periastron with respect to the line of sight (c_{b}). (a) The planetmoon barycenter has an orbital period of 30 days around a Sunlike star. The upper panel shows v_{t,in} (red solid line), v_{t} (black dashed line), and v_{t,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 Δv_{t} = o_{t},_{in}/o_{t},_{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 intransit 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 Jupitersized planet (R_{p} = 0.1 R_{s}) in the worstcase 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 intransit 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 intransit 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 Earthsized 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 worstcase orientation of periastron and a Jupitersized 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 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 planetmoon eclipses
If a planetmoon 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 Earthsized moon transiting a Sunlike 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 Earthsized 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 Earthsized body in front of a Sunlike 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 planetmoon 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 edgeon geometry, when the planet–moon orbit is in alignment with the line of sight. In this case, the maximum fraction of the planetmoon 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 planetmoon 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 planetmoon occultation occurs during ingress or egress is only 0.018% of the transit data points for an EarthMoon 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 benchmarked 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 justintime (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 numbabased.
Fig. C.1 Flux errors resulting from the smallbody approximation with linear interpolation. Left: For R_{p}/R_{s} = 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 R_{p}/R_{s} = 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 pureC version, compiled with GCC11.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 pythonnumba). 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 instructionlevel 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 smallbody approximation with linear interpolation
The smallbody approximation assumes constant limb darkening under the occulted area when calculating fluxes. It is ~ 16× faster to calculate than the exact MandelAgol 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 smallbody 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 MandelAgol model (black) is the most precise, but also slowest, reference version. The hybrid model (red) uses interpolation between the smallbody and the MandelAgol methods, and it is applicable for R_{P}/R_{S} < 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 smallbody approximation is the fastest method, but only applicable (errors below 1 ppm) for R_{P}/R_{S} < 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 Corei7 1185G. 
The interpolation method requires one to calculate the exact flux F_{x}(ɀ, k, u_{1}, u_{2}) 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_{1} and u_{2} are the quadratic limb darkening coefficients). We then calculated all other values ɀ < ɀ_{cutoff} with the smallbody approximation and subtracted the difference between these values and the linear interpolation through F_{1} and F_{2} 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 hardcoded 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 (C.1) (C.2) (C.3) (C.4) (C.5) (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 , 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: (C.7)
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 planetmoon system is irrelevant as both are black disks.
C.4 Numerical accuracy requirements
In computer hardware, floatingpoint numbers are represented as base2 (or “binary”) fractions. With a finite number of bits available for a given number, any number is an approximation. Highlevel programming languages such as python use doubleprecision floating point data types with 64bit accuracy, resulting in almost 16 significant decimal digits of precision. The standard math library follows the IEEE754 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 Sunlike 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 64bit 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 CORDIClike 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 circlecircle intersect than the one used by Pandora. The standard method of computing the overlap of two circles requires the calculation of a squareroot 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 Clanguage yield identical performance within a measurement uncertainty of 1%. Due to the lower complexity of the standard approach, we decided to keep the classical circlecircle 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 wellknown 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 GPUbased 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 Cudacores, yields a sustained throughput of 225,000 model calculations per second, including the corresponding loglikelihood evaluations. This is a speedup by a factor of 20 compared to a single core of an Intel Core i77700k 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 20fold 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
 Agol, E., Luger, R., & ForemanMackey, D. 2020, AJ, 159, 123 [Google Scholar]
 Batalha, N.M., Rowe, J.F., Bryson, S.T., et al. 2013, ApJS, 204, 24 [NASA ADS] [CrossRef] [Google Scholar]
 Buchner, J. 2021, J. Open Source Softw., 6, 3001 [CrossRef] [Google Scholar]
 Bulirsch, R. 1965a, Numer. Math., 7, 78 [CrossRef] [Google Scholar]
 Bulirsch, R. 1965b, Numer. Math., 7, 353 [CrossRef] [Google Scholar]
 Claret, A. 2017, A&A, 600, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Claret, A., & Bloemen, S. 2011, A&A, 529, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Czesla, S., Schröter, S., Schneider, C.P., et al. 2019, PyA: Python astronomyrelated packages [Google Scholar]
 Eastman, J., Gaudi, B.S., & Agol, E. 2013, PASP, 125, 83 [NASA ADS] [CrossRef] [Google Scholar]
 Feroz, F., Hobson, M.P., & Bridges, M. 2009, MNRAS, 398, 1601 [NASA ADS] [CrossRef] [Google Scholar]
 Fewell, M.P. 2006, Maritime Operations Division, DSTOTN0722, Department of Defence, Tech. Rep. [Google Scholar]
 Fukushima, T. 2013, J. Comput. Appl. Math. 253, 142 [CrossRef] [Google Scholar]
 Harre, J.V., & Heller, R. 2021, Astron. Nachr., 342, 578 [NASA ADS] [CrossRef] [Google Scholar]
 Harris, C.R., Millman, K.J., van der Walt, S.J., et al. 2020, Nature, 585, 357 [NASA ADS] [CrossRef] [Google Scholar]
 Heller, R. 2014, ApJ, 787, 14 [CrossRef] [Google Scholar]
 Heller, R., & Albrecht, S. 2014, ApJ, 796, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Heller, R. & Hippke, M. 2022, A&A, 657, A119 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Heller, R., Hippke, M., & Jackson, B. 2016a, ApJ, 820, 88 [NASA ADS] [CrossRef] [Google Scholar]
 Heller, R., Hippke, M., Placek, B., Angerhausen, D., & Agol, E. 2016b, A&A, 591, A67 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Heller, R., Rodenbeck, K., & Bruno, G. 2019, A&A, 624, A95 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Heller, R., Harre, J.V., & Samadi, R. 2021, A&A, submitted [Google Scholar]
 Hippke, M., & Heller, R. 2019, A&A, 623, A39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hippke, M., David, T.J., Mulders, G.D., & Heller, R. 2019, AJ, 158, 143 [NASA ADS] [CrossRef] [Google Scholar]
 Kipping, D.M. 2009, MNRAS, 392, 181 [NASA ADS] [CrossRef] [Google Scholar]
 Kipping, D.M. 2010, MNRAS, 408, 1758 [NASA ADS] [CrossRef] [Google Scholar]
 Kipping, D.M. 2011, MNRAS, 416, 689 [NASA ADS] [Google Scholar]
 Kipping, D.M. 2013, MNRAS, 435, 2152 [NASA ADS] [CrossRef] [Google Scholar]
 Kipping, D. 2021, MNRAS, 507, 4120 [NASA ADS] [CrossRef] [Google Scholar]
 Kipping, D.M., Bakos, G.Ä., Buchhave, L., Nesvorny, D., & Schmitt, A. 2012, ApJ, 750, 115 [NASA ADS] [CrossRef] [Google Scholar]
 Kipping, D.M., Forgan, D., Hartman, J., et al. 2013, ApJ, 777, 134 [NASA ADS] [CrossRef] [Google Scholar]
 Kipping, D.M., Schmitt, A.R., Huang, X., et al. 2015, ApJ, 813, 14 [NASA ADS] [CrossRef] [Google Scholar]
 Kipping, D., Bryson, S., Burke, C., et al. 2022, Nat. Astron., 6, 367 [NASA ADS] [CrossRef] [Google Scholar]
 Kovâcs, G., Zucker, S., & Mazeh, T. 2002, A&A, 391, 369 [CrossRef] [EDP Sciences] [Google Scholar]
 Kreidberg, L. 2015, PASP, 127, 1161 [Google Scholar]
 Lam, S.K., Pitrou, A., & Seibert, S. 2015, in ACM Digital Library (ACM Press) [Google Scholar]
 Luger, R., LustigYaeger, J., & Agol, E. 2017, ApJ, 851, 94 [NASA ADS] [CrossRef] [Google Scholar]
 Luger, R., Agol, E., ForemanMackey, D., et al. 2019, AJ, 157, 64 [Google Scholar]
 Mandel, K., & Agol, E. 2002, ApJ, 580, L171 [Google Scholar]
 Manduca, A., Bell, R.A., & Gustafsson, B. 1977, A&A, 61, 809 [NASA ADS] [Google Scholar]
 Markley, F.L. 1995, Celest. Mech. Dyn. Astron., 63, 101 [NASA ADS] [CrossRef] [Google Scholar]
 Murray, C.D., & Dermott, S.F. 1999, Solar System Dynamics (Cambridge, UK: Cambridge University Press) [Google Scholar]
 Parviainen, H. 2015, MNRAS, 450, 3233 [Google Scholar]
 Rauer, H., Catala, C., Aerts, C., et al. 2014, Exp. Astron., 38, 249 [Google Scholar]
 Simon, A.E., Szabö, G.M., Kiss, L.L., & Szatmâry, K. 2012, MNRAS, 419, 164 [NASA ADS] [CrossRef] [Google Scholar]
 Speagle, J.S. 2020, MNRAS, 493, 3132 [NASA ADS] [CrossRef] [Google Scholar]
 Teachey, A., & Kipping, D.M. 2018, Sci. Adv., 4, eaav1784 [NASA ADS] [CrossRef] [Google Scholar]
 Teachey, A., Kipping, D.M., & Schmitt, A.R. 2018, AJ, 155, 36 [Google Scholar]
 Weidner, C., & Horne, K. 2010, A&A, 521, A76 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Zechmeister, M. 2018, A&A, 619, A128 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Zechmeister, M. 2021, MNRAS, 500, 109 [Google Scholar]
 Zeng, L., Jacobsen, S.B., Sasselov, D.D., et al. 2019, Proc. Natl. Acad. Sci. USA, 116, 9723 [NASA ADS] [CrossRef] [Google Scholar]
PLATO Definition Study Report, https://sci.esa.int/s/8rPyPew
All Tables
All Figures
Fig. 1 Orbit geometry of the planetmoon barycenter (M_{b}) 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 Planetmoon 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 Sunlike star, a Jupitersized planet in a oneyear orbit around that star, and a Neptunesized 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 dotteddashed 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 m_{V} ~ 11 Sunlike star from the PLATO mission. Digital star colors are from Harre & Heller (2021). 

In the text 
Fig. 4 Simulated planetmoon 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 m_{v} ~ 11 photometrically quiet Sunlike 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 R_{s} is given in units of R_{⊙}, a_{b} in units of AU, and M_{p} in units of Jupiter masses. Due to the small number of transits, the posterior is multimodal for combinations of P_{pm} and τ_{pm}. Additional transits would constrain the system much better. 

In the text 
Fig. A.1 Transit tangential velocity of the planetmoon barycenter during ingress (v_{t},_{in}), midtransit (v_{t}), and egress (v_{t},_{eg}) as a function of the orientation of the periastron with respect to the line of sight (c_{b}). (a) The planetmoon barycenter has an orbital period of 30 days around a Sunlike star. The upper panel shows v_{t,in} (red solid line), v_{t} (black dashed line), and v_{t,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 Δv_{t} = o_{t},_{in}/o_{t},_{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 intransit 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 worstcase orientation of periastron and a Jupitersized 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 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 Earthsized body in front of a Sunlike 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 smallbody approximation with linear interpolation. Left: For R_{p}/R_{s} = 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 R_{p}/R_{s} = 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 MandelAgol model (black) is the most precise, but also slowest, reference version. The hybrid model (red) uses interpolation between the smallbody and the MandelAgol methods, and it is applicable for R_{P}/R_{S} < 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 smallbody approximation is the fastest method, but only applicable (errors below 1 ppm) for R_{P}/R_{S} < 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 Corei7 1185G. 

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