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/0004-6361/202243129 | |
Published online | 14 June 2022 |
Pandora: A fast open-source exomoon transit detection algorithm
1
Sonneberg Observatory,
Sternwartestraße 32,
96515
Sonneberg,
Germany
e-mail: michael@hippke.org
2
Visiting Scholar, Breakthrough Listen Group, Berkeley SETI Research Center, Astronomy Department,
UC Berkeley
USA
3
Max-Planck-Institut für Sonnensystemforschung,
Justus-von-Liebig-Weg 3,
37077
Göttingen,
Germany
e-mail: heller@mps.mpg.de
4
Institut für Astrophysik, Georg-August-Universität Göttingen,
Friedrich-Hund-Platz 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 star-planet-moon system is computed with a high accuracy as a nested Keplerian problem. We have optimized Pandora for computational speed to make it suitable for large-scale exomoon searches in the new era of space-based high-accuracy 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 Neptune-sized exomoon in a one-year orbit around a Sun-like 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, mV = 11 Sun-like star for practicality. We recovered the simulated system parameters with the UltraNest Bayesian inference package. The run-time of this search is about five hours on a standard computer. Pandora is the first photodynamical open-source 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 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 algorithm1 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 Mp) and its moon (with mass Mm) orbit their common center of mass (with mass Mb = Mp + Mm) 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 (t0), the orbital period of the barycenter around the star (Pb), the semi-major axis of the barycenter orbit around the star (ab), and the transit impact parameter of the barycenter (bb). In general, the transit impact parameter can be parameterized in terms of the orbital inclination with respect to the line of sight as per bb = ab tan(π/2 – ib)/Rs, where ib = π 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 Rs, that is, as (Rp/Rs), (Rm/Rs), and (ab/Rs). By default, the orbit of Mb is assumed to have zero eccentricity (eb = 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 ipm, longitude of the ascending node (Ωpm), semimajor axis (apm), orbital period (Ppm), and time of periap-sis passage (τpm). In Pandora, τpm is normalized with respect to Ppm, 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 τPpm. For eccentric orbits, the number of orbital elements increases to six, including epm and the argument of periapsis (ωpm). Finally, the moon mass Mm is required to model the motion of the planet and moon around their joint barycenter.
The star is parameterized by Rs and two limb darkening coefficients (u1u2) 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 (Ms) follows directly from the barycentric orbital mean motion nb = 2π/Pb and Mb via Kepler’s third law as 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 Mb 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 sky-projected velocity. We describe this direction by a unit vector et and we introduced a unit vector er for the (negative) radial velocity component:
(2)
(3)
The projection of u on er, which we refer to as υr, is given as
(4)
where a = <(υ, er) (Fig. 1). The cosine of this angle between υ and er can be calculated as
(5)
Plugging Eq. (5) into Eq. (4) yields
(6)
Analogously, projection of u on et, 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 point-like Mb 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 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 epm = 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 .
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, 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.
Model parameters.
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 Rp/Rs < 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 < Rp/Rs < 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 < (Rp/Rs) < 0.2 and distances from the stellar center 0 < z < 1 + (Rp/Rs) 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 Rp/s > 0.05, 48million for Rp/Rs = 0.02, and 110million for Rp/Rs < 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
(12)
while the distance between their centers (dpm) 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 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 tracker2. 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 u1 = 0.4089 and u2 = 0.2556 as per Claret & Bloemen (2011) for a star with a surface gravity of log(g) = 4.5, and an effective temperature of Teff = 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 (Mm/Mp = 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 Rp (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 Rp around Jupiter. The planet-moon orbit is slightly inclined with respect to the line of sight (ipm = 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 ≤ mV ≤ 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 sample3, but it ignores any kind of astrophysical variability. Specifically, the photon noise from this test star would correspond to an apparent visual magnitude mV ~ 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 . 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 (), 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 and
based on q1 and q2 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 Ppm/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 Mp + Mm. 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 (cb - y), with y = arctan(Rs/r). During egress, the true anomaly is (cb + 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 |vt| from Eq. (10) with |utin| 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 cb. For all curves, the stellar radius was set to Rs = Rø and the stellar mass to Ms = Mø.
In Fig. A.1(a) we consider a planet-moon barycenter with an orbital period of 30 days, which means that ab = 0.19 AU. The three curves assume eb = 0.01, eb = 0.1, and eb = 0.2, respectively (see labels). The black dashed lines refer to |vt|, the orange solid line to |«yn|, and the light blue solid line to |ot,eg|. The lower panel shows the resulting Δvt 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 eb = 0.2, roughly 2.5 × 10-3 for eb = 0.1, and approximately 2.5 × 10-4 for eb = 0.01. These maximum errors occur for cb = 90° and cb = 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 eb = 0.01, eb = 0.1, and eb = 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 (Rp = 0.1 Rs) in the worst-case scenario for these above cases, in which Pb = 30 d, eb = 0.2, and cb = 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 + eb cos(cb)) 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 Rp. 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 |
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 2Rp/(πapm). 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 Rp/Rs < 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 Rp/Rs = 0.05, the approximation can be used for 65% of the data points in a light curve, growing to 95% for Rp/Rs = 0.02, keeping small < 1 ppm errors Fig. C.1. For bodies Rp/Rs > 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 Fx(ɀ, k, u1, u2) for only two radial distances ɀ = 0 and ɀ = 0.65, the latter being determined empirically (k is the occulter radius in units of Rp/Rs, and u1 and u2 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 F1 and F2 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
(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 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
- Agol, E., Luger, R., & Foreman-Mackey, 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 astronomy-related 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, DSTO-TN-0722, 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., Lustig-Yaeger, J., & Agol, E. 2017, ApJ, 851, 94 [NASA ADS] [CrossRef] [Google Scholar]
- Luger, R., Agol, E., Foreman-Mackey, 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 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 |
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 |
Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.