Pipeline for the Detection of Serendipitous Stellar Occultations by Kuiper Belt Objects with the Colibri Fast-photometry Array^*

Collisional evolution models have been used to estimate the size-frequency distribution of the Kuiper Belt (Davis & Farinella 1997; Kenyon & Bromley 2004; Campo Bagatin & Benavidez 2012; Schlichting et al. 2013). Both models and observation suggest that the size distribution of objects in the Kuiper Belt follows a broken power law, with collisions modifying the size distribution between the power laws for large and small Kuiper Belt objects (KBOs; Kenyon & Bromley 2004; Pan & Sari 2005). While over a thousand KBOs larger than 15 km have been discovered, objects smaller than 10 km are too dim to be observed by direct detection, causing the distribution of small objects to be poorly constrained.

With an accurate model and a well-constrained estimate of q (the power-law index for the distribution of small KBOs) one could obtain insight into the material strength of these objects (Pan & Sari 2005), as well as constrain the initial planetesimal population of the Kuiper Belt (Schlichting et al. 2013) and determine whether the kilometer-sized KBO population is enough to supply the Jupiter-family comets (Dones et al. 2015). However, without knowledge of the number of small KBOs, q cannot be well constrained.

While small KBOs cannot be imaged directly, their presence can be inferred from serendipitous stellar occultations (Bailey 1976). For an occulting KBO smaller than the Fresnel scaledefined as $\sqrt{\lambda D/2}$ where λ is the observation wavelength and D is the distance between the occulter and observerand for an occulted star with a projected diameter that is also small relative to this scale, an occultation event will be dominated by diffraction effects (Section 4).⁶ With a mean observing wavelength of 550 nm and a typical KBO distance of 40 au, this characteristic Fresnel scale is 1.2 km. Therefore, when a kilometer-sized KBO passes in front of an angularly small star, a characteristic diffraction pattern can be observed in the star's light curve (Roques et al. 1987).

These Fresnel diffraction patterns are observable for only a fraction of a second, but with sufficiently fast photometry, an occultation survey is able to probe the Kuiper Belt for kilometer-sized objects. In this paper, we discuss a preliminary system for such an occultation survey.

A handful of occultation surveys have been performed in recent years in an attempt to characterize the size distribution of small objects in the Kuiper Belt. Early attempts by Chang et al. (2006) and Roques et al. (2006) each claimed to have found positive detections, but the rates found in these studies were incompatible with each other and with theoretical models (Schlichting et al. 2009), and the calculated distances to the Roques objects would not place them within the Kuiper Belt (Zhang et al. 2013). Subsequent studies have suggested the Chang detections were the result of systematic effects, not KBO occultations (Jones et al. 2008; Blocker et al. 2009).

Ground-based surveys have not yet produced any viable occultation candidates, due in part to the difficulty of performing quality fast photometry from the ground. One such survey, a dedicated array called TAOS, did not detect any KBOs in 292,514 star hours spread over seven years of operation (Zhang et al. 2013). With its sampling rate of 5 Hz, the duration of a single frame is similar to the duration of a diffraction event; the TAOS null detection is perhaps unsurprising. A future version of the project, TAOS II, will use a sampling rate of 20 Hz to increase the likelihood of a successful detection (Lehner et al. 2014).

Bickerton et al. (2009) found that the Nyquist sampling rate for detection of a KBO at 40 au in the optical is 40 Hz. Recent ground-based efforts have focused on these high sampling rates, such as occultation surveys with Megacam on MMT (Bianco et al. 2009) and IMACS on Magellan (Payne et al. 2015), which sampled at rates of 30 Hz and 40 Hz, respectively. Neither have reported any detections in their 220 star hour and 10,000 star hour campaigns, although they were able to provide constraints on the surface densities of KBOs.

A potential detection of a KBO occultation was obtained by Schlichting et al. (2009), who analyzed 12,000 star hours of Hubble Space Telescope archival data taken at a cadence of 40 Hz. A follow-up survey of 19,500 additional star hours yielded a second occultation candidate (Schlichting et al. 2012).

The most stringent constraints on the small KBO population have been placed by TAOS. TAOS established an upper limit of 3.34 to 3.82 on the index q of the small-KBO power law $\left(\tfrac{{dn}}{{dR}}\propto {R}^{-q}\right)$, with the size distribution of small KBOs normalized using surface densities of 38.0 deg⁻² (Fuentes et al. 2009) and 5.4 deg⁻² (Fraser & Kavelaars 2008) at the break in the size-frequency distribution (R = 45 km), respectively (Zhang et al. 2013). Furthermore, some surveys have produced model-independent constraints at specific KBO diameters,⁷ such as Bickerton et al. (2008); Jones et al. (2008); Bianco et al. (2009) (see Bianco et al. 2009, Figure 12) and Schlichting et al. (2009).

Recent direct observations of impact craters on Pluto and Charon by the New Horizons probe have also placed constraints on the size-frequency distribution of the Kuiper Belt. Even after consideration of geological factors, these observations are suggestive of substantially lower densities of small KBOs than the maximum permitted by the upper limits established by previous occultation surveys (Singer et al. 2016). This low cratering rate may therefore suggest a potential unreliability in the Schlichting et al. detections (Greenstreet et al. 2015).

Because of the rarity and brevity of sub-km KBO occultations, a successful KBO detection campaign must optimize both sampling frequency and star hours. The expected detection rate for KBO occultations is

$$\begin{eqnarray}&&\mu =2{bv}{\rm{\Sigma }}{\left(\displaystyle \frac{180}{\pi D}\right)}^{2},\end{eqnarray} \tag{ 1 }$$

where b[m] is the maximum detectable impact parameter, v[ms⁻¹] is the relative perpendicular velocity of the KBO, Σ[deg⁻²] is the sky surface density, and D[m] is the distance between observer and occulter (see Bickerton et al. 2009, Section 3.1 for further description of these parameters). Considering a minimum detectable KBO diameter of 500 m, the time between occultations of the same star is on the scale of years to centuries, depending on the estimate used for the size-frequency distribution.

We are developing a robotic telescope array, Colibri, that will be dedicated to KBO occultation monitoring. Colibri will consist of three 50 cm telescopes on the vertices of a triangle with sides ranging from 130 to 270 m. Each telescope system will image at $f/2.2$ onto a 1024 × 1024 electron-multiplied charge-coupled device (EMCCD) detector. The field of view of each camera will be 069 on a side. While modern 1024 × 1024 EMCCDs offer sampling rates that are 25 Hz at best over the full detector, by binning pixels 1 × 2 we will attain sampling rates of 40 Hz. The limiting magnitude of individual images is expected to be 13.5 in the broad band optical at signal-to-noise ratio (S/N) = 10 (similar to the G band pass on Gaia).

The number of star hours obtainable from Colibri is effectively unlimited, and the 40 Hz cadence will sample ∼200 ms-long occultations over a number of frames. As Colibri consists of three telescopes, candidate events can be corroborated by each independent camera and therefore a lower threshold S/N can be used to flag candidate events in each data stream, as compared to a single-camera system.

For the preliminary trials discussed in this paper, the Colibri pipeline was operated at 17 Hz using two independent EMCCD cameras, each mounted on a 75 mm $f/1.3$ video lens (D = 5.8 cm). However, when the main campaign begins in late 2017, the array will be operating with three independent 50 cm telescope systems at 40 Hz.

In Section 2, we discuss the specifications of Colibri's preliminary observation campaign. Our data pipeline and detection algorithm are described in Section 3. In Section 4 we detail the theoretical expectations for the occultation light curves, and build a set of time-domain kernels for use in matched-filter detections of occultations. Sections 5 and 6 discuss false positive rates and simulated event retrieval, respectively. Results and future work are outlined in Section 7.

Pending the completion of the Colibri array, we performed the preliminary trials of the Colibri data pipeline (Section 3) using the EMCCD Canadian Automated Meteor Observatory (EMCCD CAMO) array (Figure 1; Table 1) at Elginfield Observatory in Ontario, Canada (43°11'33''N, 81°18'57''W). As discussed in Section 1.3, EMCCD CAMO operates at half the cadence and with fewer, smaller telescope systems than the full Colibri array. The purpose of these preliminary trials, therefore, is to verify the efficacy of the detection software and the plausibility of a multi-EMCCD occultation detection system at the Elginfield site, rather than to detect KBOs.

Table 1.  EMCCD CAMO Specifications

Detectors	2
Resolution	1024 × 1024
Frame rate	16.7 fps
Exposure time	59.91 ms
Bit depth	16 bit
Pixel scale	359/px
Limiting magnitude	8 mag at S/N = 10, unfiltered V band
Camera	Nüvü HNü 1024 EMCCD
EM gain	100x (set in Nüvü software)
Optics	Navitar 75 mm f/1.3
FOV size	102 × 102

Download table as:  ASCIITypeset image

EMCCDs' onboard gain reduces the effects of read noise, resulting in higher S/Ns than typical CCDs and making them ideal for fast photometry (Giltinan et al. 2011). These devices have recently begun to be employed in occultation surveys, such as in the Bickerton et al. (2008) study, the LCOGT network (Bianco et al. 2013) and the CHIMERA photometer on the Hale telescope (Harding et al. 2016).

The EMCCD CAMO array is autonomous, with data collection automatically beginning when skies are clear and pausing when they are cloudy. Cloud cover is monitored by an independent camera directed at Polaris.

The Polaris monitoring module uses a Hi Cam HB310E analog NTSC camera digitized to 8 bit which has a 2° FOV. Every 5 minutes, the camera does a 300 frame stack and the brightest pixels in the stacked image are selected. Standard stellar photometry routines are used to find an apparent magnitude for Polaris. As long as the apparent magnitude is within 0.5 mag of the brightest magnitude that Polaris reaches, the skies are considered clear. As the EMCCDs are pointed within 15° of Polaris, the sky conditions at Polaris are a very good proxy for the conditions in the EMCCD FOV.

The automation and overall design of the CAMO system has been described in detail elsewhere (Weryk et al. 2013).

A field with high densities of both angularly small stars and KBOs is desirable for occultation surveys, although accurate photometry is difficult in overly dense stellar fields. Stellar field density is maximized in the galactic plane and KBO density is maximized in the plane of the Kuiper Belt, which is within a few degrees of the ecliptic (Brown & Pan 2004).

The EMCCD CAMO array's principal function is meteor detection and therefore it does not track siderally, but remains fixed in Az/Alt pointing at a zenith angle of 2489 and an azimuthal angle of 2286. For our observations (described in Section 2.3), this corresponds roughly to ecliptic latitudes of 45°–82° and galactic latitudes of 37°–58°. While these telescope pointings are not optimal for KBO detection, selection of an ideal field is in progress and will be implemented for the main Colibri campaign. This selection process involves a full characterization of the stars in the stellar field (distances, spectral types, projected angular diameters, etc.), which will allow the parameters of any candidate occultation to be thoroughly checked for consistency (see Section 6.2). Such a characterization process is unnecessary for the test field, as no occultation events are expected in high ecliptic latitude data. Instead, the high ecliptic latitude data taken in these preliminary trials will serve as a control sample to the main campaign, as discussed in Section 7.

The overlap between the two EMCCD CAMO fields is ∼75%. As the distribution of stars that pass the detection threshold is relatively uniform across the detector, two simultaneous time series are able to be obtained for roughly 75% of the stars detected by a single camera.

We used the Colibri pipeline to examine a total of 11,900 star hours. This includes 1200 star hours observed with a single EMCCD camera on 2016 August 1 and 10,700 star hours of observations with two EMCCD cameras between 2017 February 19 and April 1 (Table 2). Because of the incomplete (∼75%) field overlap between the two cameras, the dual-camera observations represent a total of 4000 hours on stars with two cameras observing simultaneously.

The resulting data sets comprise 14,456 5-minute light curves from the single-camera observations on 2016 August 1, and ∼48,000 5-minute light curve pairs from the dual-camera observations in early 2017. The extracted time series in Table 2 represent the number of 5-minute point-source time series recorded each night, the accepted time series represent the fraction of these that pass the acceptance criteria outlined in Section 3.1, and the dual-telescope time series represent the number of simultaneous light curve pairs obtained by the two EMCCD cameras.

The single-camera light curves are used in Sections 5 and 6.1 to determine false positive detection rates and to simulate the retrieval of occultation parameters.⁸ The dual-camera light curve pairs represent the actual KBO occultation experiment that we performed with the EMCCD CAMO setup, and in which we find one candidate event, likely a false positive (Section 6.3).

The pipeline preprocesses each 5-minute data cube to remove dark and flat fielding effects, with the dark frame created from a median combination of 167 dark exposures and the flat field created from a median combination of 80 frames taken at 5-minute intervals throughout one observing night.

We then extract stellar point sources using Source Extractor for Python (Bertin & Arnouts 1996; Barbary 2016). The pipeline also uses the Astrometry.net API for plate solving (Lang et al. 2010) and the Astropy library within Python for various file management and convolution functions (Robitaille et al. 2013).

The pipeline only considers time series where the net counts per frame in an r = 5 pix (180'') aperture are greater than 5000 after preprocessing and background subtraction (corresponding to a S/N ≈ 9 from empirical measurements), where the star does not leave the field of view within the first minute, and where the tracking algorithm does not otherwise fail. For the preliminary trials, this represented approximately 1/4 of time series extracted from EMCCD CAMO data.

There are multiple techniques for numerically searching time series for KBO occultations, such as the variability index method used by Roques et al. (2006), the rank-probability method used by TAOS (Zhang et al. 2008), and the template cross correlation method used by Bickerton et al. (2008). We designed a novel algorithm for the Colibri array, although this structure shares some similarities with the Bickerton et al. method.

Zoom In Zoom Out Reset image size

Instead of performing template matching across the entire 5-minute time series (as with the Bickerton et al. algorithm), we first convolve each light curve with a Mexican hat filter, the width of which is determined by the characteristic duration of a KBO occultation at the camera sampling rate. The duration of an occultation is ∼200 ms for a KBO that is small relative to the Fresnel scale and viewed at opposition⁹ (see Equations (3)–(5)), which corresponds to a Mexican hat filter with a width of three frames in the 17 Hz-sampled data. A KBO much larger than the Fresnel scale may generate a longer-duration occultation event, as the width of these occultations then becomes governed by the geometric regime rather than Fresnel optics (see Figure 2 and Nihei et al. (2007), Figure 3). While our Mexican hat filter is not optimized for these larger KBOs, such events are more easily detected due to the near-complete extinction of the background star (see Nihei et al. (2007), Figure 5). Indeed, our three-frame-width filter remains sensitive to objects as large as at least 5 km, as evidenced by the retrieval rates that we will discuss in Section 6.1.

Zoom In Zoom Out Reset image size

The full algorithm for the Colibri pipeline consists of four phases:

• 1.  
  The first step comprises the three-frame Mexican-hat convolution and concludes with the determination of the most significant dimming event for each 5-minute convolved time series.
• 2.  
  In the second step, flagged events with dips deeper than 3.75σ below the mean of the convolved time series are allowed to proceed to the next step in the pipeline; this 3.75σ value was empirically selected to optimize event retrieval while retaining a sufficiently low false positive rate.
• 3.  
  In the third step, the candidate event is compared to a set of kernels modeling the diffraction patterns for various parameter sets (see Section 4). With s = time series, k = kernel, P = Poisson noise on s, and σ = standard deviation of the time series¹⁰ , the noiseless kernel is a successful match to the noisy data when the following condition is true:

  $$\begin{eqnarray}&&\displaystyle \sum _{i=1}^{{len}(k)}\sigma {/P}_{i}\gt \displaystyle \sum _{i=1}^{{len}(k)}| {s}_{i}-{k}_{i}| {/P}_{i}.\end{eqnarray} \tag{ 2 }$$

  If Equation (2) evaluates true, the event proceeds to the final step in the pipeline.
• 4.  
  In this step, the events flagged in step 3 are compared between the cameras. An event which occurs at the same time in both time series is marked as an occultation candidate.

The EMCCD cameras provide GPS-based time stamping at the image collection point, allowing the simultaneity of events to be determined with high precision. However, the Colibri pipeline evaluates matches between the time series with a 150 ms (∼3 frame) tolerance, as further testing is required to establish confidence in the EMCCD time stamping at the microsecond level. This three-frame tolerance also acknowledges the slight differences in observed light curves that may arise from the differing geometries of the two camera systems.

The Colibri algorithm provides a heuristic determination of occultation candidates. To quantitatively establish the confidence level of a given detection, we must benchmark the performance of this algorithm for simulated data sets. Retrieval rates at each threshold are discussed in Section 6.1, with false positive rates determined in Section 5.

We will consider the noise in our time series in terms of two components: the Poisson-distributed shot noise resulting from the photons measured from the star (${\sigma }_{P}$), and all other noise (${\sigma }_{O}$). We can therefore write

$$\begin{eqnarray}&&\sigma =\sqrt{{\sigma }_{P}^{2}+{\sigma }_{O}^{2}},\end{eqnarray} \tag{ 6 }$$

where σ is the standard deviation of our light curve.

We can also calculate ${\sigma }_{P}$:

$$\begin{eqnarray}&&{\sigma }_{P}=\sqrt{\displaystyle \frac{N}{G}}\ast G.\end{eqnarray} \tag{ 7 }$$

In Equation (7), N is the mean counts in the time series and G is the EMCCD detector gain (G = 100 for EMCCD CAMO). Using Equations (6) and (7), we can determine ${\sigma }_{O}$ from measurements of N and σ.

During an occultation event, there is a decreased number of photons arriving at the detector from the star. To include noise in the kernels of our kernel set, we must consider this variability.

For each kernel, we replace each point with a random variate drawn from a Poisson distribution with $\lambda =\tfrac{N\ast k}{G}$, where k is the fractional flux at that point in the kernel. This variate is then multiplied by G, as our light curves are in units of counts. Finally, we apply ${\sigma }_{O}$ using a random Gaussian distribution. In Figure 7, we illustrate that our two-component noise model closely reproduces the noise effects in our observed light curves.

To investigate the Colibri pipeline's event retrieval rates, noise-convolved versions of each kernel in the kernel set were injected into 14,456 5-minute time series. These series correspond to the stars that appeared in EMCCD CAMO data taken on 2016 August 1st and passed the acceptance criteria outlined in Section 3.1.

After injection, the Colibri pipeline attempted to retrieve these simulated events. Table 3 benchmarks the performance of this method at each of the four steps described in Section 3.2.

The first five columns of Table 3 identify each kernel by its occultation parameters. As the first step in occultation detection, the pipeline selects the most significant flux dip in each 5-minute time series after the Mexican-hat convolution (Section 3.2). However, depending on the parameters of the simulation kernel, this dip does not always correspond to the simulated event. The fraction under the Correct Dip Detected column in Table 3 shows how frequently the injected event was the one identified as the most probable occultation candidate in the time series. The entries in the subsequent columns show the cumulative fractions of originally simulated events that continue to be recovered after each discriminating step in the detection pipeline.

Retrieval rates from the final dual-detection step of our event retrieval simulation range from ∼2.5% to ∼92% (Table 3), illustrating a significant dependence on the occultation parameters, including KBO size, stellar angular diameter, and impact parameter. Kernels 38, 40, and 42, the templates with the lowest successful retrieval rates in Table 3, model small KBOs (1.5 km) occulting angularly large stars (0.08 mas) with large impact parameters (750 m).

Table 3 shows that these kernels perform particularly poorly at the first threshold: identification as a candidate for further analysis using the Mexican-hat convolution. Due to scintillation and other noise, these ∼25% dips often cannot be distinguished from non-occultation fluctuations. As the survey's sensitivity to these type of events is already so low, we calculate that removing these kernels from the kernel set for a 17 Hz-sampled, 5.8 cm aperture survey would reduce the false positive rate by a greater factor than it would reduce occultation sensitivity.

As discussed in Section 1.3, the expected occultation rate varies greatly depending on the index q of the small-object power law. This rate also has a significant dependence on KBO size. Using Equation (1) and related equations from Bickerton et al. (2009) and assuming q = 1.9 (Fraser & Kavelaars 2008), we calculate that the number of occultations by 1500 m to 2750 m objects is approximately 1.7× the number of occultations by 2750–5000 m objects.¹¹

If the nine kernels with final retrieval rates less than 10% were removed from the kernel set, the pipeline's final false positive rate would decrease fourfold, from 0.002% to 0.0005%. Above, we calculated that approximately two thirds of the occultations that may be detected by this kernel set correspond to objects smaller than 2750 m. Therefore, removing these 1500 m kernels from the kernel set could result in a threefold decrease in KBO sensitivity at worst. Realistically, the sensitivity decrease would be less severe: only a fraction of kernels sensitive to 1500–2750 m objects would be removed, and the removed kernels have much lower retrieval rates (and therefore are sensitive to a smaller fraction of the events they are theoretically capable of detecting) than the other kernels.

Nevertheless, it would be beneficial to lower the minimal detectable diameter of KBOs and improve retrieval rates at these small diameters rather than simply removing the low-sensitivity kernels. The number of occultations by objects between 500 and 1500 m is 2.6x the number of occultations by objects between 1500 and 5000 m: lowering the minimal detectable diameter would greatly improve the likelihood and increase the frequency of occultation detections. We aim to decrease the minimal detectable diameter of KBOs with the main Colibri campaign, which will utilize much larger (50 cm) telescopes.

In addition, the main campaign will use 40 Hz sampling as opposed to the 17 Hz sampling of these preliminary trials. At a 40 Hz sampling rate, characteristic diffraction patterns in the time series are a strong differentiator of occultation events from random noise or scintillation (Figure 6). From the data in Table 3 and the false positive rates described in Section 5, the kernel match step may not appear to be a significant benefit to the detection pipeline, as diffraction patterns are not clearly visible in the data with 17 Hz sampling (Figures 4 and 5). We anticipate a substantial improvement in retrieval rates at the kernel matching step with the 40 Hz main Colibri campaign.

Zoom In Zoom Out Reset image size

The kernel set will need to be expanded and reoptimized prior to main campaign data collection to ensure that the 40 Hz templates with visible diffraction patterns are sensitive to the full spectrum of parameter values, including the lower KBO diameters available to the Colibri array. Based upon previous studies, a set of 300 or more kernels may be required to ensure sensitivity throughout the parameter space for a survey conducted at a 40 Hz sampling rate with a noise level of 2%–4% (Bickerton et al. 2008).

Using a Monte Carlo simulation based on the diffraction modeler described in Section 4, we are able to estimate the parameters that generate a given occultation pattern. This process allows KBO characteristics to be determined in the event of a successful detection.

Zoom In Zoom Out Reset image size

Additionally, this method may be used to filter out false positive candidates that were not removed by the main Colibri pipeline. By constraining the occultation parameters for the dips detected by each camera, discrepancies between the two curves can be analyzed quantitatively. Disparate curves can then be flagged as false positives.

Zoom In Zoom Out Reset image size

Figure 8 shows the parameter estimates for an occultation in a high-noise time series (S/N = 6). The original input parameters were a 2000 m KBO, a point source star (0 mas stellar angular diameter), 0 m impact parameter, 0 shift adjustment, and a distance of 40 au, with the kernel corresponding to these parameters displayed in Figure 9. Despite this time series exhibiting even higher noise levels than allowed by Colibri's acceptance criteria, the histograms in Figure 8 show that for each parameter except distance, the simulation results are clustered around the input value. Distance cannot be well constrained with 17 Hz data.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Some degeneracy can be observed between KBO size and impact parameter. This is unsurprising, as a larger KBO deepens the occultation curve while an increased impact parameter shallows it. In the case of Figure 8, this may account for the increase in estimated KBO size and impact parameter as compared to the simulation input. This degeneracy is therefore an effect to consider when evaluating the results of this method on non-simulated data.

From the 4000 star hours of dual-camera data obtained in the present experiment (Section 2.3), the Colibri pipeline flagged one candidate occultation of a 6.72 V-magnitude star (S/N = 7) on February 26 (Figure 10).

Using the Monte Carlo parameter retrieval method described in Section 6.2, we found that the occultation in one time series corresponded to a KBO of radius 1800 m ± 500 m, while we determined a KBO radius of 600 m ± 200 m using the second time series. The two estimates of the radius are marginally inconsistent: a t-test results in a p value of 0.026, so the estimates are inconsistent at the 97.4% level. For this reason, we classify the candidate as a likely false positive. The impact parameter estimate is the same within uncertainties for the two curves, which suggests that the size/impact parameter degeneracy is not contributing to the discrepant KBO size estimates.