Orbital Foregrounds for Ultra-short Duration Transients

We observed 3,372,044 high-probability candidates that passed our vetting criteria as described in Section 2.2. Of these, we identify 1,415,722 candidates that do not appear in multiple epochs as likely satellite flashes. This cut removes variable stars and persistent artifacts from bright stars.

Single-epoch candidates tend to occur in tracks across the sky on timescales ranging from subcadence to hours. Figure 1 shows a typical track observed by Evryscope, followed for 1 hr, along with typical 30 × 30 pixel (6.6 × 66) subtraction stamps from the EFTE pipeline. The timescale of the delay between flashes gives an angular speed of 10'' s⁻¹, or roughly one Evryscope pixel per second. However, we do not observe any streaking at any epoch, implying that the duration of each individual flash is much less than 1 second. This is consistent with the population of fast optical flashes noted in Biryukov et al. (2015) and Karpov et al. (2016), but observed in images integrated over minutes.

Zoom In Zoom Out Reset image size

EFTE uses a simplified image-subtraction technique for candidate detection. Evryscope focus and optics, and thus PSFs, are stable on month-to-year timescales (Ratzloff et al. 2020), and atmospheric seeing is dominated by optical effects under all observing conditions in each Evryscope's 132 pixels. This stability means that image subtraction within a single pointing does not require the PSF-matching techniques addressed by standard routines, such as HOTPANTS (Becker 2015) or ZOGY (Zackay et al. 2016). Similarly, because we are targeting the fastest-timescale events, we do not require widely spaced reference frames; instead, an earlier image from the same pointing is subtracted from each image, typically with a 10 minute separation.

For each single-camera science image S(x, y) and reference image R(x, y), we calculate a discovery image D(x, y) from the per-pixel change in signal-to-noise ratio. Both S(x, y) and R(x, y) are calibrated, background-subtracted, and aligned images in electron units. We measured image background levels using an interpolated and clipped mesh, as implemented in sextractor (Bertin & Arnouts 1996). The discovery image D(x, y) is given by:

$$\begin{eqnarray}&&D(x,y)=\displaystyle \frac{S(x,y)}{{s}_{S}(x,y)}-\displaystyle \frac{R(x,y)}{{s}_{R}(x,y)},\end{eqnarray} \tag{ 1 }$$

where s_R(x, y) and s_S(x, y) are noise images calculated for each image I(x, y) with measured background noise s_B²(x, y) as

$$\begin{eqnarray}&&{s}_{I}(x,y)=\sqrt{I(x,y)+{s}_{B}^{2}(x,y)}.\end{eqnarray} \tag{ 2 }$$

While this technique is not statistically optimal, it is efficient (98.5% of images are reduced in cadence) and produces stable artifacts that can be rejected via automated means. We identified all sources where D(x, y) ≥ 3 in at least three contiguous pixels for automated vetting. We perform aperture photometry for all candidates using the science image, and calibrated the results to the ATLAS All-Sky Stellar Reference Catalog (Tonry et al. 2018b).

A typical single image subtraction using this method produces O(10³) candidates. We require the ratio of negative-to-positive pixels within a 5 pixel diameter to be less than 0.3 to remove dipole artifacts from misalignments and Poisson noise peaks near bright stars, which reduces the number of candidates to O(10²). These are then processed with a convolutional neural network (LeCun et al. 1998) trained on small cutouts from S(x, y), D(x, y), and R(x, y) in both real and simulated data, reducing the final number of candidates to ${6}_{-5}^{+13}$ per image when averaged over all weather conditions and cameras.

We characterized the completeness of the EFTE survey with injection-recovery testing. The image sample contained 500 pairs with the same 10 minute separation used on-sky, and was representative of all weather and instrument conditions contained in the survey. We injected ∼1200 sources into each image ($8.0\leqslant {m}_{g^{\prime} }\leqslant 16.75$), using the normalized median of nearby stars as the injection PSF. We processed each image pair as described in Section 2.2 twice, alternating which image in the pair was used as the science frame. We define a recovery as a candidate detected within three pixels of an injection position, with the radius determined by the 99th percentile of the nearest-neighbor distance between the vetted candidates and their nearest injection position.

Figure 2 shows the measured completeness as a function of magnitude, along with the 90% confidence interval (CI), which includes measurement uncertainty, variability with observing conditions, and the vignetting effects. Completeness decreases at both the bright and dim ends of the distribution. Saturation causes poorly constrained centroids leading to nonrecovery at the bright end. The average 50% completeness limit of ${m}_{g^{\prime} }=14.2$ is brighter than the Evryscope dark-sky, single-image limit of ${m}_{g^{\prime} }=16$ due to a combination of including all-sky conditions in our analysis and using reference images with similar noise profiles to science images.

Zoom In Zoom Out Reset image size

Despite careful vetting, the median EFTE false positive rate is ∼6800 per hr. In EFTE surveys for rapidly evolving transients, further vetting is performed by crossmatching with galaxy and stellar catalogs, or by association with gamma-ray and gravitational-wave transient skymaps, combined with direct human interaction. Measuring an all-sky event rate, however, requires knowing the fraction of real events in the sample. To measure this fraction, we visually inspected 27,817 randomly selected single-epoch candidates that passed our automated vetting as described in Section 2.2.

Each candidate was assigned a binary classification, where a real classification corresponds to a candidate that is morphologically consistent with an astrophysical transient. Of the candidates classified, 2823 were classified as real. From this, we estimate the real flash fraction (RFF) in the vetted EFTE event stream to be 0.10 ± 0.03. Least-squares fits to the RFF as a function of magnitude and solar elongation produced slopes not different from zero, with an uncertainty set by the standard deviation of the residuals around the fit. Bogus candidates included:

• 1.  
  Subtraction artifacts from bright stars (50%).
• 2.  
  Optical ghosts, distinguished based on non-PSF-like shapes and presence in the reference frame (1%).
• 3.  
  Aircraft strobes, distinguished based on nearby parallel streaks (3%).
• 4.  
  Particle strikes (46%).

The final category includes both readily identifiable cosmic-ray muon tracks and signals caused by Compton recoil electrons from environmental radioisotopes, which can be PSF-like. To constrain this population, we reduced a series of 120 s darks with EFTE and searched for candidates meeting our vetting criteria. We place an upper limit on the base rate of PSF-like particle strikes of ≤0.1 per image, or ∼300 sky⁻¹ hr⁻¹, which is small compared to the all-sky orbital flash rate (Section 3.2).

Evryscope PSFs, and the resulting degree to which stellar sources are undersampled, are variable across each individual camera's FOV. At the center of the field, dim sources can have an FWHM ≪ 1 pixel. To avoid distortion of these sources, we estimate the prevalence of particle strikes in our candidate sample using the RFF, rather than directly removing them from the images with standard cosmic-ray mitigation tools, such as LA-Cosmic (van Dokkum 2001).

We used a Monte Carlo approach to model the effects of completeness and reliability as a function of magnitude (Section 2.3). We modeled the per-magnitude completeness as a bounded Johnson distribution (Johnson 1949) using maximum likelihood estimates for the mean, standard deviation, skewness, and kurtosis from injection testing. We corrected for the reliability of the event sample using the RFF described in Section 2.4, modeled per magnitude bin as a normal distribution.

For each magnitude bin, we made 100,000 draws from the fitted completeness and reliability distributions, each time calculating a corrected event rate as

$$\begin{eqnarray}&&{\rm{\Gamma }}{\left({m}_{o}\right)}_{i}=\displaystyle \frac{{r}_{{m}_{o}}{R}_{i}}{{C}_{i}},\end{eqnarray} \tag{ 4 }$$

where C and R represent values drawn from the completeness and reliability distributions, respectively. The reported mean rate is the 50th percentile of Γ(m_o)_i. The lower and upper bounds of the 90% CI were taken to be the 5th and 95th percentiles of Γ(m_o)_i.

As noted in Lyutikov & Lorimer (2016), flashes shorter than an image exposure time (typically ≫1 second) will exhibit phase blurring, diluting their flux relative to the surrounding stars by the ratio of the integration time to their intrinsic duration. For subsecond durations, consistent with the example seen in Figure 1 and in Biryukov et al. (2015), the peak of the flash light curve will be brighter than observed in long integrations. In general, the peak magnitude of the flash m_P is given by

$$\begin{eqnarray}&&{m}_{P}=-2.5{\mathrm{log}}_{10}\left(\displaystyle \frac{{T}_{{\exp }}{10}^{-0.4{m}_{o}}}{{\tau }_{f}}\right),\end{eqnarray} \tag{ 5 }$$

where τ_f is the equivalent width of the light curve and T_exp is the exposure time. We assume a flash duration of 0.4 s, based on the mode of the distribution presented in Karpov et al. (2016), when estimating the peak brightness from the observed magnitude.

Figure 3 shows the cumulative and normalized magnitude density distributions as a function of observed and estimated peak brightness. The shaded regions represent the 90% CI. For the cumulative distributions, the CI is bounded by the cumulative quadrature sum of the per-bin CI limits. We extrapolate an all-sky event rate from the observed rates in three regions: around the south celestial pole (SCP), within 10° of the equator, and across all declinations. Both the all-sky and equatorial distributions peak at m_o = 13, whereas the polar distribution peaks at m_o = 12.2, with a possible second peak near the saturation limit.

Zoom In Zoom Out Reset image size

We measured integrated flash rates for m_o < 14.25 of ${1.0}_{-0.15}^{+0.27}\times {10}^{3}$, ${4.0}_{-0.60}^{+1.40}\times {10}^{3}$, and ${1.8}_{-0.28}^{+0.60}\times {10}^{3}$ sky⁻¹ hr⁻¹ in the SCP, equatorial, and averaged all-sky regions, respectively. Assuming that the observed population is dominated by satellite glints, we expect a higher rate near the equator because equatorial orbits are confined to a narrow decl. band. The SCP region will only contain objects in polar orbits spanning the full decl. range, lowering the observed rate.

If the observed flashes are caused by reflections from satellites, none should be expected from the region of the sky covered by Earth's shadow. The solid angle subtended by the shadow depends on satellite altitude, following the geometry illustrated in Figure 4. For low-Earth orbit (LEO) altitudes (<2000 km), the shadow covers a 50° radius around the solar antipode. The angular size shrinks to a 14° radius for medium-Earth orbit (MEO), or a 9° radius for geosynchronous orbit. We evaluated the distance between each of the human-vetted candidates from Section 2.4 and the solar antipode. Shaded regions depict the angular extent of Earth's shadow for LEO, MEO, and geosynchronous satellites.

Zoom In Zoom Out Reset image size

The prevalence of flashes decreases steadily with proximity to the shadow in the region covered for LEO objects, before falling in the solid angle covered for MEO and higher orbits. Approximately 34% of the flashes occurred within 50° of the center of the shadow. Few flashes occur within the region within Earth's shadow for MEO and geosynchronous orbits, with 3.5% within 14° of the antipode and only 1.1% within 9°, which suggests that the majority of the flashes are generated by satellites in middle- and high-Earth orbit. No Evryscope observations occur within ∼80° of the Sun.

The all-sky magnitude distribution in Figure 3 shows that the instantaneous peak brightness of many flashes should be detectable to the unassisted human eye. For a typical suburban sky with limiting magnitude V ≈ 4, we predict a naked-eye-visible event rate of ${340}_{-85}^{+150}$ sky⁻¹ hr⁻¹ based on the all-sky averaged cumulative event rate of ${1.8}_{-0.28}^{+0.60}\times {10}^{3}$ sky⁻¹ hr⁻¹ (m_o < 14.25). At the darkest sky sites, with limiting magnitude V ≈ 6, this rate increases to ${740}_{-140}^{+200}$ sky⁻¹ hr⁻¹. Within the equatorial region, the expected naked-eye event rate increases to ${840}_{-220}^{+390}$ sky⁻¹ hr⁻¹ (V ≈ 4), or ${1800}_{-350}^{+500}$ sky⁻¹ hr⁻¹ (V ≈ 6), based on the cumulative equatorial event rate of ${4.0}_{-0.60}^{+1.40}\times {10}^{3}$ sky⁻¹ hr⁻¹.

The time resolution (Cornsweet 2014) and biochemical adaptation time (Yeh et al. 1996; Dunn & Rieke 2006) of the human eye are comparable to orbital flash durations. We expect the naked-eye detectability of these flashes to vary strongly with observer ability.

The impact of orbital flashes on narrow-field imagers is negligible. For a 10 × 10' FOV, we predict approximately 1.2 flashes per 1000 hr of exposure time based on the all-sky averaged event rate, or 2.7 flashes per 1000 hr of exposure time based on the equatorial region event rate. While these are relatively rare events, they could account for occasional flashes that have been remarked on by amateur astronomers and some rejected single-image detections in tiling sky surveys.

The flash rate measured here also implies a high coincidence rate for multimessenger events. The Canadian Hydrogen Intensity Mapping Experiment (CHIME) is able to localize FRBs to the nearest arcminute (The CHIME/FRB Collaboration et al. 2018), reducing the expected flash count within the error radius to a likely acceptable $3.8\times {10}^{-5}$ hr⁻¹ based on the all-sky averaged event rate. Wide-angle surveys like Evryscope can be expected to have a false alarm rate (FAR) of 1 per 3 yr for apparent FRB optical counterparts due to orbital flashes. For FRBs localized to the equatorial region, the expected event rate and FAR increase to 8.5 × 10⁻⁵ hr⁻¹ and 1 per 1.3 yr, respectively.

In contrast, gamma-ray bursts (GRBs) from the Fermi Gamma Burst Monitor have a median 90% localization area of 209 square degrees (Goldstein et al. 2020), leading to an expected flash count of 9.1 hr⁻¹ based on the all-sky averaged cumulative event rate. With simultaneous coverage and physical constraints on the timescale of the optical component, the expected flash count drops to 0.15 minute⁻¹. For GRB localizations within the equatorial region, the rates increase to 20 hr⁻¹, or 0.33 minute⁻¹.

Events with larger localization regions or weakly constrained early lightcurves, like gravitational wave events from LIGO/Virgo, will be more heavily impacted. For a typical 1200 sq. degree sky map, the expected flash is 53 hr⁻¹ assuming the all-sky averaged event rate, or 120 hr⁻¹ within the equatorial region. The resulting FAR will increase linearly with both localization area and trigger rate.

We predict that point-like flashes will occur in images from Vera C. Rubin Observatory at a rate of 6.4 × 10⁻⁴ flashes per image (15 s of 3.5 square degrees), based on the all-sky averaged event rate. Assuming O(1000) pointings, we expect Rubin to observe on the order of 1.3 flashes per night from this population. In the worst case scenario of all pointings being confined to the equatorial region, this rate increases to 1.4 × 10⁻³ flashes per image, or 2.8 flashes per night. Due to the 15 s exposure time of Rubin, the average observed magnitude of the flash distribution will be shifted to m_g ∼ 10.7, 100× brighter than the g-band bright limit of m_g = 15.7 (LSST Science Collaboration et al. 2009). Our observed magnitude distribution drops off at the faint end, but there is the potential for a second distribution to exist beyond the Evryscope depth limit, as is seen for geosynchronous debris observed by the DebrisWatch survey (Blake et al. 2020).

We estimate that a worst-case scenario 10 millisecond flash from a geosynchronous satellite will produce a 015 streak, likely indistinguishable in Rubin's 02 pixels. Longer-duration flashes, like the 0.4 s events we have assumed here, would produce obvious 6'' streaks. However, fully constraining the expected morphology of orbital flashes in Rubin data will require modeling of their duration via orbital characteristics of the population.