Gaussian process (GP) regression is a nonparametric Bayesian approach that has been used successfully in various astronomical domains, especially in time-domain astronomy. The most common applications are the smoothing of data for interpolation and the detection of periodicities. The ability to create unbiased data-driven models without a predefined physical model can be a major advantage over conventional regression methods. Prior knowledge can be included by setting boundary conditions or constraining hyperparameter values, while unknown hyperparameters are optimized during the conditioning of the model. We have adapted and transformed previous approaches of GP regression and introduce three new applications for this regression method, especially in the context of stellar occultations: the modeling of occultation light curves, the correction of public JPL ephemerides of minor planets based on publicly available image data of the Zwicky Transient Facility, and the detection of natural satellites. We used data from observations of stellar occultations to validate the models and achieved promising results in all cases, and thus we confirmed the flexibility of GP regression models. Considering various existing use cases in addition to our novel applications, GP regression can be used to model diverse data sets addressing a wide range of problems. The accuracy of the model depends on the input data and on the set boundary conditions. Generally, high-quality data allow the usage of loose boundary conditions, while low-quality data require more restrictive boundary conditions to avoid overfitting.
1.Introduction
Today, the astronomer can pick from various regression models and algorithms to fit a model to a data set. For big-data problems, the choice of the regression model is usually time driven. Due to the large amount of data, even simple models can lead to the same result as more sophisticated models. Hence, a computationally efficient model is the best choice (Ivezić et al. 2020). For small data sets (<1000 samples), however, different regression models can lead to different results. This instantly raises the question of which model is best for the specific use case. In general, a regression problem has the form of
In this equation, M represents the regression model used to estimate the expected value of the dependent variable y, given the independent input variable x and the model parameters θ. The difference between y and its expected value M(θ, x) is captured by the term 
, which accounts for the inherent noise in measurements. It is important to note that the regression model M relies on user-defined assumptions.
For instance, the regression model may be an analytic parametric function, where the selection of the model parameters θ is guided either by prior knowledge or by optimization through a metric that measures the goodness of fit. In either case, selecting the parametric function assumes that the chosen model parameters adequately describe the data set. While including prior knowledge is important to keep the model physically plausible and within well-understood boundary constraints, some degrees of freedom allow the model to account for hidden features. More importantly, the result is an unbiased, data-driven model.
Gaussian process (GP) regression has become a popular choice to generate such models after it was extensively described by Rasmussen & Williams (2005). GP regression is a nonparametric Bayesian approach most efficient for small data sets, as the time complexity to create a GP model based on a data set with a sample size N is 
. The key feature of GPs is that they provide meaningful predictions including precise uncertainty estimations, even for small data sets (Ivezić et al. 2020). In astronomical applications, GP models have proven to achieve good results especially with time-domain data sets owing to their flexibility (Aigrain & Foreman-Mackey 2023). A common application is the modeling of light curves, and thereby (i) detecting trends, correlated noise, and systematic errors (e.g., Aigrain et al. 2016; Luger et al. 2021; Sun et al. 2023); (ii) interpolating the data to generate a smooth curve (de Oliveira et al. 2023); and (iii) filling gaps and making forecasts (e.g., Goumiri et al. 2022, 2023). Furthermore, the ability of GP models to detect periodicities has been used to model radial velocities (e.g., Camacho et al. 2023; Stock et al. 2023) and to detect features in the residuals of existing models (Czekala 2018).
One area in astronomy to which GPs have not yet been applied are stellar occultations. Stellar occultation observations, during which the flux of light from a distant star gets temporarily blocked by a body crossing the line of sight to the observer, can provide highly accurate constraints for the physical characterization of bodies in our solar system, surpassing any other means except direct spacecraft encounters (e.g., Elliot et al. 2010; Braga-Ribas et al. 2014; Buie et al. 2020; French et al. 2023; Sicardy 2023).
This paper is structured as follows. We provide a brief overview of the GP regression algorithm in Section 2. We transfer the approaches of existing GP models to create new models and describe three new applications specifically in the context of stellar occultations in Section 3. We present our results in Section 4, and we conclude with a discussion in Section 5.
2.Gaussian Process Regression
Various descriptions of GPs exist, and Rasmussen & Williams (2005) in particular provide a detailed description in their book. In this work, we focus on the advantages and disadvantages of GP regression instead of the math behind the algorithm. In general, a GP is completely defined by its mean function m(x) and covariance function k(xi , xj ) as
where x is the independent input variable, such as the time of an observation, and M is the model defined by a GP. The GP model is not a single predictive function, but it is a multivariate Gaussian distribution over functions as shown in Figure 1.
The covariance function shapes the model functions and thereby the Gaussian distribution over functions. It is the most important term of the GP and is mostly referred to as the kernel. The kernel defines the similarity between the data points, that is, the covariance. Various predefined kernels exist, all designed for specific use cases, such as modeling periodic or linear signals. The parameters of a kernel are called hyperparameters since they control not a single function but the distribution over an ensemble of functions. For example, the Gaussian kernel (GK) is defined as
where the length scale L is the hyperparameter of the kernel. The common optimization approach of the hyperparameters θ of any kernel is the maximization of their posterior probability defined as
by simply maximizing the marginal likelihood p(y∣θ, M). For GPs, this marginal likelihood is a Gaussian distribution 
itself, defined as
allowing a solution in closed form. This method is called type 2 maximum likelihood (ML-II) optimization (Rasmussen & Williams 2005). The prior kernel matrix Kxx consists of the covariances between all observed values. For convenience, the logarithmic marginal likelihood is used for optimization in ML-II. Its first term measures the ability of the model to describe the observed values y. The second term is a complexity penalty, as it is only dependent on the covariance matrix Kxx . The last term is a normalization constant to guarantee that the marginal likelihood is a probability distribution, only depending on the number of samples n. The usage of the marginal likelihood alone is sufficient to optimize and compare kernels, as it natively penalizes overfitting. The hyperparameters θ that maximize the marginal likelihood are used to define the conditioned model (Duvenaud 2014). With the conditioned model, the final prediction of the future value y* given the observed values y is the conditional probability
which again is a Gaussian distribution itself. The best approximation of the true value y* is 
, and its uncertainty is 
(Rasmussen & Williams 2005). The implementation of this GP regression algorithm is available in the high-level GaussianProcessRegressor 6 framework of the Python package scikit-learn (Pedregosa et al. 2011).
2.1.Model Generation
We use the scikit-learn package and its predefined kernels to create the models. We call the predefined kernels "basic kernels" (BKs), as they can only account for a single, specific purpose. We use five different BKs as listed in Table 1. The GK, 7 also called the radial basis function kernel, enables smooth interpolation between data points. The GK is the most commonly used kernel, as it can fit any arbitrary signal while depending only on one hyperparameter. The downside of having only one hyperparameter—the length scale L (see Equation (3)), comparable to the bandwidth of a filter—is that the kernel cannot be used alone to extrapolate data. Forecasting measurements beyond the optimized length scale would always return the mean value.
Table 1.Predefined Basic Kernels in the scikit-learn Package and Manually Compounded Regression Kernels
| Type | Name | Label | Structure |
|---|---|---|---|
| BK | Gaussian kernel | GK | ⋯ |
| BK | Periodic kernel | PK | ⋯ |
| BK | Linear kernel | LK | ⋯ |
| BK | Constant kernel | CK | ⋯ |
| BK | Noise kernel | NK | ⋯ |
| RK | Polynomial kernel | ⋯ | ∏ LK |
| RK | Locally periodic kernel | ⋯ | GK × PK |
Note. The first section lists the BKs modeling the signal. The second section contains the NK, which is a BK solely used to avoid overfitting. For better understanding, we use self-explanatory names from Duvenaud (2014) for the predefined BKs instead of the names given in the scikit-learn package. The third section shows two exemplarily created RKs that are based on the predefined BKs. The RK structures were also adapted from Duvenaud (2014).
Download table as: ASCIITypeset image
If the desired application involves forecasting data, other kernels are mandatory (Duvenaud 2014). Three kernels are suitable in particular: (i) the periodic kernel (PK), 8 also called the exponential-sine-squared kernel, which can detect and model periodic signals; (ii) the linear kernel (LK), 9 also called the dot product kernel, which models linear relationships in the data; and (iii) the constant kernel (CK), 10 which is only called a kernel for consistency, as it is simply a constant value. For example, the CK can be used to account for a fixed offset of the model. The last BK is the noise kernel (NK) 11 —also called the white kernel—as it is based on a Gaussian white-noise model; it is used to account for noise in the data and to avoid overfitting.
In addition, we can arbitrarily combine the BKs to create more complex kernel structures. The composition of kernels follows regular algebraic rules. Multiplications are executed before additions, and added terms are independent of their order (Duvenaud 2014). We call the composition of multiple BKs a regression kernel (RK), as the kernel that we will eventually use to model the data is always a composition of the NK and at least one additional BK. The choice of the RK structure is crucial to model the data correctly. Thus, the only terms of the RK that should be fixed are those that describe already-known features in the data. If no features are known and the RK cannot be defined in advance, we develop an automated model generation procedure to find the best RK structure.
2.2.Automated Model Generation
The automated model generation procedure generates the RK for a specific data set iteratively. We use the ML-II optimization as provided by the Python package scikit-learn to condition and compare different RKs. Figure 2 outlines our developed model generation algorithm. We provide a list of free and fixed BKs, where fixed kernels represent the base of the RK, and the free kernels are the possible additions to the RK. Even if no features are known, we can be certain that our measurements contain noise, so we always set at least the NK as a fixed kernel. As we discuss later in Section 3, the list of fixed kernels varies depending on use case and availability of prior knowledge.
The RK-generating procedure starts with a primary loop iterating through the input list of free kernels 
. The first free kernel x from the input list will be added to the fixed kernel to form the basis for this iteration step. Subsequently, the hyperparameters get optimized by applying the ML-II method. The resulting new compounded kernel serves as the new RK. The procedure then enters a secondary loop, during which the next BK y will be determined, as well as whether it will be added or multiplied to the existing RK. For that, the algorithm iterates again through the input list of free kernels 
, adds and multiplies each free kernel as a term to the existing RK, and adds the newly constructed RKs to a new list M. In the case of four free kernels to choose from, the execution of the secondary loop thus results in eight new possible RKs. The best model z will be selected from the list M, which is the RK that returns the highest marginal likelihood 
during its conditioning. The secondary loop repeats until the marginal likelihood does not increase further by extending the RK with yet another free kernel. The RK with the highest marginal likelihood 
will be the final kernel of this iteration of the primary loop, and it will be added to the output list 
. The entire RK-generating process continues until the primary iteration through the list of free kernels 
is complete. Each iteration of the primary loop results in one RK with a corresponding 
. In the case of four free kernels, the output list 
contains four conditioned RKs. The RK that we select for our model is the one with the highest 
among the RKs in 
.
3.Applications of Gaussian Process Regression
GP regression is applicable to a variety of regression problems. Here we present three examples where GPs are a particularly good choice.
3.1.Light-curve Modeling
Our first example of GP regression is the modeling of light curves. GP regression has been successfully used by Goumiri et al. (2022) to predict gaps in light curves and by Aigrain et al. (2016) to detrend light curves. We combine these methods to model stellar occultation light curves by solar system bodies that have an atmosphere and to extract their key features, such as the half-light time at which the light curve reaches the 50% normalized stellar flux level—an important metric for astrometric referencing and the global evolution of an atmosphere (Elliot & Young 1992). We use the half-light times to mark the beginning and the end of the occultation, that is, the immersion and the emersion, respectively, because they are easily comparable to previously observed occultations of the same body. The standard approach is to fit a predefined atmospheric model to the measured light-curve data. To create such a model, all effects that govern the light curve—such as diffraction due to the atmosphere's composition and vertical structure, or the potential presence of aerosols—must all be well understood, considering that erroneous or insufficient assumptions will propagate into the results of the modeled light curve. Planetary atmospheres can have complex structures and properties, and often only limited information is available to model the occultation light curve unambiguously. As an alternative, we propose to model the light curve using the GP regression method instead, as no prior assumptions are necessary to fit the model.
For our light-curve model, we are interested in the interpolation between the data points, not the extrapolation outside the range of measurements. Thus, there is no need for the automated model generation approach previously described in Section 2.2, and it is sufficient to use the composition of the GK and the NK as our RK instead. The first step is to fit the RK to the entire data set; this fit serves as a first estimate of the modeled light curve. We can derive the upper and lower baselines from this model and isolate the time interval during which the occultation has occurred. If the light curve shows a trend in its upper baseline, it can be removed in an optional second iteration. This is achieved by fixing the hyperparameter length scale L of the GK to the number of samples of the upper baseline, as we are aiming at modeling global trends instead of local features. We can subsequently renormalize the data with this optional "trend model." In the final iteration, a GP model is fit to the data within the time interval of the occultation. This optimized model is then applied to the entire light curve.
The modeled light curve allows us then to extract the half-light times with subsampling precision and to properly derive robust uncertainties from the GP model fit of the estimated times from measurements that will always exhibit noise.
3.2.Ephemeris Correction Model
The next example is more complex, and it requires some preliminary steps before we can eventually generate an ephemeris correction model (ECM) with GP regression. An ECM is used to estimate the ephemeris offset to a reference ephemeris at an epoch of interest. We use the publicly available ephemerides provided by the Jet Propulsion Laboratory's Solar System Dynamics group (JPL SSD) 12 and JPL's Navigation and Ancillary Information Facility (NAIF) 13 as the reference for our ECMs.
Our motivation to create ECMs is to improve the prediction of stellar occultations; a successful prediction depends on the precise knowledge of the position of the star and the solar system object (SSO). With the release of the Gaia star catalogs (Gaia Collaboration et al. 2023, 2016), our knowledge about the positions and proper motions of stars has been significantly improved, leading to more accurate predictions of stellar occultations (Ferreira et al. 2022). However, ephemerides with milliarcsecond-level accuracy are required to accurately predict at which locations a stellar occultation can be observed, and ephemerides of most SSOs—particularly fainter and/or more distant objects—still lack that level of accuracy. Aside from their low brightness, trans-Neptunian objects (TNOs) have only been observed over a fraction of their long-period orbits; their ephemerides can have significant errors. These systematic offsets to the JPL ephemeris must be modeled and taken into account to successfully predict and observe an occultation event (e.g., Ortiz et al. 2017; Person et al. 2021, 2013).
The JPL ephemerides are predominantly based on observations from the Minor Planet Center (MPC) database. 14 Observers measure positions of known SSOs from their images themselves and report them to the MPC; images are not archived, preventing it from revisiting reported positions or rereducing pre-Gaia observations. The MPC publishes all observations in its database, listing the reported position, instrumental magnitude, and observation date. Quality metrics—such as astrometric uncertainties—that could be used to weight or omit measurements when calculating orbit solutions are not available; fit residuals can be used a posteriori to isolate problematic observations (Farnocchia et al. 2015).
We thus need additional image data that would allow us to conduct position measurements that would consequently help us to improve ephemerides of interest. The usual approach is either to predict the ephemeris offset after conducting dedicated, long-term observations of SSOs—often collecting data for years in advance, and particularly when coordinating a larger observing campaign, as frequently as resources allow in the weeks leading up to an event—or to fit an entirely new orbit solution (e.g., Desmars et al. 2015). However, both approaches require a significant investment of time and resources, as well as a significant time allocation at suitable telescopes, which often have limited, institution-dependent access.
Instead, we use the public data of the Zwicky Transient Facility (ZTF) wide-field sky survey, acquired with the 48-inch Samuel Oschin Telescope at Palomar Observatory (Masci et al. 2018). During this work, calibrated images of the ZTF survey got published bimonthly at the Infrared Science Archive (IRSA); the release schedule was altered to 4-month intervals at the end of 2023. 15 The ZTF survey started in 2018 March, and as of today, more than 50 million single-exposure images are available. The Moving Object Search Tool (MOST) at IRSA allows us to filter the database and to easily access the images containing the SSO whose ephemeris we want to improve.
We use the TNO and dwarf planet (134340) Pluto as an example to generate an ECM. We aim at improving Pluto's public JPL ephemeris by using a data set with superior astrometric accuracy compared to data from the MPC, and we create such a data set from publicly available ZTF images as shown in Figure 3. As of 2023 December, the entire ZTF data set of Pluto observations consists of 233 images, but only 30 of these observations have been reported to the MPC. The submission of ZTF observations to the MPC is part of an automated pipeline, whose scientific purpose is to detect transient events instead of observing known SSOs. Hence, only a small fraction of the observations gets eventually reported to the MPC (F. Masci 2024, personal communication). For all observations, the calibrated images and their auxiliary tables listing the extracted and plate-solved sources, identified with SExtractor (Bertin & Arnouts 1996) and SCAMP (Bertin 2006), are available in FITS format at IRSA. We exclude images taken at high air mass or during bad seeing conditions, defining a threshold of S < 2'', where S is the FWHM of the stellar images assuming a Gaussian fit. We also exclude all images in which field stars blend with Pluto, impeding astrometric measurements of high accuracy. With 
(Wilson 2001), and with Z being the zenith angle, we set the minimum angular distance between nearby stars and Pluto to 
to force increasing minimum angular distances at lower elevations.
3.2.1.Astrometry
The astrometric solution provided by the ZTF is referenced to Gaia DR1, which does not include proper motions. Thus, it has a systematic error in right ascension (R.A., or α) and declination (decl., or δ) that increases as we move farther away from the Gaia DR1 epoch (J2015.0). To reduce this systematic error, we calculate a new astrometric solution referenced to Gaia DR3—now considering proper motions—using the MATLAB astronomy and astrophysics toolbox (MAAT) developed by Ofek (2014). The MAAT provides source extraction, centroiding, and plate-solving algorithms to compute precise and robust astrometric solutions. The performance of the algorithms was tested and validated with ZTF data. The reference stars used for plate solving are nonsaturated stars distributed over the entire image, and the rms of the reference star residuals describes the overall accuracy of the astrometric solution (Ofek 2019).
The astrometric uncertainty of an image can have various reasons, such as centroid offsets caused by optical aberrations. While a thorough astrometric analysis of each image is beyond the scope of our work, we estimate the uncertainty of Pluto's measured position by assessing the residuals of calculated positions of nearby stars with comparable brightness. We refer to this subset of stars as "evaluation stars." Per our definition, evaluation stars are located within a radius of 1000'' around Pluto and within a magnitude range of ΔmV = ±1 mag to ensure a sufficient sample size for statistical inference. We use the transformation polynomial from Riello et al. (2021) and Carrasco & Bellazzini (2022) to convert stellar magnitudes measured in the Gaia mean G band to the Johnson V band, in which JPL Horizons reports the apparent magnitude of an SSO. Eventually, we map the measured positions of the evaluation stars to their respective epoch-propagated positions from Gaia DR3.
For an exemplary image, 16 we visualized the residuals between the evaluation star positions based on the MAAT plate-solve and the Gaia DR3 catalog positions in Figure 4 (right panel). For comparison, we also evaluated the original ZTF astrometric solution by computing the evaluation star residuals to the Gaia DR3 catalog positions (Figure 4, left panel). We fitted normal distributions to the residuals to quantify the astrometric errors based on the evaluation star sample. The mean values (μα , μδ ) of each normal distribution represent the systematic error, and the standard deviations (σα , σδ ) represent the random error of an astrometric solution. While the random error of MAAT solutions is comparable to ZTF solutions—indicating that our MAAT plate solves are of comparable quality—the systematic error is notably smaller, presumably due to referencing our solution to Gaia DR3. After verifying this trend throughout the entire ZTF data set of Pluto observations, we are using the positions measured by MAAT to generate our ECMs. We expect that the random error estimated from the evaluation star residuals is representative of the astrometric uncertainty of the measured target position (in this example, Pluto).
3.2.2.Input Data
To obtain a high-quality data set, we exclude observations with an estimated astrometric uncertainty exceeding 50 mas. This threshold represents a compromise, as filtering too aggressively leaves us with a low number of data points, but incorporating lower-quality measurements in the data set increases the uncertainties of the ECM. As more ZTF observations become available in the future, it could be an option to lower this threshold. Further, we exclude observations where less than 97 evaluation stars are available to estimate the astrometric uncertainty. A minimum sample size of 97 is required to estimate the 95% confidence interval with a ±10% margin of error (Guthrie 2020).
Compared to the data set of the MPC, we have only a small data set available, and we do not aim at computing a new orbital fit. Instead, we calculate the residuals (Δα, Δδ) between the observed positions (αZTF, δZTF) and the publicly available JPL ephemerides (αJPL, δJPL) for the moving object as
The observed residuals (Δα, Δδ) and estimated measurement uncertainties (σα , σδ ) are the input data for the GP regression model. We use the automated model generation as described in Section 2.2 to find the best RK and to condition the model. We compute an individual RK for each residual, Δα and Δδ. Thus, the ECM contains two RKs to predict the offset (Δα*, Δδ*) to the reference ephemeris at a future epoch of an occultation.
3.3.Detection of Natural Satellites
Numerous TNOs have at least one reported natural satellite (e.g., Grundy et al. 2019; Noll et al. 2020). With the exception of high angular resolution observations, a ground-based optical telescope cannot resolve such a satellite; the measured position thus represents the center of light of the entire TNO system. If the automated model generation procedure detects a periodicity in the observed residuals—which is accounted for as a PK term in the final RK—and the detected periodicity coincides with the known orbital period of a satellite, we need to correct the observed center of light to successfully predict a stellar occultation of the main body. The optimized, conditioned PK describes the influence of the satellite, which we need to remove from the mean value to create a data set with corrected positions (Δαcorr, Δδcorr). This approach is adapted from Aigrain et al. (2016), who used GP regression to detrend light curves. After removing the periodic signal from our observations, we start a second iteration of the automated model generation process using the corrected positions.
If the object does not have a known satellite and the detected periodicity does not align with its reported rotation period, two conclusions can be drawn from the detection of a new periodicity: (i) the reported rotation period lacks accuracy owing to various factors, such as albedo variations across its surface or an irregular shape of the body, or (ii) at least one previously undiscovered satellite exists.
4.Results
To demonstrate our first example, we created a set of synthetic light curves to simulate different atmospheric profiles and noise levels (Figure 5). We use this set of synthetic light curves to validate the GP modeling approach to extract half-light times. Next, we use ZTF image data to create ECMs for (134340) Pluto, (136108) Haumea, and (2060) Chiron. Data from a stellar occultation by Pluto observed on 2022 August 23 UT act as a real-world test case to verify our Pluto ECM and to demonstrate the GP modeling approach on real light curves. For both Haumea and Pluto, our ECMs correctly model the orbital periods of their primary satellites purely based on the data at hand. For Chiron, we compare the modeled ephemeris offset to stellar occultations observed in 2018 and 2019.
4.1.Synthetic Light Curves
Synthetic light curves with known properties pose an ideal test case to validate the GP approach to model light curves. We simulated two occultation light curves, LC1 and LC2, assuming an isothermal profile of an atmosphere free of aerosols, a scale height of 65 km, and different half-light radii resulting in different half-light times. We sampled them at a frequency of 2 Hz and added Gaussian noise to each sample. Each test data set has a unique signal-to-noise ratio (S/N), defined as the arithmetic mean over the standard deviation of the baseline, and ranging from S/N = 1 to S/N = 800. Exemplary reconstructions of three synthetic light curves are shown in Figure 5; the estimated half-light times for a range of S/Ns are listed in Table 2. To visualize the error of the GP model as a function of S/N, we show the differences between the estimated half-light times and the simulated half-light time in Figure 6.
Table 2.Estimated Half-light Times from the GP Modeled Synthetic Light Curves
| LC1 | Immersion (mm:ss.sss) | Emersion (mm:ss.sss) |
|---|---|---|
| Truth | 03:25.489 | 04:54.511 |
| S/N 1 | 03:40.537 (±44.09 s) | 04:50.687 (±29.55 s) |
| S/N 3 | 03:27.035 (±9.391 s) | 04:52.531 (±6.544 s) |
| S/N 10 | 03:24.203 (±2.897 s) | 04:54.064 (±3.247 s) |
| S/N 20 | 03:25.991 (±1.598 s) | 04:53.969 (±1.598 s) |
| S/N 40 | 03:25.652 (±0.799 s) | 04:53.929 (±0.849 s) |
| S/N 80 | 03:25.731 (±0.400 s) | 04:54.214 (±0.400 s) |
| S/N 120 | 03:25.587 (±0.255 s) | 04:54.369 (±0.255 s) |
| S/N 160 | 03:25.397 (±0.190 s) | 04:54.504 (±0.195 s) |
| S/N 200 | 03:25.487 (±0.150 s) | 04:54.509 (±0.150 s) |
| S/N 400 | 03:25.486 (±0.080 s) | 04:54.494 (±0.080 s) |
| S/N 800 | 03:25.502 (±0.045 s) | 04:54.504 (±0.045 s) |
Note. The 1σ errors are given in parentheses.
Download table as: ASCIITypeset image
As expected, for lower S/Ns, the uncertainties of the half-light times are larger. For an S/N ≤ 40, the 1σ errors exceed 1 s. Furthermore, at lower S/N, the GP model tends to estimate half-light times slightly closer to the occultation center. Consequently, the immersion half-light time estimate occurs after the actual simulated half-light time, and the emersion half-light time estimate occurs before the simulated one. Nevertheless, half-light times could be estimated well within the 1σ errors in all cases. For an S/N ≥ 80, the 1σ errors are already below 0.5 s.
4.2.Pluto
Our ECM for Pluto is referenced to the latest JPL NAIF ephemerides, which are the JPL development ephemerides DE440 and the Pluto satellite ephemeris kernel version #58 (PLU058). We computed individual RKs for Δα and Δδ based on our data set of ZTF observations, from which we derived Pluto's astrometric position through plate solves obtained with MAAT, and estimated uncertainties individually for each observation as described in Section 3.2.1. For Pluto, our final data set consisted of 182 samples.
4.2.1.Detection of Charon
We started with the automated model generation procedure as described in Section 2.2. The resulting model and its 95% confidence region are shown in Figure 7. The oscillation of the observed residuals to Pluto's reference ephemeris is a clear indication that the observed center of light shifts owing to the superposition of fluxes from multiple bodies. The optimized RKs for Δα and Δδ both include a periodic term with a periodicity and 1σ uncertainty of p ≈ 6.387 ± 0.008 days. The detected periodicity coincides extremely well with Charon's orbital period of 6.387 days determined from data from the New Horizons mission (Stern et al. 2018).
In theory, not only Pluto's primary satellite Charon but all of its five moons contribute fluxes that shift the observed center of light, and their orbital periods should thus be detectable by the GP model as well. In reality, however, the astrometric uncertainty of the measured center of light is too large to detect secondary orbital periods of smaller, fainter satellites. To estimate the limits of the GP model's ability to detect periodicities in data sets, we created another synthetic data set, resembling the same sample size and random, Gaussian distribution of astrometric uncertainties of our existing ZTF data set. The (Δα, Δδ) observed with respect to the reference ephemerides caused by a single satellite can be described by a sine function. The sine function is dependent on the brightness difference and angular separation between Pluto and the satellite. The sine function has a periodicity set to the orbital period of the satellite, which is 6.387 days for Chiron and 38 days for Hydra. The resulting center of light is the sum of the sine functions. Using these two sine functions, we sampled synthetic measurements of (Δα, Δδ) at the epochs of the ZTF data set. Finally, we added Gaussian noise to each sample. The results are shown in Figure 8. Pluto and Charon have a maximum angular separation of 1000 mas and a ΔmV = 2 mag magnitude difference. The simulation concurs with the actual data that the model can detect Charon's orbital period. For Pluto's second-brightest and outermost satellite, Hydra, the simulation confirms that a magnitude difference of ΔmV = 6 mag results in a periodic variation that is undetectable by the GP model.
To predict a stellar occultation by Pluto itself, we correct the observed center of light as described in Section 3.3 to obtain Pluto's estimated position. Consequently, the remaining periodicity for the center-of-light-corrected ECM is the interval between two oppositions of Pluto and Earth, which is p = 1.004 yr as illustrated in the right panels of Figure 7.
4.2.2.Occultation on 2022 August 23 UT
We consider Pluto's stellar occultation on 2022 August 23 UT to validate our ECM, which predicted an offset to the JPL DE440 and PLU058 ephemeris of (Δα, Δδ) = (6 mas, 4 mas) at the time of the occultation. As shown in Figure 9, the predicted shadow path was over Northern Africa and Europe. The event was successfully observed using the fast-readout camera RISE (Steele et al. 2008) at the 2 m Liverpool Telescope (LT; Steele et al. 2004) at Roque de los Muchachos Observatory on La Palma, Spain, and using a QHY174M-GPS camera at the Trebur 1.2 m Telescope (T1T) at Michael Adrian Observatory in Trebur, Germany. RISE uses a fixed "V+R" filter; the observations at the T1T were conducted without a filter and with an aggressive focal reduction of about ×0.37. Table 3 lists the latitude and longitude of both stations. Full details of these observations and their discussion in the context of the recent evolution of Pluto's atmosphere can be found elsewhere (A. Sickafoose et al. 2023).
Table 3.Location of the Sites That Observed the Stellar Occultation by Pluto on 2022 August 23 UT
| Location | N Latitude | E Longitude | Altitude |
|---|---|---|---|
| (° ' '') | (° ' '') | (m) | |
| LT | 28 45 44.8 | −17 52 45.2 | 2363 |
| T1T | 49 55 31.6 | 08 24 41.1 | 103 |
Download table as: ASCIITypeset image
LT provides readily reduced (=bias-, dark-current-, and flat-field-corrected) data to the observer. 17 At the T1T, dark frames of equal exposure time and sensor temperature were acquired right after the observation concluded; the images were then debiased and dark corrected by subtracting a corresponding median dark frame. A flat-field correction could not be applied.
Light curves were extracted through an iterative process using aperture photometry, which is implemented in our own Python tool FLUXI based on the publicly available Python packages astropy (Astropy Collaboration et al. 2022), numpy (Harris et al. 2020), and photutils (Bradley et al. 2021). In a first pass over the image series, centroids of all relevant stars are measured in each image, and their average relative positions to each other are determined. That enables us to construct a mask of apertures with fixed locations relative to each other. In the passes that follow, aperture radii are systematically varied to conduct a curve-of-growth analysis (see, e.g., Howell 1989), determining the aperture radius that enables the highest S/N for optimal light-curve extraction. Rather than placing apertures individually on each star based on its centroid measurement in each image, the constructed aperture mask is fitted via least-squares optimization to the centroid list of each image, ensuring consistent aperture placement. The flux f measured with the optimal aperture radius is then normalized to fnorm using the arithmetic mean of the baseline before and after the event. A more detailed discussion of our photometric reduction approach can be found in K. Schindler et al. (2024, in preparation).
Due to refraction in Pluto's atmosphere, residual flux from the occulted star still reaches the observer during the occultation, in addition to flux from Pluto and its satellites and other sources of spurious signal. To be able to measure the residual flux from the occulted star, we estimate the so-called "background fraction," defined as 
, by measuring the flux ratio r = fBody/fStar from images acquired before and/or after the occultation while both have sufficient angular separation. We then use 
to calculate the background-corrected, normalized flux 
, which we are using to create a GP model based on a GK and NK over three iterations as described in Section 3.1.
For the LT data set, we utilized images acquired the night before and 2.5 hr after the event to conduct a linear fit on bf versus air mass. This allowed us to obtain bf = 0.897 ± 0.006 at the air mass of the occultation. Due to the extremely low elevation and unfavorable weather conditions, T1T was unable to observe Pluto separated from the star, so we applied the same background fraction as derived from the LT data.
With the GP regression method, we can model the light curve from both sites without an atmospheric model. The modeled light curves are shown in Figure 10. The T1T light curve has a very low S/N, causing immersion and emersion time uncertainties of about 20 s. This required us to constrain the shadow radius to Pluto's half-light radius R = 1213.55 km to reconstruct the shadow path. This half-light radius was derived from a comparable stellar occultation by Pluto in 2022 June (A. Sickafoose et al. 2023). The resulting circle fit is shown in Figure 11, and we derived an offset to the JPL ephemeris of 
, based on the coordinates of the circle's center in the sky plane 
. The measured and predicted occultation parameters are listed in Table 4.
Table 4.Measured and Predicted Times for the 2022 August 23 Pluto Occultation at the Trebur 1.2 m Telescope and the 2 m Liverpool Telescope
| T1T | Immersion | Emersion | Δta | ρb |
|---|---|---|---|---|
| (UT) | (UT) | (s) | (km) | |
| Observed | 21:19:46.4 ± 21.0 | 21:20:33.5 ± 17.1 | 47 | 1162 |
| ECM c | 21:19:47 | 21:20:24 | 37 | 1158 |
| JPL d , e | ⋯ | ⋯ | ⋯ | ⋯ |
| LT | Immersion | Emersion | Δta | ρb |
|---|---|---|---|---|
| (UT) | (UT) | (s) | (km) | |
| Observed | 21:21:28.8 ± 4.4 | 21:23:19.8 ± 5.0 | 111 | 498 |
| ECM c | 21:21:27 | 21:23:20 | 113 | 494 |
| JPL d | 21:21:23 | 21:23:07 | 104 | 550 |
Notes.
a The duration of the observed occultation. b The closest approach of the observed chord to the center of the shadow. c Prediction based on our ECM with Pluto radius set to 1213.55 km. d Prediction based on the JPL ephemeris for Pluto with its radius set to 1213.55 km. e Trebur was not in the shadow path of the prediction.
Download table as: ASCIITypeset image
4.3.Haumea
We use Haumea's JPL orbital solution #116 in combination with JPL DE440 as the reference ephemeris for our Haumea ECM. Again, we exclude observations with a random astrometric error exceeding 50 mas. Our final data set for Haumea consists of 467 data points from the ZTF between 2018 May and 2023 May, marked by red circles in Figure 12. The estimated 1σ uncertainties as described in Section 3.2.1 are indicated by the error bars. The dark-blue line represents the mean prediction of the created ECM, and the light-blue area represents its 95% confidence region. With the automated model generation, we found for both Δα and Δδ models an optimal RK structure consisting of an LK, PK, and NK. The overall linear trend of the ephemeris offset is modeled by the LK.
The optimum value and its 1σ uncertainty for the hyperparameter of the PK, the periodicity p, was found at p ≈ 49.0 ± 0.5 days. This coincides with the orbital period of Haumea's primary moon Hi'iaka, previously estimated as p ≈ 49.3 days (Ortiz et al. 2017). As for the Pluto system, the measured ZTF positions represent the center of light of the entire Haumea system. Haumea's secondary moon, Namaka, has a maximum angular separation of 700 mas from Haumea and a ΔmV ≈ 4.5 mag difference in apparent magnitude. Thus, Namaka's influence on the measured center of light is less than the measurement noise. Consistent with our simulation for the Pluto system as shown in Figure 8, we can only detect the orbital period of Haumea's primary satellite.
To predict stellar occultations by Haumea, we correct the astrometric positions from ZTF observations by removing the periodic influence of the primary satellite, analogous to our Pluto ECM. The linear trend of the ephemeris offset remains in the ECM. Unfortunately, we could not find any stellar occultations by Haumea observable from Earth until 2024 December.
4.4.Chiron
For Chiron, more than 700 ZTF observations are available at IRSA, but due to its faintness (mV ≈ 18 mag) and thus low S/N in ZTF images, the estimated astrometric uncertainties are much bigger than in the previous two examples. The usage of the entire set of observations would result in occultation predictions with large error bars. Yet our previously applied uncertainty threshold of 50 mas would exclude most of the data. We thus exclude only observations with an astrometric uncertainty exceeding 70 mas to create a sufficiently large data set. This provides us with 138 ZTF observations between 2018 May and 2023 May. JPL ephemeris #149 and JPL DE440 are the reference ephemerides for our Chiron ECM.
We revisit occultations by Chiron on 2018 November 28 UT observed by A. A. Sickafoose et al. (2023), on 2019 September 8 UT observed by Braga-Ribas et al. (2023), and on 2022 December 15 UT observed by Ortiz et al. (2023) to compare (a) a prediction based on the JPL ephemeris alone, (b) a prediction based on our ECM, and (c) the reconstructed event geometry ("postdiction").
For the occultation on 2018 November 28 UT, our ECM predicts an offset to JPL #149 of (Δα, Δδ) = (24 mas, − 37 mas). The prediction based solely on Chiron's JPL ephemeris is provided in Figure 13(a); the prediction made with our ECM is shown for comparison in Figure 13(b). We used the immersion and emersion times from the 74-inch telescope at the South African Astronomical Observatory (SAAO) as provided by A. A. Sickafoose et al. (2023) to calculate the chord length. Since we have only one chord for this occultation, we cannot estimate Chiron's size from the data, so we use the equivalent radius of R = 105 km (Lellouch et al. 2017) to postdict the shadow path. Based on the equivalent radius and the single SAAO chord, we postdict an ephemeris offset of (Δα, Δδ) = (14 mas, − 36 mas) at the time of the occultation. Figure 13(c) shows the postdicted shadow path for the 2018 November 28 UT event.
We tried to postdict the shadow path of Chiron's occultation on 2019 September 8 UT in the same way, but the immersion and emersion times of the stations provided in Braga-Ribas et al. (2023) revealed some inconsistencies. We thus abandoned the effort to estimate the shadow size and to postdict the shadow path for this event. Although we cannot reconstruct the location of the actual shadow path, the shadow path predicted with our ECM, shown in Figure 14(b), covers all observing sites that reported a positive observation for this event. In comparison, the prediction solely based on the JPL ephemeris has a more northerly shadow path as seen in Figure 14(a): none of the stations are located inside the predicted path, and three out of the four stations that reported positive observations are even located outside the prediction's 3σ confidence region.
For the 2022 December 15 UT event, Ortiz et al. (2023) did not provide any immersion and emersion times. Once again, we cannot generate a postdiction and instead only compare the JPL and ECM predictions based on the reported positive observations of the event. The only solid-body observation was reported by the Wise Observatory. The Kottamia Astronomical Observatory reported an observation of a secondary feature indicating additional material around Chiron, and the 
National Observatory reported a negative observation. Figure 15(a) shows the prediction based on the JPL ephemeris with all observing stations located outside the predicted shadow path. The prediction based on our ECM has a more southerly shadow path, resulting in predicted solid-body observations by the Wise Observatory and the Kottamia Astronomical Observatory as shown in Figure 15(b).
5.Discussion
The GP regression modeling was successful for all created synthetic light curves. We used different atmospheric conditions to simulate synthetic occultation light curves with varying levels of random noise and validated the GP modeling approach by fitting a GP model to each test data set. Without including prior assumptions in our models, we could estimate the half-light times well within 1σ errors in all cases. We thus expect to derive the half-light times from real occultation light curves with the same precision. The RK used for the light-curve models can be compared to a smoothing function, where the optimized hyperparameter represents the degree of smoothing. This model can be applied to data sets where unbiased feature extraction is more important than the underlying model itself. In our case, we are only interested in the half-light times to reconstruct the geometry of an occultation event.
Subsequently, we utilized observations of an occultation by Pluto on 2022 August 23 UT to model the extracted light curves, derive the half-light times, and reconstruct the shadow path. Direct comparison of our modeled light curves with atmospheric fitted light curves showed consistency within 1σ error bars (A. Sickafoose et al. 2023). We then compared our ECM-based predictions with the JPL ephemeris. All half-light times predicted using our ECM fell within the 1σ errors of the half-light times derived from light curves, whereas the prediction using solely the JPL ephemeris fell within the 2σ uncertainty interval of the measurement (see Table 4). Based on a circle fit to the LT and T1T chords, our ECM-predicted ephemeris offset differs by less than 1 mas from the observed ephemeris offset, indicating that our ECM improved the JPL ephemeris successfully.
Next, we demonstrated the abilities of GPs to detect orbital periods of natural satellites. Pluto and Charon are both well above the limiting magnitude of the ZTF survey, and we were able to precisely measure Charon's orbital period from the ZTF data. We did the same for the Haumea system; however, Haumea and its primary moon Hi'iaka are much fainter than Pluto and Charon. Our best-fitting model has a periodicity p ≈ 49.0 ± 0.5 days, compared to p ≈ 49.3 days as stated by Ortiz et al. (2017).
As with all modeling approaches, the uncertainty of a GP model depends on data quality. This is also evident from the ECM for Chiron; its apparent magnitude of mV ≈ [18 − 18.5] mag is closer to the limiting magnitude of the ZTF, whose median 5σ (i.e., S/N = 5) limiting magnitudes were determined as 
, 
, and 
over one lunation (Bellm et al. 2019). While Chiron is well above the detection threshold, its low S/N limits the accuracy at which its position can be measured. As a result of the large uncertainties of Chiron's position, the uncertainties of predictions made with our ECM exceed those of a prediction based on the JPL ephemeris. As uncertainty correlates with sample size, we anticipate that model uncertainty would decrease somewhat by adding future ZTF observations. Still, we showed that occultation predictions with our ECM are more precise than using solely the JPL ephemeris of a body. For Chiron's occultations on 2018 November 28 UT, 2019 September 8 UT, and 2022 December 15 UT, the predicted shadow path based on the GP model accurately covered all observing sites that were successfully observing the occultation. In comparison, the predicted shadow path based on the JPL ephemeris was shifted by about 1000 km for the occultations in 2018 and 2019 and by about 500 km for the occultation in 2022. For the occultation in 2019, parts of the reconstructed shadow path were even located outside of the 3σ error interval of the JPL ephemeris, suggesting that the uncertainties provided by JPL are underestimated.
While the uncertainties of GP models increase quickly with increasing random errors in the underlying data, we always found the best estimate to be close to the true value, despite its sometimes large confidence intervals. It is crucial for good results to remove all known systematic errors from the data before computing the final model; being entirely data-driven, GPs have no feature to compensate for such bias. In our work, we did that by fitting a detrending model to light curves before modeling the actual occultation signature and through our own astrometric rereduction of ZTF images with MAAT in the Gaia DR3 reference frame before creating our ECMs.
The use of higher-quality data improves the performance of GP regression, as well as other regression models. However, while conventional regression methods could still produce similar results with the right boundary conditions, we believe that GP regression is a superior method, as models do not depend heavily on boundary conditions.
In addition, our research highlights the potential of mining publicly available data. We tried to achieve the highest astrometric accuracy based on the input data available for our ECM, but as of today, the list of possible targets to create ECMs from public data is limited. The 48-inch Samuel Oschin telescope utilized by the ZTF to acquire 30 s exposures allows us to get meaningful astrometric positions for targets of mV ⪅ 18.5 mag. For even fainter targets, we estimated the astrometric 1σ uncertainties to exceed 100 mas for the majority of observations. The usage of low-quality input data then requires additional boundary conditions of hyperparameters of the RK. For example, by allowing the PK to model high frequencies, it is likely to model an artificial periodicity from very noisy measurements, and the resulting ECM would not provide a meaningful improvement to the public JPL ephemeris anymore.
We expect the effectiveness of GP regression models to continue to improve as new large surveys emerge. As the number of publicly available data sets with deeper limiting magnitudes increases, the list of targets that we can probe for natural satellites or for which we can improve their ephemeris grows. Of particular interest could be additional data releases from Pan-STARRS and, of course, the LSST survey, which is anticipated to commence in 2025.
6.Conclusions
GP regression is becoming an increasingly important tool for astronomers because of its flexibility to describe multiple data sets and the ability to create meaningful models even with small data sets. The results of GP regression are calculated as a normal distribution, which means that the modeled expected values are natively provided with a Gaussian uncertainty. The magnitude of the modeled uncertainty strongly depends on the quality of the input data. To keep the uncertainty intervals sufficiently small to draw meaningful conclusions, it is helpful to preprocess the data to remove systematics. If heteroscedastic data like our astrometric position measurements from ZTF data are used, we recommend excluding samples with significant random errors to ensure a high-quality model fit.
In addition to the input data, the accuracy of the GP model depends on its RK. Prior knowledge about the data set can be included in the model either by manually adding terms to the RK or by fixing hyperparameter values. The radial basis function (GK) is among the most flexible predefined kernels to model data, but it cannot be used alone to extrapolate data, such as to predict future ephemeris offsets. Instead, it is a powerful kernel to interpolate data, such as to model light curves. If no prior knowledge is available to the astronomer, we present an automated approach that iteratively finds the best-fitting model for the data.
In general, GP regression models are flexible and can be fitted to diverse data sets to solve both interpolation and extrapolation problems. As GP regression is a data-driven method, the key to a precise prediction is to fully understand the data and its limitations and to include prior knowledge—if available—as boundary conditions in the RK.
Given the advantages of GPs over conventional regression methods for certain data sets, we expect that their use will continue to expand across various astronomical domains in the future.
Acknowledgments
We thank Johannes Ohlert, Jürgen Wolf, Felipe Braga-Ribas, and Jose Luis Ortiz for their observations of stellar occultations, which assisted in validating the method presented in this work.
The Liverpool Telescope is operated on the island of La Palma by Liverpool John Moores University in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofisica de Canarias with financial support from the UK Science and Technology Facilities Council. The Trebur 1.2 m T1T Telescope is operated by the Astronomie Stiftung Trebur.
SOFIA, the Stratospheric Observatory for Infrared Astronomy, is a joint project of the Deutsches Zentrum für Luft- und Raumfahrt e.V. (DLR; German Aerospace Center, grant: 50OK0901, 50OK1301, 50OK1701 and 50OK2002) and the National Aeronautics and Space Administration (NASA). It is funded on behalf of DLR by the Federal Ministry for Economic Affairs and Climate Action based on legislation by the German Parliament, the State of Baden-Württemberg, and the University of Stuttgart. SOFIA activities are coordinated on the German side by the German Space Agency at DLR and carried out by the German SOFIA Institute (DSI) at the University of Stuttgart, and on the U.S. side by NASA and the Universities Space Research Association (USRA).
A portion of this work was supported by the MIT-Germany University of Stuttgart Seed Fund (misti.mit.edu) and by the NASA Solar Systems Observations grant 80NSSC21K0432.
Results presented in this work are based on observations obtained with the 48-inch Samuel Oschin Telescope at the Palomar Observatory as part of the Zwicky Transient Facility project. ZTF is supported by the National Science Foundation under grant Nos. AST-1440341 and AST-2034437 and a collaboration including current partners Caltech, IPAC, the Weizmann Institute for Science, the Oskar Klein Center at Stockholm University, the University of Maryland, Deutsches Elektronen-Synchrotron and Humboldt University, the TANGO Consortium of Taiwan, the University of Wisconsin at Milwaukee, Trinity College Dublin, Lawrence Livermore National Laboratories, IN2P3, University of Warwick, Ruhr University Bochum, and Northwestern University and former partners the University of Washington, Los Alamos National Laboratories, and Lawrence Berkeley National Laboratories. Operations are conducted by COO, IPAC, and UW.
This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.
Last but not least, this research has made use of NASA's Astrophysics Data System Bibliographic Services.
Software: Astropy (Astropy Collaboration et al. 2022), MAAT (Ofek 2019), Matplotlib (Hunter 2007), NumPy (Harris et al. 2020), Photutils (Bradley et al. 2021), scikit-learn (Pedregosa et al. 2011), SciPy (Virtanen et al. 2020), SpiceyPy (Annex et al. 2020).
