Détails: Mis à jour : 13 mai 2022

After a lightcurve is obtained, it is necessary to convert it into the astrometric data of the asteroid, which will then be used to refine the asteroid’s orbit. Here we compare the results obained by the open source SORA astrometry module, to the approach based on the bayesian model sampler DIAMONDS.

We need to know not only how to properly determine the RA/Dec position, but also accurately determine their respective uncertainties, which have three main sources in single-chord events, two of which are highlighted in Fig. 1:

“Along-section” (AL or AT), tied to the event’s timing;
“Across-section” (AC), tied to how well the observed chord fits into the asteroid’s shape model (or, if there isn’t one, an elliptical or spherical model);
The uncertainties associated with the star’s astrometric data at the moment of the event, which depend on the star’s position and proper motion uncertainties.

Fig.1 - The AL axis is defined by the asteroid’s motion vector, and the AC is orthogonal to the AL. Different shape models can lead to different uncertainties for the same registered chord.

Many tools have been developed for this purpose, and the ACROSS team created two different routines, based on different approaches:

Least-Squares Fit, based on the SORA library;
Bayesian Inference, based on the DIAMONDS tool.

Here, we will briefly explain each method, and compare their performances (Fig. 3).

Bayesian Inference Method - BIM (DIAMONDS)

The Bayesian Inference Method is based on Bayes’s Theorem:

θ is the parameter vector, which is comprised of any number of parameters in occultation models, usually at least 4: the standard star+asteroid flux, the star’s flux (drop caused by the event), and two parameters for the timing, either ingress/egress or duration/centre epoch. We adopted the latter pair.

D is the set of data from the lightcurve.

M is the model applied to the data.

L is the likelihood function, which represents how the data is sampled.

π is the prior data on each parameter.

The model used for the occultation includes what is called a “super-gaussian”, where the gaussian’s peak is extended to instead be a plateau:

Here, F(t) is the lightcurve model, F0 is the original Flux, FA is the star flux, S is the super-gaussian, σ is its standard deviation, which controls the function’s decay, t0 is the centre epoch and L is the duration of the event. S takes the following form:

The priors can take any form, and we choose to use gaussian priors, centered on values that are estimated by scanning through the data or the lightcurve image.

The Bayesian Inference method is applied by using the tool DIAMONDS (Corsaro et al., 2014), and has been tested for single-chord events (Ferreira et al., 2020).

One of the results from DIAMONDS is a file with statistical data on the fit, from which we take the median and standard deviation of each of the four parameters.

Once we have the timing parameters, and their uncertainties, it becomes possible to determine the astrometric position and uncertainty of the asteroid, by grabbing the RA/Dec position and uncertainty of the star from the best available catalogue, which should be most of the time Gaia’s EDR3, known for being very precise, particularly when it comes to the proper motion (see Gaia documentation listed at the end).

We also need to know the asteroid’s ephemerides, specifically the RA/Dec velocity vector. This will determine the AL/AC orientation (AC is orthogonal to AL). It will also allow us to convert the timing uncertainty into a position uncertainty, expressed in arcseconds in MPC.

The AC uncertainty is more difficult to determine, especially in single-chord events. Since there are always at least 2 possible solutions to the chord, the uncertainty in this axis can be very large, regardless of the timing precision. An event with length close to the maximum expected should have a low uncertainty, while a grazing event should have an uncertainty almost as large as the diameter of the object! For more information on how this component is calculated, please read (Herald et al., 2020).

Finally, we can calculate the RA/Dec uncertainties of the asteroid’s position by making an error propagation with the contributions by the AC, AL and star components in each axis. The star’s uncertainties are already in RA/Dec format, while the AC/AL components have to be rotated based on the angle of the velocity vector.

Least-Squares Fit - LSF (SORA)

Stellar Occultation Reduction Analysis (SORA) is a Python-based, object-oriented library for optimal analysis of stellar occultation data. The package is divided into modules, each of which concerns an important aspect of an occultation (Body, Star, Observer etc.) and consists of functions to extract and process data for the purpose of fair analysis of the phenomenon (Gomes-Júnior et al., 2022). More information about the library can be found on their website and their GitHub repository (see bibliography). Using SORA’s functions and exploiting its capabilities along with the support from other libraries and databases, we were able to build a pipeline, through which we can analyze occultation data and extract the astrometric position of the asteroid at the time of the event, for multiple cases.

Analysis preparation:

This first stage of the pipeline consists of pre-analysis actions like data inputs and object instantiations. For each observation we have two files, one containing the Lightcurve - two columns first containing time instants and the second the corresponding flux at that instant - and one containing the rest of the observation details, like asteroid number, star name, site coordinates etc. . The format of the file containing the light curve is recognized by SORA and thus the file is used directly as input to instantiate the Lightcurve object.

The first of the participants of an occultation event is the asteroid. We instantiate the Body object using the asteroid’s number and specify ‘horizons’ for the ephemerides. SORA automatically grabs information from the Small Body Database (SBDB) and JPL Horizons (both links listed below). Continuing with the second participant, the target star, when instantiating the Star object, we use Astroquery (link at the end) to access the source id of GaiaEDR3 and ask SORA to retrieve data from this catalog. For the observer, we just enter the coordinates and the name(s) of the observer(s) as an EarthLocation object, provided by the Astropy library (The Astropy Collaboration et al 2018 AJ 156 123).

At this point we have the star, the observer, the body, and the reference time close to the instant that all three will line up, so we instantiate the Occultation object and SORA automatically recomputes the prediction of the event and provides mandatory, for the analysis, information.

Light curve analysis:

The Lightcurve object, which we mentioned earlier in its initialization, will contain a list of the light flux and one with the time instants, and since we have not done any processing yet, these are the raw data. We must ask SORA to detect the occultation using the occ_detect() function . This is done, in order to have a first look at the occultation, as well as to obtain a ‘mask’ of the event, in case it serves next.

The next step is to normalize the light flux, so that our data is compatible with the library that works with relative light flux and eliminates any data trends. This is done with the function normalize(), using a polynomial fit. There might be some exceptions that need special handling, where we might have to preselect some parameters (e.g. polynomial degree) with which the normalization is done.

Finally, a model can be fitted to the data to calculate the disappearance and reappearance instances of the star as accurately as possible. For this purpose, SORA provides the function occ_lcfit(), which executes a Monte Carlo chi square fit for occultation's lightcurve. The result of this procedure is a ChiSquare Object that contains the fitted model, which can be visualized along with the lightcurve, the immersion/emersion instances plus their uncertainties and the fitted parameters. Out of these, we are mostly interested in the timing of the event.

Astrometry reduction:

First, we need to add the instances of the observation to the Occultation object, SORA automatically “translates” timings into on-sky positions, thus creating the chord. The next step is to fit (generally) an ellipse to the available chord. For this, SORA provides the function fit_ellipse() and since we assume circular asteroids and we only have one chord, we keep oblateness and position angle constant and at zero and we let the center of the ellipse (f,g) vary. We also let the radius vary, but only inside a small interval defined by the radius uncertainty, which makes the analysis of cases where the chord is slightly larger than the reported diameter possible. SORA will then try to fit a number of ellipses with parameters inside the defined intervals and keep only the best fits, i.e. those with the smallest chi-square (Fig. 2). Before saving the results, we must determine the new astrometric position using the function new_astrometric_position() which does the calculations using the fitted parameters. We use this function to get the RA/Dec coordinates at the center epoch of the event with the geocenter as a reference, along with their uncertainties . Furthermore, we get the ephemeris offset both in distance and angular units.

Fig. 2 - The on-sky Observed chord and its uncertainties along with the fitted ellipses. The ellipse (black) is the best fit, i.e. lowest chi-square, and all the possible ellipses within 3-sigma (gray).

Fig. 3 - Uncertainties on the RA/Dec position of the asteroid at the centre epoch of the event for both methods. The distributions are similar.

Bibliography

Papers:

Corsaro, E. & De Ridder, J. 2014a, A&A, 571, A71;

Corsaro, E. & De Ridder, J. 2014b, DIAMONDS: high-DImensional And multi-MOdal NesteD Sampling, Astrophysics Source Code Library;

Ferreira, J. F., Tanga, P., Machado, P., & Corsaro, E. 2020, Astronomy and Astrophysics, 641, A81;

Gomes-Júnior, A., Morgado, B., Benedetti-Rossi et al. 2022. SORA: Stellar occultation reduction and analysis. Monthly Notices of the Royal Astronomical Society, 511(1), pp.1167-1181.

Herald, D., Gault, D., Anderson, R., et al. 2020, Monthly Notices of the Royal Astronomical Society, 499, 4570–4590;

Price-Whelan, A., Sipőcz, B., Günther et al. 2018. The Astropy Project: Building an Open-science Project and Status of the v2.0 Core Package. The Astronomical Journal, 156(3), p.123.

Robitaille, T., Tollerud, E., Greenfield, P., et al., 2013. Astropy: A community Python package for astronomy. Astronomy & Astrophysics, 558, p.A33.