Kuiper Belt Occultation Predictions

Software from the MegaPipe data pipeline (Gwyn 2008) was used to produce the point-source catalog and to astrometrically calibrate each of the MegaPrime images. Starting with images already preprocessed by the Elixir pipeline (Magnier & Cuillandre 2004a), for each image we produced a source catalog with positions in pixel coordinates using SExtractor (Bertin & Arnouts 1996). SExtractor's parameters were set such that only fairly bright sources were detected; the detection criteria are set to flux levels 5 sigma above the sky noise in at least five contiguous pixels. The catalogs were further cleaned of cosmic rays and extended sources, leaving only point sources. Source centroids were found using SExtractor's simple centroid method. More complicated methods such as Gaussian or PSF fitting were found to provide no noticeable benefit.

Each of the pixel coordinate catalogs were matched to an external astrometric catalog, either the SDSS or 2MASS (Skrutskie et al. 2006). The SDSS DR9 (Ahn et al. 2012) provides a superior source density, and very small astrometric errors (0.07–0.09''). Where it was not available, 2MASS was chosen as the reference catalog (astrometric error 0.08–0.09''). We did not use UCAC (Zacharias et al. 2004) because it is fairly shallow; most of the sources are saturated in MegaPrime images, leaving too few sources for accurate astrometric calibration. The USNO catalogs were also considered: they go deeper and have a higher source density. The astrometric errors (0.4–0.6'') on each source are larger than with the 2MASS catalog. Empirically, it was found that using the 2MASS catalog gave the smallest astrometric residuals.

Initially, the pixel coordinate catalogs were matched to the external astrometric catalog on a chip-by-chip basis. The initial Elixir (Magnier & Cuillandre 2004b) preprocessing done by CFHT provides an initial astrometric calibration, which is typically accurate to better than 1''. The matching is therefore relatively simple. Sophisticated techniques (such as the quad-matching method used by Lang et al. 2010) are not required. During the initial matching process, any catalog source that was more than 1'' away from the nearest observed source was ignored. Once the catalogs were matched, the transformation was computed. For the initial match, we used a second-order polynomial in x and y. This transformation was used to refine the matching of sources in the images to sources in the external catalog, and the transformation was recomputed. The second transformation was computed slightly differently. The MegaPrime distortion map can be adequately described by a polynomial with second and fourth order terms in measured radius, r, measured from the center of the mosaic. This is given by

where R is the true radius, and a₁ and a₂ are distortion coefficients. The coefficients of this polynomial were determined for all 36 chips simultaneously. In addition, a linear distortion map was computed for each chip. The combination of a global, nonlinear transformation and 36 local, linear transformations sufficiently describes the distortion, such that additional complexity in the map does not detectably improve the solution. The global transformation takes care of most of the distortion caused by the MegaPrime optics, while the linear transformation takes care of any non-coplanarity of the detector array as well as the effects of differential refraction. We avoid the usually adopted method of determining the full distortion map which uses a third-order polynomial transformation for each chip. This utilizes up to 20 parameters per chip, or a total of 720 parameters for all 36 chips. This number of parameters is uncomfortably close to the number of sources available for astrometric calibration. In most MegaPrime fields of view one typically finds 2000 suitable sources, but the number can be as low as 1000.

Once transformations had been determined, they were used to convert the pixel coordinates in the individual source catalogs to R.A. and decl. The catalogs from each image were merged to produce a master catalog. For each source from each catalog, all the catalogs are checked for matching sources. A match occurs if a source in one catalog lies within 2'' of a source in another catalog. The matching sources must also have measured magnitudes within 1 mag of each other. Once all the matching sources had been found, their measured positions from the different catalogs were averaged. To avoid confusion in the matching, if two sources in the same catalog lie within 4'' of each other, both sources are discarded.

Next, the astrometric calibration of each input image was repeated using the merged master catalog as the astrometric reference. Using these new calibrations, the source catalogs for each image are converted from pixel coordinates to R.A. and decl. These catalogs are merged as described previously to produce a new merged master catalog. This second catalog is used to calibrate the images a third time. We refer to this step as "merge-by-catalog".

The images were photometrically calibrated by one of three methods:

• 1.  
  If the images overlapped the SDSS, it was used for calibration. To account for slight differences in photometric systems due to detector and filter differences, the SDSS photometry was transformed to the MegaPrime system as follows:

  

  This results in an absolute photometric calibration accurate to about 0.01 mag.
• 2.  
  If the image lay outside the SDSS, the Elixir photometric calibration was used for data taken on photometric nights. The absolute photometric calibration in this case is slightly worse, typically 0.03 mag.
• 3.  
  Images taken on nonphotometric nights were calibrated using parts of the image which overlapped with images photometrically calibrated with one of the previous methods.

The photometrically calibrated catalogs were merged using a simliar method to the astrometric catalog merging described above to produce a merged master photometric catalog. The final photometric calibration was done using this catalog. Merging the photometric catalogs and recalibrating in this way does not significantly improve the external photometric calibration, but ensures that the internal image-to-image zero-point calibration is typically better than 0.005 mag.

The astrometrically and photometrically calibrated images were then combined using SWarp (Bertin 2004). SWarp is a program that resamples and stacks multiple images onto a projection defined by the image world coordinate systems. The background was computed on a 128 pixel grid using a median filter and removed. The astrometric distortion was removed and the photometric scaling was applied. The images were resampled to a common grid using a Lanczos 3 pixel kernel. The scaled, resampled pixels are combined using a median. We refer to this step as "merge-by-pixel". SExtractor was then run on the resulting image to produce the final catalog of point sources which we refer to as the Master Point Source Catalog (MPSC). Finally, the individual images were re-calibrated using the MPSC as an astrometric reference.

After each iteration of the astrometric calibrations, the image-to-image astrometric residuals were checked. After the initial match to the external catalog, the astrometric residuals range from 0.08 to 0.1'' rms, typically slightly higher if the 2MASS catalog was used as an external reference, slightly lower if the SDSS was used. After using the second master catalog, the one generated using the "merge-by-catalog" method, the astrometric residuals are ∼0.06''. The residuals get marginally lower if the merge–recalibrate–merge cycle is repeated. After using the final "merge-by-pixel" master catalog the astrometric residuals are typically 0.03–0.04'' between two individual images. The residuals between the individual images and the corresponding master catalog are typically 60–70% smaller than residuals between two individual images.

Figure 2 shows the astrometric residuals between two input images after matching to the "merge-by-pixel" master catalog. The top left plot shows the astrometric residuals as a vector field. The lengths have been greatly exaggerated. The bottom left plot shows the residuals in R.A. and decl. as a scatter plot. Histograms of the residuals in both directions are also plotted. The two plots on the right show the residuals in R.A. as a function of decl. and the residuals in decl. as a function of R.A. The residuals are seen to be on the order of 0.02'', while there are a few misidentifications (indicated by the outliers in the scatter plots and by longer than usual lines in the vector field), no large systematic shifts are apparent.

Zoom In Zoom Out Reset image size

When the KBO of interest fell in one of the CFHT images, its position was measured with respect to the MPSC. For astrometry of the KBOs measured in the CFHT data, we adopt an uncertainty of 0.04'', typical of the astrometric residuals found after the "merge-by-pixel" step, the last step in producing the MPSC.

Along with the target astrometry provided by the MegaPrime data, additional astrometry of the targets was measured from images taken with the GMOS detector (Gemini program GN-2011B-Q-60). Each target was visited multiple times. During each visit, a pair of images was taken with small ∼30'' dither between pairs. All exposures were taken in r^'-filter and exposure times were tuned such that the resultant photometric signal-to-noise ratio (S/N) was ∼40–50. This ensured that the resultant astrometric error was dominated by the match between the images and the MPSC and not by the quality of the object measurement. The target's observed positions are presented in Table 2.

After standard image reductions (bias removal and flat fielding) the astrometric plate solutions for the GMOS data were found by matching those images to the MPSC. The match was done with the use of SCAMP (Bertin 2006). SCAMP is a program that produces astrometric solutions which match one point source catalog to a reference catalog. During the matching process, each image was treated individually, providing as independent individual astrometric measurements as possible. The plate solution of GMOS was found to be adequately described by a simple CD matrix solution; utilizing higher-order terms revealed no improvement in the astrometric solution. Experience has shown that SCAMP typically produces unreliable astrometric uncertainties. That is, the root mean square (rms) astrometric scatter of a solution is typically erroneous in one of two ways. Either the reported rms is small despite the fact that the astrometric solution is clearly erroneous, or the rms value is found to be significantly smaller than the scatter in repeat images of the same targets. Repeated measurements of a KBO at a given epoch have demonstrated that the plate solutions provided by SCAMP result in a ∼0.08'' scatter in the astrometry of the tracked KBOs. Therefore, we forgo use of the quoted SCAMP RMS values, and adopt an astrometric uncertainty of 0.08'' for all GMOS measurements.

The second step in occultation predictions is the generation of accurate ephemerides. As discussed in § 1, a common approach is what we refer to as the constant offset method, in which the nominal ephemerides of the targets in question are offset by a value representative of the typical differences between the observed and predicted positions of the object.

One merit of the constant offset method comes as a result of the use of a nominal ephemeris for the target object. By use of a previously determined ephemeris, all past reported astrometry and an orbit determination from those data are automatically considered in the predictions. While it is certainly true that past astrometry are typically not of the quality required for occultation predictions, many measurements spanning a multiyear baseline provide modestly accurate ephemerides that, with small astrometric corrections, can be suitable for occultation predictions.

The accuracy of predictions made by the constant offset method primarily suffer from one major issue. In reality the difference between the nominal-offset ephemeris and the true ephemeris is not actually constant. The difference between the two ephemerides depends primarily on the quality of the nominal ephemeris. A small error in the nominal ephemeris can result in both a periodic shift (with period of a year) and a nearly linear shift between the nominal-offset ephemeris, and the true ephemeris (see Fig. 1). Even for the best ephemerides, the amplitude of the periodic offset can be of order ∼0.02''. As a result, without frequent tracking and updates to the offset, predictions with this method can be unreliable.

A different technique involves adjustment of the orbital elements of the nominal ephemeris to match available astrometry. That is, with appropriate tweaks to the nominal orbital elements, the ephemeris itself can be corrected and used directly for occultation predictions. We consider the orbital element adjustment approach and its comparison with the constant offset method here.

We start with the orbital elements provided by the astDys catalog and make use of Orbfit¹⁶ in calculating ephemerides. To determine the appropriate orbital element corrections, we adopt a maximum-likelihood approach. Specifically, given adjustments to the semimajor axis, eccentricity, inclination, longitude of ascending node, argument of perihelion, and mean anomaly δa, δe, δI, δΩ, δω, δM, we adopt the log-likelihood

where α_i and β_i are the ith observed R.A. and decl. of the object, and and are the R.A. and decl. predicted from the adjusted ephemeris of the object in question. The terms σ_i,α and σ_i,β are the astrometric uncertainty of the ith measurement (recall that we adopt 0.04'' and 0.08'' uncertainty in the MegaPrime and GMOS observations, respectively).

In equation (3), P(δa,δe,δI,δΩ,δω,δM) is the log-prior on the orbital elements. For this, we turn to the element errors in the nominal ephemeris. We however, choose not to apply any priors on the orbital element angles. This decision is motivated by the consideration of the coordinate origin of the astrometric catalog (2MASS or SDSS), which we first consider when generating the MPSC. Each catalog will have a slightly different coordinate origin. The ephemerides provided by astDys reference a different zero point (usually that defined by the USNO A catalog). As a result of the differing origins, the orbital angles Ω, ω, and M, are incorrect at the level of accuracy required. Similarly, as the nodal angle has changed, so has the orbital inclination. It is only the semimajor axis and eccentricity that are not affected by adjustment of the reference. We adopt Gaussian priors on a and e with standard deviations equal to the uncertainties in those elements provided with the astDys ephemerides.

When determining the orbital element corrections with equation (3), we only consider the MegaPrime and GMOS observations we gathered ourselves and matched to the MPSC. These observations only span 1–2 years, and as a result partially avoid potential zonal errors which affect the astDys orbital elements. The disadvantage of this approach is that past observations are not directly used by our method. Rather, those data are only used as a prior, a choice which reflects how much prior information we can safely extract from the nominal astDys elements without producing unreliable results.

To converge on a maximum likelihood, we utilize MCMC Hammer (EMCEE, Foreman-Mackey et al. 2012). MC Hammer is an affine invariant Markov Chain Monte-Carlo sampler that rapidly converges on the maximum likelihood solution in a minimum number of steps. For the routine, we adopt 100 walkers, and utilize a 100 step burn-in phase. During the maximization, EMCEE determines the autocorrelation time, τ—the number of steps required such that further samples are independent and distributed as the posterior likelihood function. Thus, the routine must be run for multiple τ to ensure that the likelihood space is accurately evaluated. For all orbital elements, τ was found to be ∼5–8. Thus, our routine was run for 150 steps after the burn-in phase, or roughly 20 or more autocorrelation times. Additional steps were found to not improve the resultant likelihood evaluations.

As an endproduct of this routine, we can extract confidence intervals on each of the orbital elements. From these confidence intervals, we can determine the astrometric uncertainty on the target for any future date and time. No direct method of positional astrometric uncertainty can be generated from the constant offset method.

To compare the performance of the constant offset and element adjustment methods, we set up a theoretical situation. For this example, we consider a theoretical Eris, and assigned as its elements the elements reported by astDys for the real Eris. Using Orbfit, we then simulated the historical observations of the theoretical KBO on the dates of all reported astrometry for the real Eris—the dates were extracted from the Minor Planet Center.¹⁷ We note that Orbfit includes planetary perturbations, a necessary consideration for the decades-long arcs of the KBOs we consider. Noise was added to the simulated observations to account for typical random uncertainties in reported astrometry. A Gaussian distribution with standard deviation 0.3'' in both R.A. and decl. was used. We also attempted to include zonal errors in the simulated observations. This was done by scattering the observations by a small value every time the theoretical KBO moved more than 3°; we adopted a Gaussian distribution with width 0.08'' in both R.A. and decl. In a similar vein, we generated fake tracking observations (like those we present from MegaPrime and GMOS) adopting 0.04'' astrometric uncertainties.

For the theoretical KBO, a set of orbital parameters was determined from the simulated historical observations with Orbfit (ignoring the simulated tracking observations), the same routine used for the astDys elements. We then utilized our maximum likelihood approach to determine (δa, δe, δI, δΩ, δω, δM), the six element corrections which provided the best match to the simulated tracking observations. We also evaluated the best offset which minimized the residuals between the simulated tracking observations and the nominal ephemeris.

The results of our simulation are shown in Figure 1, where we compare the results of both ephemeris correction processes. This figure clearly demonstrates one of the difficulties with occultation predictions. Both methods suffer from a 365 day period oscillation in the difference between the predicted and actual ephemerides. In addition, the constant offset method suffers from an increasing error in decl. with time.

One advantage of the element correction approach is clear. The approach can produce more accurate ephemerides further away from the last epoch of observation than the constant offset method. The quality of our element correction approach however, depends critically on the quality of the astDys ephemeris. In addition, without sufficient baseline in the astrometry used in equation (3), the resultant element offsets can be incorrect; simulations suggest that at the very least, 1 year baselines are required. The element correction approach is also sensitive to discrepant tracking astrometry which can result in ephemerides that quickly deviate away from the true ephemeris. Because the constant offset method averages all tracking astrometry, this method is much less sensitive to discrepancies in the tracking observations.

Which of the constant offset method or the element adjustment approach produces the most reliable occultation predictions will be different for each object, and the determination of which essentially requires future observations. Eventually, the element adjustment approach should surpass the constant offset method. How long is required is not easily determined a priori, but simulations suggest a ∼2 year baseline or longer will allow the element adjustment approach to surpass the constant offset method. Both methods however, suffer from different effects, and as a result, comparison of occultations predictions produced with both methods should be considered until the results of the element adjustment approach have been demonstrated to be superior on an object by object basis. Thus, we choose to present the results of both methods when reporting candidate occultations.

It should be noted that, along with the ephemeris uncertainty, there is also uncertainty in the stellar astrometry. From the data used in generating the MPSC, each star is observed multiple times, with a positional accuracy of ∼0.04'' in both R.A. and decl. The quoted stellar position represents a global mean of each measurement. As we will discuss further below, it appears that the uncertainty on the mean position does not decrease as the square-root of the number of times the star is observed, but rather, appears roughly a factor of 2 worse than that expectation.

During the generation of the MPSC, we make no effort to account for the effects or proper motion, or chromatic differential refraction. The first effect, that of proper motions, was ignored during our analysis because proper motions were not determined for either of the reference catalogs we used. As a result, the stars used in the first step of the astrometric calibration will have moved slightly since they were observed by 2MASS (1997–2001) or the SDSS (2000 to present). Estimates of the proper motion of sources in the 2MASS catalog are available from the PPMXL catalog (Roeser et al. 2010). Experimentation with that catalog revealed that the proper motions measured by PPMXL were only accurate enough to improve the positions of the brightest sources in the 2MASS catalog, those which were primarily saturated in our observations. No noticeable improvement in positional accuracy of fainter sources from which our astrometric solutions were derived, was found. As a result, we chose to adopt the cataloged positions for the first stage of our astrometric calibrations and make no effort to correct for proper motions.

Tholen et al. (2013) point out that not considering proper motions can result in zonal errors. In the small fields we observe to produce the MPSC, these zonal errors would manifest themselves as a nearly constant systematic offset between the frame of the refence catalog (2MASS or SDSS) and the frame of the MPSC. The offset will be nearly constant across the small fields covered by the MPSC for each object.

To determine if ignoring proper motions could result in inaccurate predictions, we performed a test of the predictions for the object 2002 TC302, in which the positions of stars in the reference catalog (2MASS) were scattered to mimic the effects of proper motions. The distribution of proper motions was extracted from the PPMXL catalog in a 2°-wide patch around the position of the KBO. The R.A. and decl. positions of each source were then randomly scattered according to the distribution of R.A. and decl. proper motions. Analysis of the PPMXL proper motions suggest that ∼15% of sources have moved more than 1'' since the 2MASS catalog was created, and hence would be ignored during the MPSC generation. After the sources were randomly scattered, this number increased by a factor of . The full occultation predictions routine (MPSC generation, determination of KBO astrometry, and ephemeris correction) was then applied using the randomized reference catalog. As predicted, compared to the astrometry referenced to the non-scattered 2MASS catalog, a systematic offset in both the MPSC source positions and that of the KBO were found, with amplitude of nearly 0.1'' R.A., and 0.05'' decl. In addition, the astrometry from the scattered catalog had a slope of ∼0.05'' per degree of R.A. in both the R.A. and decl. axes. The effects of the offset and slope produced a shift of ∼300 km in the predicted shadow tracks of 2002 TC302. Even the most accurate predictions have ∼1000 km uncertainties in the shadow path, which is entirely a result of ephemeris uncertainty. As ephemeris uncertainty grows rapidly with time away from the last tracking observation, the prediction errors induced by not accounting for proper motions will always be small compared to the ephemeris error.

The reason why not accounting for proper motions has such little effect on the predictions can be easily understood. Recall that during the MPSC generation, an external catalog is only used in the first stage of generating astrometric solutions. An internal catalog is used in later stages to refine those solutions and generate a self-consistent MPSC. As a result of this internal reference to the observed sources, small variations in the solution caused by proper motions are reduced at each subsequent iteration of our calibration routine, mitigating these issues. In general, as shown by our test with 2002 TC302, this systematic offset does not significantly affect the quality of the astrometric solution or the KBO ephemeris as the astrometry of both the MPSC and the KBO are equally affected by the offsets.

While proper motions in general do not affect the quality of the astrometric solution, they will cause a degradation of some individual predictions with time away from the observations used to generate the MPSC. Take for example, a star with proper motion of 0.1'' year^-1—roughly 15% of all stars have at least this proper motion. At the time of observation of that star for MPSC generation, its position is known to the precision of the MPSC. At typical KBO distances, the high proper motion of the star will have moved the shadow track more than 1000 km from that predicted in just 6 months. Clearly, advanced monitoring of the target stars is warranted for particularly profitable events.

Like proper motions, we make no effort to account for the effects of chromatic differential refraction (CDR), primarily due to the single band observations we acquired for MPSC generation. Our observations were taken over a broad range of airmasses, but less than 2 in all cases. In the standard conditions on Mauna Kea (Cohen & Cromer 1988), the reddest stars, with (V - I) ∼ 3, will experience a shift of amplitude ∼0.05'' at an airmass of 2 compared to its relative position when observed at zenith (Stone 1984; Monet et al. 1992). For most stars, (V - I) ≈ 0.8 resulting in a shift is of order 0.01''. The reddest KBOs eg. Quaoar and Sedna with (V - I) ∼ 1.3, will experience a slightly larger shift. The primary consequence of CDR is a fundamental limit on the precision of the MPSC; unless observations of the MPSC sources are taken at random airmasses, their astrometry cannot be more precise than the shift cause by CDR at their average observed airmass. For most stars, and our observations, we conservatively estimate the amplitude of this effect at up to a ∼0.02'' shift of individual stars compared to the astrometry of the KBO. Improvements beyond this limit cannot be made, unless either the effects of CDR are accounted for, requiring multiband observations, or observations that are restricted to a small range in airmasses.