Modelling the Light Curves of Transiting Exomoons: a Zero-order Photodynamic Agent Added to the Transit and Light Curve Modeller

Despite the ever-growing number of exoplanets discovered and the extensive analyses carried out to find their potential satellites, only two exomoon candidates, Kepler-1625b-i and Kepler-1708 b-i, have been discovered to date. A considerable amount of effort has been invested in the development of algorithms for modelling, searching, and detecting exomoons in exoplanetary light curves. In this work, we incorporate moon handling capabilities into the state-of-the-art and publicly available code, the Transit and Light Curve Modeller (TLCM). The code is designed for the analysis of transiting exoplanet systems with the inclusion of a wavelet-based noise handling algorithm. Here we present an updated version of TLCM that is capable of modelling a coplanar planet-moon system on an elliptical orbit around its host, accounting for mutual eclipses between the two bodies (and neglecting perturbative effects) -- a so-called photodynamic model. The key benefit of this framework is the ability for a joint analysis of multiple planet-moon transits. We demonstrate the necessity of this software on a case study of Kepler-1625b. Similarly to prior works, we conclude that there is no firm evidence of an exomoon in that system, by showing that temporally correlated noise can mimic apparent lunar transits.


INTRODUCTION
Since the discovery of the first extrasolar planet orbiting main sequence stars (e.g.Mayor & Queloz 1995), the field of exoplanet research grew rapidly and became increasingly important in astronomy.Based on our understanding of the Solar System, the quest to find satellites of these planets (called 'exomoons') began to expand shortly thereafter (e.g.Sartoretti & Schneider 1999).Although there are > 5000 known exoplanets to date, only two exomoon-candidates have been identified (Kepler-1625b-i;Teachey & Kipping 2018) and (Kepler-1708 b-i; Kipping et al. 2022).Heller et al. (2019) offer an alternative no-moon explanation for the signal observed by Teachey & Kipping (2018), involving the presence of time-correlated noise and the properties of Bayesian inference.In Kreidberg et al. (2019), the authors also argue that the there is no conclusive evidence for an exomoon in the transit light curves of Kepler-1625b.One can therefore argue that the transiting exomoon candidate proposed by Teachey & ★ E-mail: xilard1@gothard.hu(SzK) Kipping (2018) is the result of temporally variable noise mimicking the transit of a moon.There are two distinct sources of this noise type (also commonly called red noise): astrophysical and instrumental.The former of these include the presence of stellar pulsations, flares, granulation, and spots, which may produce light variations similar to the ones expected from transits of solid bodies (e.g.Silva-Valio et al. 2010).The latter group includes cosmic rays causing excess flux, temperature variations of the detectors/telescope assembly, pointing drifts etc (see e.g.Csizmadia et al. 2023, for more discussion on this topic).Time-correlated noise can manifest on a wide range of timescales.The long-term trends can easily be accounted for by linear or parabolic baselines (Teachey & Kipping 2018).Kreidberg et al. (2019) suggests that instrumental noise sources (of shorter time-scales) may be behind the apparent exomoon transit.Furthermore, Fox & Wiegert (2021) found six exomoon candidates (KOIs-268.01, 303.01, 1888(KOIs-268.01, 303.01, .01, 1925(KOIs-268.01, 303.01, .01, 2728.01, 3320.01) .01, 3320.01)based on Transit Timing Variations (TTVs).Kipping (2020) disproved these specific exomoon candidates from observational aspects, while Quarles et al. (2020) suggest that the systems mentioned above could not be stable from a dynamical perspective.The field of exomoon research is intertwined with current and upcoming spacebased telescopes, such as CHEOPS (Simon et al. 2015;Benz et al. 2021;Ehrenreich et al. 2023), PLATO (Rauer et al. 2014) and Ariel (Tinetti et al. 2021;Szabó et al. 2022), the number of potentially discoverable exomoons is therefore expected to rise.
In recent decades, a number of methods were proposed to detect and characterise exomoons.These include, but are not limited to, the analysis of TTVs (Simon et al. 2007;Kipping 2009), the Rossiter-McLaughlin effect (Simon et al. 2010), direct modelling of the transits (e.g.Kipping 2011), finding radio signatures caused by the interaction of the moon and its host planet (e.g.Noyola et al. 2014Noyola et al. , 2016)), spectrophotometry (Jäger & Szabó 2021) and direct imaging (Kleisioti et al. 2023).As expressed above, some of these methods were used to identify exomoon candidates.The long-term stability of the planet-moon systems is also favoured from an observational point of view.A considerable effort was made to assess the dynamics of the formation and evolution of such systems in general (Domingos et al. 2006;Spalding et al. 2016;Alvarado-Montes et al. 2017;Sucerquia et al. 2019), and the satellites of Kepler-1625b/Kepler-1708b (Moraes & Vieira Neto 2020;Moraes et al. 2023).The topic of habitability is also often discussed in relation to the exomoons (e.g.Kipping et al. 2009;Kaltenegger 2010;Dobos et al. 2022).
Following the success of the transit method in the detection of exoplanets, several algorithms have been put forward to model the light curve of a transiting planet-moon system (e.g.Luger et al. 2017;Kipping 2021;Saha & Sengupta 2022;Gordon & Agol 2022).At the time of writing however, there is only one open-source software available that is capable of fitting the light curves of transiting planet-moon systems out-of-the-box (Hippke & Heller 2022).Recently, time-correlated noise (either of stellar or instrumental origin) has been demonstrated to mimic exomoon light curve components, and is recognised to be a major source of false alarms (e.g.Ehrenreich et al. 2023).Therefore, there is a clear need for algorithms with a combination of noise-handling capabilities and the photodynamic effects.Red noise can arise from instrumental and astrophysical sources, such as stellar spots, microflares, pulsation etc. (see Csizmadia et al. 2023;Kálmán et al. 2023, for more detailed discussions).To filter out the false-positive moon signals that are consistent with noise level fluctuations, a light curve inversion algorithm that accounts for red noise is needed.
In this paper, we introduce the photodynamic update of the Transit and Light Curve Modeller (TLCM; Csizmadia 2020), which uses the well-established wavelet-based noise filtering algorithm of Carter & Winn (2009).Its benefits in the precision and accuracy of exoplanet parameters in transit light curves are well established (Csizmadia et al. 2023;Kálmán et al. 2023).The newly integrated model also includes the ability to model mutual transits by the planet and its satellite.In Section 2, we describe the calculation of the position of the planet and its satellite.In Section 3, we provide example light curves with arbitrary parameters and a specific case of Kepler-1625bi (Teachey & Kipping 2018).In Section 4, we draw our conclusions.

METHODS
The TLCM software is designed to model a wide range of astrophysical effects that can be observed in light curves of transiting exoplanets, including transits, occultations, phase curve variability, ellipsoidal effect and Doppler beaming (Csizmadia 2020;Csizmadia et al. 2023).Special emphasis is given to the handling of timecorrelated noise: the wavelet-based noise filtering of Carter & Winn (2009) is used to remove (or reduce) the amount of red noise present in light curves (Csizmadia et al. 2023;Kálmán et al. 2023).This noise filtering maximises the likelihood of the combination of the noise and transit models, thus 'forcing' the signal into the shape with the highest likelihood.
To add the transits of exomoons to TLCM, we make the following assumptions.The common center of mass (CCM) of the planet-moon system has an elliptical motion around the host star, and both of them are on circular orbits around the CCM.Furthermore, the planet and its satellite orbit the CCM in the same plane as the CCM orbits the star.
These assumptions help us avoid high degrees of degeneracy (Simon et al. 2010;Hippke & Heller 2022), and to limit the number of free parameters to five to describe the moon (Mandel & Agol 2002;Csizmadia 2020).These are the scaled semi-major axis of the moon,  Moon / Planet , the moon's relative radius  Moon / Planet , its relative mass,  =  Moon / Planet , its orbital period  Moon , and Δ -a phase shift between the epoch of the CCM and the transit of the satellite.The complete transit light curves are computed via numerical integration of the two bodies, whose sky-projected positions are calculated via the equations described below.The model presented here is also capable of accounting for mutual eclipses of the planet and its satellite, as shown on Fig. 1.
We describe these orbits by introducing two phase angles: where  0 is the epoch at which the CCM is collinear with the center of the stellar disk and the observer,  CMC and  Moon are the orbital periods of the CCM around the star and of the satellite around the planet, respectively.The mean anomaly at  0 , denoted by  0 , is calculated via Eqs.(12-15) of Csizmadia (2020).The mean anomaly can then be written as  =  CMC +  0 .Let  CMC and  CMC denote the eccentricity and the argument of the periastron (measured from the tangential plane of the sky) of the CCM, respectively.Using Kepler's equation, the eccentric anomaly, solved iteratively, can be written as  CMC =  +  CMC sin  CMC .The true anomaly, , can then be calculated as In the orbital plane of the planet, the position of the common center of mass at any point in time is where  CMC is the semi-major axis of the CCM, and  ★ is the stellar radius.With the introduction of the moon-to-planet mass ratio,  =  Moon / Planet , and the semi-major axis of the orbit of the moon,  Moon , the position of the planet in the orbital plane can be written as  The position of the moon in the orbital planet of the planet is given by In the equations above, Δ is the phase shift of the epoch of the satellite, and  Planet is the planetary radius.The sky-projected position of the exomoon and the exoplanet are then calculated from Eqs. (7 -12) via the equations described in Sect.2.2. of Csizmadia (2020).
From the known sky-projected positions, the limb-darkened transits of the two bodies are calculated via the commonly used numerical integration tools to account for the light blocked by the planet and the moon.The resulting light curve therefore consists a 'planet component' and a 'moon component' superimposed on each other (Simon et al. 2007).Owing to the assumed coplanarity, the projected orbital trajectories of the planet and the moon will only coincide if the inclination of the orbital plane of the CCM (with respect to the line of sight) is 90 • .This is highligted on Fig. 2.
Furthermore, there is a check for the occurrence of mutual transits (i.e. the moon passing in front of, or behind the planet).If such an event is detected, a correction is applied so that the mutually eclipsed area is calculated only once.A mutual transit event can occur if and only if all of the following conditions are fulfilled: Here  denotes the sky-projected mutual distance of the different objects.When all three conditions are simultaneously fulfilled, the arcs of the circular boundaries between the star, planet and moon are determined to calculate the mutually eclipsed area.This is carried out in a numerical way.A 101x101 pixel image is centered around the moon and planet and we check which points are inside the skyprojected images of the star, planet and moon.The mutually eclipsed area is enclosed by either a full circle (in the limit if the moon is completely behind or in front of the planet, or, conversely if the planet is smaller than the moon), or by two or three arcs (see Fig. 2).The user can choose between 6 × 6, 10 × 10, 20 × 20, 48 × 48 or 96 × 96 integration points which are distributed over the mutually eclipsed area according to the Gauss-Legendre quadrature rules.According to our experiences, a 6×6 grid is enough to get a few ppm precision.The intensity behind the mutually eclipsed area is then calculated with the Gauss-Legendre quadrature and the limb darkening is taken into account.This correction flux is added accordingly to the light curve, because the simple transit events of the moon and the planet would result in a duplicated subtraction of the flux behind the mutually eclipsed area.An illustration of how this works is given by Figs. 1  and 2.
We note that the current implementation of the photodynamic model is highly simplified, as the orbital plane of the moon is considered to coincide with the one of the CCM.Furthermore, the two bodies are also considered to be on Keplerian orbits (hence the term 'zero order').Future updates will focus on the inclusion of the possibility of misaligned orbital plane for the satellite, as well as the inclusion of various perturbative effects, which may lead to variations in periodicity, transit depth and duration etc. (see e.g.Borkovits et al. 2022, for a more detalied discussion on these effects in triple star systems).

A schematic configuration
An example light curve, computed via the methods described in Sect.2, is shown on Figure 1.This system was set up with the following exaggerated, arbitrary parameters:  CMC / ★ = 50,  Moon / Planet = 10.2,Δ = 0.4,  Moon / Planet = 0.53,  Moon = 2.5789 days,  Planet / ★ = 0.51,  = 0.3,  CMC = 0.35,  CMC = 124 • , and  CMC = 21.5089days.For the purposes of this demonstration, we assumed uniform surface brightness of the stellar disk (i.e.no limb darkening) and we enlarged the mass and radius of the moon.Furthermore, we added a synthetic noise model with autocorrelated effects, described by an Autoregressive Integrated Moving Average (ARIMA) model (e.g.Wilson 2016).An   ( , , ) process consists of an autogregressive model of the order  and a moving average model of the order , with a -order integration.In time series analysis, integration can be thought of as the inverse of differencing, a smoothing operation that generates from the time series () the series  ′ () by representing the change between each consecutive element of (), thus: The autoregressive and moving average processes can be represented as Here  1,2 are constants,   and   are the coefficients of the two processes, while the noise terms  are drawn from a Gaussian distribution with zero mean (i.e.white noise).In our case (Figure 1), we set up the model with the parameters   ([1, −0.5], [1], [0.5, −1])1 .Figure 1 clearly shows that the algorithm is capable of handling different angular velocities for the planet and its satellite (manifesting in different transit durations for the two objects).

A realistic configuration
We also explore a realistic planet-moon case.To that end, we adopted the parameters of the planet and its candidate satellite from Teachey & Kipping (2018), with arbitrary limb darkening parameters, while restricting the orbital plane of the moon to coincide with that of the planet.The left and middle columns in Fig. 3 show two transits with arbitrarily positioned moons.In the rightmost column, we demonstrate a false-positive detection in a no-moon case.
The top panels show the synthetic input transits with a white noise component only (with a 200 ppm point-to-point scatter, compared to the ∼ 270 ppm transit depth of the moon).The three columns present two systems with moons (A and B) and a moonless case (C).The middle row shows the same input signals, but with a red noise component added as well.To be realistic, here we introduced the red noise via an   ( [0.6819], [1], [−0.3603, −0.5747]) process, cloning the noise of the residuals of the original measurement of Teachey & Kipping (2018).The bottom row presents a comparison between the "moon" solution and the synthetic planet+moon+red noise.We observe that in cases A and B the presence of the moon is convincing in the synthetic measurements even by visual inspection.On the other hand, in scenario (C) where the moon is absent, we encounter a false positive detection due to the similarity of the rednoise-induced signal to that of the moon's.
One possible way to exclude false detections like this is to schedule observations that are long enough so that a noise filtering algorithm (e.g. the wavelet-filter of Carter & Winn 2009) is able to discern the actual noise properties in detail.
The true advantage of this filtering technique lies in the fact that it does not require prior knowledge of the exact source or realisation of the time-correlated noise, since it operates by constructing a time series model in a wavelet basis.It is therefore efficient in removing the red noise down to the time scales of minutes (Csizmadia et al. 2023).It is also important to stress that the test carried out with the TLCM thus far (although they were limited to transiting exoplanets) indicate that the radius of the transiting object can be recovered with high precision and accuracy even at low S/N (Csizmadia et al. 2023;Kálmán et al. 2023).This would also be one of the most important parameters for a potential exomoon candidate.
Fig. 3 has two important lessons for us: • When chasing a moon, the measurement of the null signal (noise only) before and after the transit has to be as long as or longer than the part of the presumed signal itself, to enable the algorithms to distinguish between noise and signal.
• The threshold of a secure detection depends on the length of the noise sampling as well as the actual geometrical configuration.This has to be tested even at the observation planning phase, and has to be taken into account in reasoning the feasibility of any kind of measurement in this field.
We note however that the possibility of dense rings around Kepler-1625b, or even the as yet unconfirmed Kepler-1625b-i (Sucerquia et al. 2022) can not be excluded, although exploring that possibility is beyond the scope of this paper and will be the subject of a subsequent one.In cases A (left column) and B (middle column), the transits of the moon are present in the synthetic photometry (points), case C (right column) presents the transit of an exoplanet only but leads to false detection because of the correlated noise in combination with short observations.The top row represents the ideal case where only white noise is present on the light curve.The middle and bottom rows also have an ARIMA noise model superimposed upon them.The synthetic photometry in the bottom row is the same as in the middle row, but the solid red line represents a planet-and-moon scenario.We emphasise that these figures are not results of light curve fitting.

DISCUSSION AND SUMMARY
Motivated primarily by the lack of confirmed exomoons and the work of Heller et al. (2019), we developed a complex photodynamical model, capable of describing a planet-moon system on an elliptical orbit around a host star in the presence of red noise.We integrated this model into the TLCM code (Csizmadia 2020).One of the primary features of the TLCM algorithm is the ability to handle time-correlated noise via the wavelet-based routines of Carter & Winn (2009).
It is therefore expected that by combining this noise handling method and the photodynamical model presented here, we will be able to () reject false positive exomoon candidates (Heller et al. 2019), () uncover exomoons 'hidden' by red noise, and () characterise real exomoons with sufficiently high precision and accuracy, as it is commonly done for exoplanets (Csizmadia et al. 2023).
The results shown in the paper can be summarised as follows: • The Exomoon Agent embedded in the TLCM algorithm is introduced.With this development, users can now search for planet+moon solutions in addition to planet-only cases, and will be able to chose the best model family based on the e.g.Akaike and Bayesian Information Criteria.
• As a case study, we investigated possible transits of a planetmoon system with parameters taken from Teachey & Kipping (2018).We added correlated noise to the simulations using ARIMA clones of the actual residual of the observations.
• We observed that correlated noise components can mimic the characteristic pattern of exomoons, leading to false positive detections.Our results are in accordance with Heller et al. (2019).
• The single way to avoid these false detections is to take long enough observations.They must include the pre-and post-transit areas of null signal to teach the algorithm the noise properties in details, that can lead to a reliable mitigation.The threshold of a suspected detection depends on the sampling of the observed data, which must be calibrated by accounting for the length of the 'red noise only' (or out-of-transit) light curves.
• With the currently available instruments (including those that are in advanced stages of development such as Ariel), we will likely never be able to claim with full confidence that an exomoon has been detected by modelling a single transit.This is partly because the proper separation of time-correlated noise of astrophysical and instrumental origin itself can be challenging in and of itself, even when the instrumental effects are well understood as in the case of CHEOPS (Hoyer et al. 2020).In order to separate the temporally correlated noise from the signal of a transit, the desirable strategy of observing exomoons with HST, JWST or even Ariel would have to rely on multiple transits, spanning years (due to the relatively long orbital periods of the exomoon-hosting planets).This would, ideally, allow observations of the target system in multiple configurations (Fig. 2), thus providing a deeper understanding of the probability of false alarm.A strong claim to the presence of an exomoon would have to be the change in occurrence of its transit from observation to observation (relative to the transit of its host).Ideally, TTVs would also have to be observed, and a thorough follow-up would be needed.
• While an understanding of systematic noise sources is beneficial, the wavelet-based noise handling of TLCM can reasonably be expected to allow us do distinguish between an actual moon and a false detection.This has been the inspiration behind the Exomoon Agent upgrade of TLCM.
The software is available publicly for further tests and uses.The thorough validation of the parameter-retrieval capabilities is beyond the scope of this Letter, and it remains in the user's responsibility to look for the stability of solutions in the actual modelling.In an upcoming publication, we will test the parameter stabilities in simulated data sets and data from real measurements.It is expected however, that the Exomoon Agent can be used to reanalyse existing data.
We expect that the current chase for exomoons will remain with throughout the operations of the current space telescopes.In the search for ever smaller satellites, orbiting planets whose hosts are ever farther away, there will be a clear need for data from larger and more sensitive future space-born telescopes, such as the Large UV Optical Infrared Surveyor (LUVOIR; The LUVOIR Team 2019).

Figure 1 .
Figure1.Simulated light curve of an arbitrary transiting planet-moon system (blue dots and solid red curve).Left.The moon remains hidden behind the planet during the transit.Middle.A mutual transit event occurs, requiring flux correction (solid blue curve).Right.The transit of the satellite occurs after the transit of the exoplanet.Note that, for demonstrative purposes, the sizes of both bodies are enlarged, and we assume no limb darkening of the host star.

Figure 2 .
Figure 2. Possible mutual eclipse cases with an illustration of the Gauss-Legendre quadrature integration points (blue) distributed during the mutual eclipses.The three discs represent, in decreasing order of diameter, the star, the planet and the moon.The projected orbits of the planet and the moon are shown with dashed black and magenta lines.The dashed blue line is used to visualise the impact parameters of the two bodies.The limb darkening of the stellar disc is generated via solar-like parameters according to the routines ofHarre & Heller (2021).The mutually eclipsed area (by the planet and the moon) is outlined with solid red lines.Case A: the planet fully occults the moon (or the moon totally transits the planet, i.e. the moon can be behind or in front of the planet).Case B:The planet and the moon exhibit partial mutual eclipse but both totally eclipses the star.Case C: During the partial eclipses of the planet and/or the moon, three arcs may be needed to find the mutually eclipsed area.Case D: no mutual eclipse.

Figure 3 .
Figure 3. Demonstration of time-correlated noise mimicing transits of an exomoon.In cases A (left column) and B (middle column), the transits of the moon are present in the synthetic photometry (points), case C (right column) presents the transit of an exoplanet only but leads to false detection because of the correlated noise in combination with short observations.The top row represents the ideal case where only white noise is present on the light curve.The middle and bottom rows also have an ARIMA noise model superimposed upon them.The synthetic photometry in the bottom row is the same as in the middle row, but the solid red line represents a planet-and-moon scenario.We emphasise that these figures are not results of light curve fitting. )