The Advantages of Global Photometric Models in Fitting Transit Variations

Estimation of planetary orbital and physical parameters from light-curve data relies heavily on the accurate interpretation of Transit Timing Variation (TTV) measurements. In this letter, we review the process of TTV measurement and compare two fitting paradigms—one that relies on making transit-by-transit timing estimates and then fitting a TTV model to the observed timings, and one that relies on fitting a global flux model to the entire light-curve data set simultaneously. The latter method is achieved either by solving for the underlying planetary motion (often referred to as “photodynamics”), or by using an approximate or empirical shape of the TTV signal. We show that, across a large range of the transit S/N regime, the probability distribution function of the mid-transit time significantly deviates from a Gaussian, even if the flux errors do distribute normally. Treating the timing uncertainties as if they are distributed normally leads, in such a case, to a wrong interpretation of the TTV measurements. We illustrate these points using numerical experiments and conclude that a fitting process that relies on a global flux fitting, rather than the derived TTVs, should be preferred.

Interpreting photometric light curves is increasingly common in the age of multiple ground-based and space observatories such as Kepler (Borucki et al. 2010) and TESS (Ricker et al. 2010), which led to the discovery and characterization of thousands of exoplanets.The transit method provides information regarding the orbital period, phase, and several normalized geometrical properties: the planetary size, semi-major axis, and impact parameter -all relative to the stellar radius (e.g.Seager & Mallén-Ornelas 2003).If the stellar mass is known, the absolute semi-major axis can be inferred, and if the stellar radius is known, the absolute planetary radius can be inferred.The planet's mass, though, cannot be obtained using periodic transits alone.The eccentricity of the orbit is related to the transit duration (e.g.Carter et al. 2008), but is degenerate with the orbital orientation.To infer the planetary mass and orbital shape, perturbations among the planets need to exist in the system.In this case, Transit Timing Variations (TTV, Agol et al. 2005;Holman & Murray 2005) arise; this phenomenon was found abundantly in the Kepler population (Holczer et al. 2016;Ofir et al. 2018).
Multiple authors have extensively studied the functional form of the TTV signal.Of specific importance was the introduction of the fundamental super-period TTV (Lithwick et al. 2012), which enabled a fast linear inversion of the sine-like TTV signal to masses and eccentricities, up to the degeneracies involved.These degeneracies can be removed if a synodic "chopping" signal is visible (Deck & Agol 2015) or if second-order TTV modes are visible (Hadden & Lithwick 2016).Another avenue for using TTVs to infer planetary properties is to use full N-body integrations to predict transit times (e.g.Deck et al. 2014), driven by a non-linear fitting process that compares the predicted times with the time estimates obtained from the data.Such methods were used to interpret the planetary parameters in many systems (Hadden & Lithwick 2014, 2017;Jontof-Hutter et al. 2014, 2016, andmany more).
Light-curve data that contains TTVs can be interpreted in two different methods.One method is to perform a global fit on the light-curve photometric data, either by using a dynamical model of the planetary motion (often called photodynamics, e.g.Carter et al. 2011;Grimm et al. 2018;Judkovsky et al. 2022a), or by using an empirical TTV functional form (e.g.Ofir et al. 2018).The second method is to estimate the individual properties of each transit, typically the mid-transit times, thus transforming the light-curve flux data into timing data, and then to interpret these to estimate the planetary properties via the measured pattern of the TTVs (e.g.Hadden & Lithwick 2016;Jontof-Hutter et al. 2016;Hadden & Lithwick 2017).
In this work, we compare the performance of these two methods at various single-transit SNR values.We do not present a full statistical analysis of the errors, but instead, we aim to illustrate the consequences of utilizing typical practices in light-curve and TTV fitting, highlighting the applicable methods and their advantages or shortcomings.We achieve this goal by performing a numerical experiment to show that the inaccurate propagation of the flux errors into transit timing errors can yield wrong results in the TTV analysis.This elaborates on and clarifies the arguments presented in Ofir et al. (2018, §2.3).In addition, we stress that the following analysis does not aim to show that full photodynamical modeling outperforms a transit-by-transit timing analysis in a particular exoplanetary system or a particular modeling code but rather to show that performing a fit on the photometric data outperforms a transit-by-transit analysis, regardless of the underlying model from which the TTV shape is derived.

TRANSIT TIMING FROM FLUX
A prerequisite for the method that relies on using the individual times of mid-transit is the translation of the light-curve flux data to a set of individual transit time estimates, usually done using a non-linear fitting algorithm.The timing values and errors are then used to fit the TTV model parameters, which are functions of the planetary physical parameters.An underlying assumption within this process is that a Gaussian can approximate the posterior probability density function (PDF) (which, given a uniform prior, is equal to the likelihood of the data) as a function of the transit time.In other words, the assumption is that the propagation of the flux errors to timing errors can be done by summarizing the PDF properties in two numbers: the expectation value and standard deviation.This method was used, for example, in several ground-breaking works dealing with exoplanetary parameter estimation, such as Hadden &Lithwick (2016, 2017, andothers).See formulation in Hadden & Lithwick (2016, Appendix B.1).As we show below, this assumption might not hold even for perfectly normally-distributed flux values for a large range of transit SNR values, where the individual transit SNR is defined as where δ is the transit depth, σ F is the typical uncertainty on the flux, N pts is the number of data points within a single transit event, and R p /R * is the planet-to-star radius ratio.The transit SNR is related to the measurement precision of individual mid-transit time.A deeper transit would have more distinctive ingress and egress, enabling a more accurate timing measurement.The transit timing error is a complicated function of the planetary parameters, and in practice, it is obtained from a non-linear fitting process.We are not aware of a closed-form formula for this error.However, Holczer et al. (2016) provides a simplified empirical benchmark, which states that the timing error is roughly 100 min/SNR (with some spread around this value).
As we show below, the PDF of the mid-transit timing is not Gaussian.This deviation grows as the SNR decreases.Adequately propagating the flux errors to timing errors requires not only the mean and standard deviation of the PDF but the entire PDF shape.Propagating the expectation value and standard deviation alone and using them in the TTV fitting as if they represent a Gaussian can yield wrong results, especially when modeling shallow transits.
Assuming that the noise on the flux distributes normally, a global flux fit that integrates the information of all transits simultaneously would not suffer from this problem.In the language of the SNR definition given above, when fitting the entire light curve simultaneously, the effective SNR increases with the square root of the number of transits.
Real photometric data may not distribute normally and suffer from red noise, further complicating the light-curve analysis in both methods.Here, we neglect this additional complication and assume that the flux noise distributes normally.We note that if the full PDF of each transit was measured and propagated to the TTV model then, in principle, the transit-by-transit method would have access to the full information needed to produce an equivalent fit to that of the global flux model (up to the appearance of transit variations other than TTV).The issue that makes the transit-bytransit technique problematic is that the PDF of the individual transit time does not distribute normally, even if the flux data does.We demonstrate this argument in the following numerical experiment.

NUMERICAL EXPERIMENT
In this section, we mimic the process of interpreting a TTV-containing light curve using two methods: In the first method, we estimate the individual transit times and then use them to fit a TTV model; in the second method, we use a global model to fit the light-curve flux directly.Importantly, to isolate these factors we assume perfect knowledge of all other elements, such as the properties of the plant, star, or TTV model shape.

Data Generation
The simulated flux data is generated as follows.First, we calculate a sequence of nominal times of mid-transit for a planet with an orbital period P = 10.13 days and with a phase that sets the first nominal transit time at T mid,0 = 1 day, along the time span t = 0..1500 days.We obtain a series of 148 nominal times of mid-transit.Then, we add to these nominal transit times a sinusoidal TTV of the form TTV = A cos ωt + B sin ωt, where A = 0.01 day and B = 0.02 day are constants that together set the amplitude and phase of the TTV, ω = 0.015 Rad/day is the TTV's angular frequency, and t, in this case, is the nominal transit time for which the TTV is calculated.Adding the TTV to each nominal transit time yields a set of actual transit times.We then assume a circular orbit with a semi-major axis-to-star ratio a/R * = 15 and an impact parameter b = 0.3 and use the obtained orbital geometry to construct a full light-curve using the Mandel-Agol formula (Mandel & Agol 2002) with limb-darkening parameters u 1 = 0.36 and u 2 = 0.28 (Claret 2000).All the constants were chosen at values commonly found in Kepler data.This process results in a light curve without any measurement errors.For simplicity, we neglect the effect of binning (Kipping 2010), although binning must be addressed for real data.The next step in the flux data generation process is to add noise.We assume a white noise model in which the noise on each data point is distributed as N (0, σ 2 F ), that is normally with a zero mean and σ F standard deviation.
We repeat the experiment several times, varying σ F over a wide range.In reality, the noise spectrum is not white, and a red component is present; we elaborate on this later.

Data Interpretation
Having obtained noised light-curve data, we interpret it using the abovementioned methods.The interpretation process aims to find the best-fitting TTV frequency ω and coefficients A, B. In the first method, in which the transit times are estimated one by one, we first treat each transit individually.We keep only data points within 0.5 days from the nominal time of mid-transit.We assume perfect knowledge of the planet-to-star radius ratio, the semi-major axis, and the impact parameter.We optimize only for the transit time using the Levenberg-Marquardt (LM) algorithm.This algorithm provides a best-fitting value and an estimated variance on the (only) fitted variablethe mid-transit time.The starting value for the fitted variable is the nominal time of the individual transit under consideration.This starting point should be sufficiently close to the true value to enable convergence because our TTV amplitude is significantly smaller than the transit duration.Having obtained a set of mid-transit times, we subtract from them the linear fit to these times and obtain the simulated measured TTV signal.We then perform a second non-linear fit using the LM algorithm, where this time, the fitted variable is the TTV frequency ω.Because for each ω value, the best fitting A, B can be found using linear regression, we can optimize the single variable ω to obtain a best-fitting TTV model.In this case, the starting point for fitting ω is set to the true underlying ω value.This would be the case, for example, when the perturbing planet is known because then the TTV super-period can be calculated (Lithwick et al. 2012).This process results in a best-fitting model of the parameters ω, A, B. We note that for our finite data span, a simultaneous fit for a line and a sinusoid would, in principle, yield a better match than a sequential fit of a line and then a sinusoid; however, given that the super-period is shorter than the data span, the induced error is much smaller than the TTV amplitude.
In the second method, we perform a non-linear fitting for the flux data, simultaneously fitting ω, A, B using the non-linear LM algorithm over three variables.This process results in estimates for the values of A, B, ω, and their errors.

Results -an Illustrative Example
We highlight the results demonstrating the difficulty of propagating flux estimates to time estimates -a limitation arising from the non-Gaussian distribution of the time-of-mid-transit for an individual event.In Figure 1, we show examples of the fit quality of three individual events, each taken from a different numerical experiment with a different single-transit SNR value.This figure demonstrates that for an SNR value significantly larger than unity, there is a clear global minimum in the χ 2 surface which is similar, but not equal to, a parabola, which means that the likelihood distribution is close to, but not exactly, a Gaussian.As the SNR approaches unity from above, the shape of this minimum deviates more and more from a parabola.When the SNR value is smaller than unity, other local minima can appear, overprinting the true transit signal by noise.The black dashed line denotes the true time of mid-transit.Bottom row: For each event, we show the χ 2 surface over different values of time of mid-transit (blue) and a best-fitting parabola for the points within the 1 σ confidence interval.The small inserts show a zoom-in on the region where the parabolic fitting was performed, with a circle indicating the parabola vertex.In the high SNR case (left-hand side), there is a clear depression in the χ 2 surface, with a shape similar, but not equal to, a parabola around the best fitting value.As SNR is reduced towards unity (middle), there is still a minimum near the true time of mid-transit, but its shape deviates from a parabola even more than in the high-SNR case.In SNR values lower than unity, the transit shape is erased by the noise (note the relatively shallow vertical span of the entire panel), and other local minima begin to appear in the χ 2 surface, thus leading to a completely wrong estimate of the time of mid-transit.
In Figure 2, we show the fit for the TTV frequency for a high-SNR and a low-SNR case, assuming a normal probability distribution of each individual mid-transit time.The resulting frequency value corresponds to the true TTV frequency in the high-SNR case but significantly differs from the true frequency in the low-SNR case.It is also shown that without prior knowledge of the starting point for ω, a random frequency would be preferred based on a Lomb-Scargle periodogram.This figure demonstrates the end effect of treating the PDF of each individual event timing as a Gaussian: the TTV fit that relies on the individual time estimates and their associated errors may result in a wrong TTV frequency.
In Figure 3, we show the ∆χ 2 surface as a function of the TTV frequency ω for both methods: fitting individual mid-transit times and fitting a global model over the flux.Along this surface, the parameters A, B are kept constant at their true values.Both methods yield a clear global minimum around the true TTV frequency in the high-SNR case.The transit-by-transit method generates a wider minimum and generally a different likelihood surface, explained by the incorrect error propagation of the individual timing estimates.In the low-SNR case, on the other hand, significant local minima appear due to the noise.These minima exist in both the timing and the flux methods, but they are deeper in the timing method, while in the flux method, the global minimum is still near the true frequency.We attribute the selection of the wrong frequency in the timing method to the incomplete propagation of errors, a process that assumed a normal distribution of each transit mid-time, whose incorrectness becomes larger and larger as the SNR decreases.Furthermore, we note that a wrong ω value would lead to wrong A, B values that would bias the χ 2 even more, thus making the erroneous global minimum even deeper.Fitting a global flux model does not suffer from Figure 2. The best-fitting TTV solutions for two numerical experiments, each with a different singletransit SNR value.The blue error bars are the TTV estimates from fitting individual transit times, and the solid orange curve is the best-fitting TTV model, obtained from an LM search beginning at the correct ω value.The yellow curve on the right-hand-side plot is the best-fitting TTV model in the Lomb-Scargle peak frequency.At the high-SNR case (left-hand side), the true TTV frequency is apparent in the estimated TTV obtained from a transit-by-transit analysis, leading to a correct estimation of the TTV frequency ω up to the estimated error.In the low-SNR case (right-hand side), the noise has erased the transit events, thus leading to wrong or highly uncertain estimates of the individual mid-transit times.The estimated TTV frequency, in this case, is far from the true value.
this problem: Because, in our experiment, the flux values are distributed normally, and the global minimum in the χ 2 surface is also kept for low-SNR values since the integration of all data points is done appropriately, in the χ 2 definition sense.

Results -a Population Analysis
As explained above, in the transit-by-transit method, the inaccurate error propagation results in an inaccurate estimation of the TTV parameters.We illustrate this in another way by running multiple experiments as described above and checking the distribution of the results among them.For each experiment, we calculate the deviation of the best-fitting parameters from the true value, denoted by ∆ω, ∆A, ∆B, and the error estimation in each parameter, denoted by σ ω , σ A , σ B .If the estimates distribute normally, then we expect that along multiple experiments with different random noise values, the quantities ∆ω/σ ω , ∆A/σ A , ∆B/σ B will distribute as N (0, 1).
In Figure 4 we illustrate the non-Gaussian distribution of the parameters which is attributed to the non-Gaussian PDF of the mid-transit time for each individual event.In this figure, we show the distributions of the quantities ∆ω/σ ω , ∆A/σ A , ∆B/σ B among 15,505 experiments for both the transit-by-transit the global fit method driven by the LM algorithm.We show this for two SNR values.In the high-SNR case, the global fit generates a distribution of the parameters strikingly close to a Gaussian, while the transit-by-transit method generates a broader distribution, clearly indicating a non-Gaussian distribution of the estimated parameters.In the low-SNR case, the distribution deviates from a Gaussian for both methods, albeit to a lesser degree for the global fit method.Interestingly, the transit-by-transit method generates an underestimate of the TTV amplitude.We tested and empirically found that this phenomenon occurs consistently for various Figure 3.The χ 2 surface in the timing method and the flux method for the high-SNR case (left-hand side) and for the low-SNR case (right-hand side).The ∆χ 2 surface (relative to the minimum) in the method using individual mid-transit times fitting is shown in blue; the ∆χ 2 surface for the global flux method is shown in orange.A dashed black line indicates the true TTV frequency.In the high-SNR case, both methods generate a clear global deep minimum around the true TTV frequency.In the low-SNR case, additional, deeper minima are found around other frequencies in the timing method.Local minima exist in the flux fitting method, but the global minimum remains around the true TTV frequency.
values of A, B. That is, the total TTV amplitude √ A 2 + B 2 is underestimated, even of high transit SNR.The systematic underestimate of the TTV amplitude would directly translate to an underestimate of planetary mass.This finding could contribute the previously observed underestimate of TTV masses with respect to RV discussed in Steffen (2016) and Mills & Mazeh (2017).modeling bias also result in an increase in the mass estimates as more data is accumulated, instance, in the mass estimates of the TRAPPIST-1 system between Grimm et al. (2018) and Agol et al. (2021).
In the appendix, we show plots showing the spread of estimated values along the sample experiments.Top row: the high-SNR case.In this case, the global fit method generates parameter estimates that nicely fit a Gaussian, while the transit-bytransit method generates a broader The black dashed indicates a normal distribution with unity standard deviation.Bottom row: the low-SNR case.In this case, both distributions deviate from a Gaussian; the transit-by-transit method generates a much broader distribution, for which the peak is almost erased, while the global fit remains well-centered.

Figure 1 .
Figure 1.Three individual events of different SNR values.Top row: The blue error bars are the synthetic data points, the solid orange line is the best-fitting model, the vertical black dashed line is the true time of mid-transit, and the pink transparent stripes are the 1,2 and 3 σ confidence intervals of the time of midtransit, as retrieved by the LM algorithm (on the left panel they are too narrow to be distinguishable).The black dashed line denotes the true time of mid-transit.Bottom row: For each event, we show the χ 2 surface over different values of time of mid-transit (blue) and a best-fitting parabola for the points within the 1 σ confidence interval.The small inserts show a zoom-in on the region where the parabolic fitting was performed, with a circle indicating the parabola vertex.In the high SNR case (left-hand side), there is a clear depression in the χ 2 surface, with a shape similar, but not equal to, a parabola around the best fitting value.As SNR is reduced towards unity (middle), there is still a minimum near the true time of mid-transit, but its shape deviates from a parabola even more than in the high-SNR case.In SNR values lower than unity, the transit shape is erased by the noise (note the relatively shallow vertical span of the entire panel), and other local minima begin to appear in the χ 2 surface, thus leading to a completely wrong estimate of the time of mid-transit.

Figure 4 .
Figure4.The distribution of the deviation of the parameters estimates from the true values in terms of error estimates, calculated among a set of 15,505 experiments with different noise realizations, shown for two methods: the global fit method and the transit-by-transit method.Top row: the high-SNR case.In this case, the global fit method generates parameter estimates that nicely fit a Gaussian, while the transit-bytransit method generates a broader The black dashed indicates a normal distribution with unity standard deviation.Bottom row: the low-SNR case.In this case, both distributions deviate from a Gaussian; the transit-by-transit method generates a much broader distribution, for which the peak is almost erased, while the global fit remains well-centered.

Figure 5 .
Figure 5.The mean deviation of the parameter estimates from the true value and the normalized deviation in error estimate as a function of transit SNR.Top row: the dots indicate the mean parameter estimate, and the error bars indicate the spread of the parameter estimates.The black dashed line indicates the true parameter value.Bottom row: the square root of the mean square of the deviation-to-error-estimate ratio.Color legend: blue=transit-by-transit, yellow=global flux fit, purple=transit-by-transit with a fixed ω value, green=global flux fit with a fixed ω value.See the text for further explanations.