© ESO 2020

1 Introduction

The Draconids are an annual meteor shower whose parent body is the Jupiter-family comet 21P/Giacobini-Zinner. The shower usually has very low activity with a zenithal hourly rate (ZHR) of ~1 (Jenniskens 2006). It produced large meteor storms in 1933 and 1946 and strong outbursts in a number of other years (Egal et al. 2019). Many of these outbursts were not predicted beforehand. The occasionally high intensity and unpredictability of the shower have made it a focus of research, particularly as it can pose a significant impact risk to spacecraft in the near-Earth environment (Beech et al. 1995; Cooke & Moser 2010; Egal et al. 2018).

In recent years, the shower has produced several notable outbursts. In 2011 an outburst was predicted in advance (Watanabe & Sato 2008; Maslov 2011; Vaubaillon et al. 2011) and was well observed using both radar (Kero et al. 2012; Ye et al. 2013a) and optical methods (Trigo-Rodríguez et al. 2013; Borovička et al. 2014; Koten et al. 2014; Šegon et al. 2014; Vaubaillon et al. 2015). That year the outburst reached a ZHR of 350, and had an average mass index of 2.0 ± 0.1, but this latter varied between 1.84 and 2.30 over 1h during the peak (Koten et al. 2014).

In contrast, the 2012 outburst was not predicted and was only well observed by the Canadian Meteor Orbit Radar (CMOR) (Brown & Ye 2012). The shower produced a meteor storm at radar sizes (ZHR ≈ 9000 ± 1000), but visual observers reported a ZHR almost two orders of magnitude lower, ZHR ~ 200, suggesting a high mass-distribution index, that is, the stream was rich in small meteoroids. Unfortunately, due to the timing of the peak (maximum over central Asia) and unfavorable weather conditions elsewhere, no optical orbits associated with the outburst were secured. The 2012 outburst was also peculiar in that modeling suggested a very high (>100 m s⁻¹) meteoroid ejection velocity from the parent comet was needed (Ye et al. 2013b).

The 2018 outburst was predicted by various authors (Kresák 1993; Maslov 2011; Ye et al. 2013b; Kastinen & Kero 2017; Egal et al. 2018), but the predicted activity varied from weak (ZHR 10–20; Maslov 2011) to possible meteor storm levels (Kastinen & Kero 2017). The most recent work by Egal et al. (2018), which accurately reproduced most historic Draconid activity, predicted a peak ZHR of ~ 80 at 00:00 UTC on October 9, 2018.

In this work we analyze 1.5 h of optical observations from Southwestern Ontario just after the peak¹ of the 2018 outburst, from 00:00 UTC to 01:30 UTC. We also develop a novel method of population and mass index estimation, and compute these indices using our observations. Finally, we compare model-predicted radiants with our multi-station observations and compute the shower flux.

2 Instruments and observations

Optical observations were conducted using low-light-level video systems in Southwestern Ontario, Canada, operated by the University of Western Ontario Meteor Physics Group. The observations started shortly after local sunset between October 8 and 9 at around 00:00 UTC and ended at 01:30 UTC due to cloudy weather – thus only a short period after the peak was observed.

The double station observations were made using four Electron Multiplying Charge Coupled Device (EMCCD) Nüvü HNü 1024 cameras² with Nikkor 50 mm f/1.2 lenses. The systems were operated at 32 frames per second, had a limiting stellar magnitude of + 10.0^M, and a field of view of 15×15 degrees.

Despite clear skies in the region, none of the dozen Southern Ontario Meteor Network (SOMN) all-sky cameras (Brown et al. 2010)observed a single Draconid during all of October 8–9. These all-sky cameras have an effective limiting meteor magnitude of − 2^M, qualitativelysuggestive of an absence of larger meteoroids in the falling branch of the outburst. The EMCCD cameras did observe many Draconid meteors, which we use in the following analysis. In this section we describe the EMCCD hardware and give details of the data-reduction procedure.

2.1 Electron Multiplying Charge Coupled Device systems

The NüvüHNü 1024 EMCCD cameras are the most recent addition to the Canadian Automated Meteor Observatory (CAMO) systems (Weryk et al. 2013). These 16-bit cameras use electron multiplying technology to increase the number of electrons accumulated in the registers before they are amplified, digitized, and read out, which greatly increases the signal-to-noise ratio (S/N). In combination with a fast 50 mm f/1.2 lens, the single-frame limiting magnitude at 32 frames per second for stars is around + 10.0^M, and + 7.5^M for meteors. The native resolution of the camera is 1024 × 1024 pixels, but the frames are binned to 512 × 512 to make the higher frame rate possible, which allows better trajectory precision. The field of view is 14.75 ° × 14.75°, and the pixel scale is 1.7 arcmin pixel⁻¹. All cameras were run at a gain setting of 200. These systems are more sensitive and have less noise than the previous generation of systems which relied on Gen III image intensifiers (Campbell-Brown 2015). The cameras each have external Global Positioning System (GPS) timing directly encoded at the frame acquisition stage. Figure 1 shows a sample meteor captured with the EMCCD camera and Fig. 2 shows the light curve of the same meteor. We highlight the agreement between the sites in range-corrected magnitudes (normalized to 100 km). In Appendix A we give more examples of multi-station Draconid light curves observed by the EMCCD systems on Oct 9, 2018.

Two cameras are deployed at each of the two CAMO sites in Southwestern Ontario. Cameras 01F and 01G are located at Tavistock (43.26420° N, 80.77209° W, 329 m), and 02F and 02G at Elginfield (43.19279° N, 81.31565° W, 324 m). The distance between the stations is about 45 km. Both pairs are pointed roughly north, one pair (the G cameras) at the elevation of 45°, the other (F cameras) at 65°. There was no volume overlap between the F and G cameras. One camera from each site had the Draconid radiant inside the field of view. The EMCCDs captured a total of 92 double station Draconids ranging in peak magnitudes from + 1^M to + 6.5^M, of which only 68 had complete light curves from at least one camera.

2.2 Datareduction and calibration

Using the tools available in the open-source RMS (Raspberry pi Meteor Station) library³, automated meteor detection (described in detail in Vida et al. 2016) was run on all collected data and both the astrometric and photometriccalibrations were performed manually. The automated detections were manually classified and only Draconids were extracted for further analysis.

The astrometric and photometric calibrations were done using the Gaia DR2 star catalog (Gaia Collaboration 2018) as the Gaia G spectral band matches the spectral response of the EMCCDs well. We applied flat fields and bias frames to the video data using the method described in Berry & Burnell (2000). The average residuals of the astrometric calibration were ~ 1∕10 pixels or ~ 10 arcsec, and the accuracy of the photometric calibration was ± 0.1 mag. Figure 3 shows the photometric calibration, and Table C.1 details the astrometric calibration for the 01F EMCCD camera.

Fig. 1 Stack of raw frames showing a + 1M Draconid recorded with camera 01G on October 9, 2018 at 00:01:48 UTC. The meteor was first detected at an apparent magnitude of + 7.0M. V -band magnitudes of a selection of stars are shown as well. The vertical gradient of the image background is intrinsic to this specific camera and is compensated for by flat field correction during data reduction. The two thinner and fainter lines left of the meteor are satellites.

Fig. 2 Absolute Gaia G-band magnitude light curve of the meteor in Fig. 1. The error bars are ~ 0.1 magnitude. Observations from camera 01G (Tavistock) are shown in blue, and observations from camera 02G (Elginfield) are shown in red.

Fig. 3 Photometric calibration of the 01F camera. The uncalibrated magnitude is simply − 2.5 log10Is, where Is is the background-subtracted sum of pixel intensities of all star pixels (LSP in the figure). The sensor response was linear (Gamma = 1.0). The magnitudes of stars used for calibration ranged from + 2.3M to 9.9M.

3 Maximum likelihood estimation method of computing population and mass indices

Population and mass indices are essential for computing meteor shower fluxes and investigating their evolutionary history. It is usually assumed that the cumulative number of meteor peak magnitudes and the cumulative logarithm of the number of meteoroids above a certain mass threshold follow single power-law distributions.

Brown & Rendtel (1996) define the population index r as: $$r=\frac{N(M_v+1)}{N(M_v)},$$(1)

where N(M_v) is the number of meteors with peak magnitude M ≤ M_v.

Assuming that dN is the number of meteoroids with a range of masses between m and m + dm, McKinley (1961) defines the mass index s as: $$\text{d}N\propto m^{-s}\text{d}m.$$(2)

The population and the mass index are the slopes of the best fit lines on a log–log (magnitude or logarithm of mass vs. logarithm of the number of meteors) cumulative distribution of magnitudes or masses. These lines can be simply parameterized as: $$f(x)={10}^{ax+b},$$(3)

where x is either the magnitude (in case of the population index) or the logarithm of the mass (mass index). We note that we assume that the cumulative distribution is normalized to the [0, 1] range.

Becauseall observation systems have some limiting sensitivity, the number of observed meteors is finite and detection efficiency drops off rapidly at magnitudes close to the limiting sensitivity. Thus there is only one part of the parameter space in which the power-law assumption is valid. As a result, operationally a power law is only fit up to some limiting sensitivity x_min. A similar approach is used in other fields as well, although ways of estimating the value of x_min differ (Corral & González 2019). Pokornỳ & Brown (2016) suggested fitting two parameters x_min, x_max which describe the range of magnitudes or radar echo amplitudes for which the linear approximation is valid. Pokornỳ & Brown (2016) argue that the addition of the upper boundary x_max is necessary as small-number statistics may skew the power law at larger amplitudes or masses. They used MultiNest (Feroz et al. 2009)to fit their model, an advanced Bayesian regression method which can also perform robust error estimation.

In many studies (e.g., Brown et al. 2002), the meteoroid mass or population index is found by constructing a cumulative histogram and performing linear regression on binned data. This approach has many serious problems, as detailed in Clauset et al. (2009), the most important of which is that the result may change with histogram binning, especially when data sets with a small number of meteors are used. Molau et al. (2014) proposed a more advanced method capable of estimating the population index and the meteor flux at the same time. The drawback of this latter method is that it requires a running estimation of the meteor limiting magnitude from observations. Furthermore, it is constructed such that estimating the exact value of the population index is rather subjective, and no details of the error analysis are provided.

3.1 Description of the new method

To estimate the population and the mass index robustly, in a statistical manner, and with a reliable error estimate, we use the maximum likelihood estimation (MLE) method to fit a gamma distribution to the observed distribution of magnitudes and/or masses. This approach was inspired by the work of Clauset et al. (2009) where they question the validity of the power-law assumption in many real-world applications.

First, we tried to fit a Gumbel distribution to the mass distribution following Blaauw et al. (2016), but we found that it always produced fits with smaller p-values than the gamma distribution; under the null hypothesis the data and model distributions are identical.

Our choice of the gamma distribution for our probability distribution has a theoretical basis. If we assume that meteors appear at a rate λ per unit time, then the time between meteors follows an exponential distribution with rate λ while the total number of meteors in a given time span t follows a Poisson distribution P(λt). A gamma distribution Gamma (n, λ) is the sum of independent exponentially distributed random variables, and describes the length of time needed to observe a number of meteors n (Akkouchi 2008; Leemis & McQueston 2008).

By replacing the time in this analogy with magnitude or mass, the theory still applies. Given a meteor rate per unit magnitude or mass (i.e., a population or mass index), a gamma distribution describes the distribution of magnitudes or masses one expects given a certain total number of meteors. In the special case where the rate is a constant in our chosen time interval, the sum can be simplified to an Erlang distribution. Because historical observations show that the mass index is not constant even on short timescales for certain showers (Koten et al. 2014), we use the more general gamma distribution.

We can show the above argument is correct quantitatively. We note that in the paragraph above we implicitly assume an exponential distribution, not a power law. The power law assumption for the distribution of meteoroid masses used in meteor astronomy comes from the foundational theoretical paper by Dohnanyi (1969) who assumed that when asteroids and/or meteoroids collide, the resulting fragment mass distribution follows a power law (Clauset et al. 2009): $$p(x)=\frac{s-1}{x_{\text{min}}}{\left(\frac{x}{x_{\text{min}}}\right)}^{-s},$$(4)

where s is the scaling parameter (differential mass index in our example) and x_min is the lower bound of the power-law behavior. The power-law assumption continues to propagate through the analytical derivation of Dohnanyi (1969), resulting in a population-wide differential mass distribution index of s = 1.834, assuming collisional steady-state conditions. Dohnanyi (1969) justified the use of the power law by appealing to the experimental work of Gault et al. (1963), who showed that the cumulative mass–size distribution of fragmented rock from craters formed by hyper-velocity impacts in the lab follows a power-law, a result confirmed in more recent experiments and simulations (e.g., Kun & Herrmann 1996).

In reality, the number of meteors naturally tapers off as the limiting system magnitude is approached. This can be modeled as a power-law distribution with an exponential cutoff. Its probability density function p(x) is given by Clauset et al. (2009): $$\begin{array}{l}p(x)=Cf(x)\\ f(x)=x^{-s}{\text{e}}^{-\lambda x}\\ C=\frac{{\lambda }^{1-s}}{\Gamma (1-s)},\end{array}$$

where λ is the decay rate, f(x) the basic functional form, C the normalization constant such that ${\displaystyle {\int }_0^{\infty }C}f(x)\text{d}x=1$, and Γ is the gamma function. If we define a shape parameter as k = 1 − s, it can easily be shown that the resulting probability distribution is the gamma distribution: $$p(x)=\frac{{\lambda }^k{x}^{k-1}{\text{e}}^{-\lambda x}}{\Gamma (k)}.$$(8)

Figure 4 shows a comparison of the cumulative density functions (CDFs) of a power law and those of a gamma distribution. Both distributions follow the same trend at smaller values, but the exponential tapering overwhelms the power-law component in the gamma distribution at around x ~ 10⁻¹ in this example. To estimate the power-law scaling parameter from the gamma distribution, one can compute the slope of the gamma CDF before the tapering becomes significant. In theory, one may derive the slope directly as s = 1 − k, but in practice (especially when one works on a small data set) we find that this does not give robust estimates ofpopulation and mass indices. It appears that the fit sometimes prefers larger values of s which it compensates for by using larger values of λ. For that reason, we develop a more robust method described below.

We use the gamma distribution MLE fit from the functions in the Python scipy library⁴. We refer readers to the documentation of the scipy library for details. The Python source code of our new method can be found in the Western Meteor Python Library⁵ under wmpl/Analysis/ FitPopulationAndMassIndex.py.

Figure 5 shows an example of the gamma distribution fit to meteor magnitude data. After the fit is performed, the inflection point in the probability density function (PDF) is found (the minimum of the first derivative of the PDF). This inflection point is taken as the limiting sensitivity at which data are complete; this means that the number of meteors has started decreasing at this point purely due to sensitivity losses. In our approach, the inflection point is used as an anchor; it indicates that the optimal position where the slope of the population or mass index should be estimated is at brighter magnitudes or larger masses fromthat point.

We operationally define the optimal reference point x_ref for: the population index, as one magnitude brighter than the inflection point; and for the mass index, as 0.4 dex larger in mass than the inflection point – a 0.4 dex in mass roughly corresponds to the difference of one meteor magnitude (Vida et al. 2018).

Assuming an 8-bit camera is used to capture meteor data, it gives log_2.512 256 ≈ 6 magnitudes of dynamic range. As shown by the examples below, the data are often not complete for the faintest three magnitudes, leaving only three magnitudes before reaching the point of saturation. One magnitude brighter than the completeness magnitude is a rough midpoint between the point of completeness and sensor saturation, allowing for some leeway in case the saturation point is shifted because of different meteor angular speeds.

The mass index s or the population index r can then be computed from the slope a of the log survival function S (i.e., complementary cumulative distribution function) of the gamma distribution at the reference point x_ref as: $$\begin{array}{l}{a}_s={10}^{-\frac{\text{d}}{\text{d}x}\log S_s(x_{s,\text{ref}})}\\ s=1+a_s\\ a_r={10}^{-\frac{\text{d}}{\text{d}x}\log S_r(x_{r,\text{ref}})}\\ r={10}^{a_r},\end{array}$$

where S_s and x_s,ref are fit to logarithms of mass, and S_r and x_r,ref are fit to magnitudes. Given a generic cumulative density function F(x), the survival function S(x) is simply: $$S(x)=1-F(x),$$(13)

and we compute the derivative in a numerically. The intercept b in Eq. (3) can be computed as: $$b=\log S(x_{\text{ref}})+a{x}_{\text{ref}}.$$(14)

The Kolmogorov-Smirnov test is used to compute the p-value which we use as a measure of goodness of fit.

Next, we compute the effective meteor limiting magnitude or mass l_m used in flux estimation as a point where the line f(x) = 10^ax+b is equal to the normalized cumulative value of 1, i.e., F(x) = 1: $$l_m=-b/a.$$(15)

We find that this approach is equivalent to the method of Blaauw et al. (2016) who take the mode of the fitted Gumbel distribution as the effective limiting magnitude. We suggest that estimating the limiting meteor magnitude directly from meteor data is preferable to estimating it using the stellar limiting magnitude and then applying angular velocity corrections, as our approach naturally includes all observational and instrumental biases.

Finally, the uncertainty of measured slopes is estimated using the bootstrap method (Ivezić et al. 2014, p. 140). The input data set is resampled with replacement, with the sample size being equal to the size of the input data set, and the fit is repeated on the sample. The procedure is repeated 1000 times and population or mass indices are computed for each run. Slopes that are outside 3 sigma of the mean are rejected, and the resulting standard deviation is reported as the uncertainty.

Fig. 4 Comparison of the cumulative density function for a power law and a gamma function.

Fig. 5 Distribution of peak R-magnitudes of sporadic meteors from the CAMO influx camera (histogram) and the gamma distribution fit (blue line is the main fit, black lines are all Monte Carlo runs). The jaggedness of the histogram is due to the limited precision of the magnitudes in the data file. The yellow shaded area around the inflection point represents the uncertainty in the inflection point location, and the green shaded area the uncertainty around the reference point.

3.2 Testing the method

Figure 5 shows the distribution of peak magnitudes of 2604 manually reduced double station sporadic meteors captured by the Canadian Automated Meteor Observatory’s (CAMO) influx camera (Weryk et al. 2013). The data set covers a period between June 2009 and August 2014, and only those meteors which had complete light curves (i.e., both the beginning and the end was visible inside the field of view of at least one camera) were used in the analysis. We note that imposing the visibility criterion was essential, as including meteors with partially observed tracks produced a bimodal magnitude distribution which greatly skewed the results. Although the CAMO influx camera is not used for the analysis of the Draconids, we use this high-quality data set used in previous work (e.g., Campbell-Brown 2015) to demonstrate the robustness of the method. The system has a meteor limiting magnitude of + 5.5^M, but using our approach we find that the data set is complete only up to magnitude + 3.2^M, indicated by the inflection point, while the population index is estimated at magnitude + 2.2^M.

Figure 6 shows the cumulative distribution of peak magnitudes and the line representing the tangent to the log cumulative distribution function of the fitted gamma function. The estimated slope of the line, that is, the population index, is r = 2.55 ± 0.06, which is consistent with other studies which have measured the expected yearly average sporadic population index from video data (e.g. Molau 2015). We note that the estimated slope follows the data histogram well. If the slope was to be estimated at the inflection point, the population index would be r = 2.20 and the slope would significantly deviate from the data histogram.

Figures 7 and 8 show the distribution of masses and the cumulative histogram of masses for the same CAMO data set. The computed mass index of s = 2.18 ± 0.05 is very consistent with expected values from other authors (e.g., Blaauw et al. 2011a). We note that Pokornỳ & Brown (2016) used the same data set, but their algorithm fitted a line to the histogram on the part after the inflection point, slightly reducing the value of the mass index to s = 2.09.

Due to the large number of data points, the estimate of the mass index of sporadic meteors can be accurately estimated, but often in practice, only tens of shower meteors are available. To demonstrate the robustness of the method, we estimate the mass index from 38 Geminids (all observed in 2012 within 0.5° of solar longitude around the peak) from the CAMO influx data set (Figs. 9 and 10). The fit is much more uncertain, but the estimated value of the mass index for the Geminids is s = 1.70 ± 0.14, which is consistent with previous work (Jones & Morton 1982; Arlt & Rendtel 2006; Blaauw et al. 2011b).

In the examples above, the p-values derived using the Kolmogorov-Smirnov test are larger than 0.05, which means that the null hypothesis that the model fits well to the data cannot be rejected with a 95% confidence level in any of the cases, even for the Geminids with a small number of meteors. This suggests that our theoretical assumptionsare valid.

We investigated how the population and the mass index of sporadic meteors change throughout the year using the new method. The results are shown in Figs. 11 and 12. The data were divided into ten bins of equal number of meteors (around 260 in each bin) but of variable ranges of solar longitudes. Large variations in both population and mass index can be seen at these optical meteoroid sizes, variations that are larger than the fairly stable values in radar sizes (Blaauw et al. 2011a; Pokornỳ & Brown 2016) but following the same trend. The variations and features in the population index are comparable to those derived from visual data by Rendtel (2004), except after λ_⊙ > 250° where the number of CAMO observations is very low due to bad weather conditions in winter and large bins of solar longitude had to be used. In Fig. 12, the mass index estimated using our novel method is plotted asa solid line, and the mass index computed from the population index using the classical relation s = 1+ 2.3logr (McKinley 1961) is shown as a dashed line. Comparing the two, it can be seen that the peak magnitude alone is not a good proxy of the total meteoroid mass or the mass distribution for sporadics.

Finally, we apply the method to EMCCD data collected on the night of August 13, 2018 (λ_⊙ = 140°). The meteorswere detected using a hybrid cluster–template-matching matched filter algorithm (Gural 2007, 2016), which detects even very faint meteors close to the noise floor. The method will be described in a future paper. In total, 134 double-station Perseids were observed with the two pairs of cameras during 6.5 h of observation. The meteor limiting magnitude was +5.0 for these fast Perseid meteors, and only those meteors which had complete light curves obtained by at least one camera were used in the analysis. We measured the population index for the Perseids to be r = 2.43 ± 0.29, and the mass index to be s = 1.90 ± 0.12, comparable to long-termvisual measurements at the same solar longitude (Rendtel 2014). In the same night, there were 283 sporadics whose population index was r = 3.56 ± 0.50 and mass index s = 2.30 ± 0.14, significantly higher than that of the Persieds, as expected.

Fig. 6 Cumulative distribution of peak magnitudes of sporadic meteors captured by the CAMO influx camera. The red line indicates the slope of the fitted population index.

Fig. 7 Distribution of masses of sporadic meteors from the CAMO influx camera (histogram) and the gamma distribution fit (blue line). Masses were computed using a fixed dimensionless luminous efficiency τ = 0.7%.

Fig. 8 Cumulative distribution of masses of sporadic meteors from the CAMO influx camera. The red line indicates the slope of the fitted mass index, while the shaded black lines represent the Monte Carlo fits of the resampled set as described in the text.

Fig. 9 Distribution of masses of 38 Geminid meteors from the CAMO influx camera and the gamma distribution fit. See Fig. 5 caption for the explanation of plot details.

Fig. 10 Cumulative distribution of masses of 38 Geminid meteors from the CAMO influx camera. The line indicates the slope of the fitted mass index. See Fig. 5 caption for the explanation of plot details.

Fig. 11 Annual variation in the population index for sporadic meteors using data from the CAMO influx camera.

Fig. 12 Annual variation in the sporadic mass index (solid line) as derived using the CAMO influx camera. The dashed line is the mass index derived from the population index using the classical relation s = 1 + 2.3logr.

4 Mass and population index of the 2018 Draconids

We apply the new population and mass index estimation method to the Draconids observed on the night between October 8 and 9, 2018, using the EMCCD cameras. Meteor detection and calibration was done using RMS software (Vida et al. 2016). Instead of computing photometric masses, we computed the integrated luminous intensity and used it as a proxy for mass.

Following Vida et al. (2018), the meteoroid mass m can be computed as: $$\begin{array}{lll}\hfill & \hfill & I(t)=P_{0m}{10}^{-0.4M_A(t)}\hfill \\ \hfill & \hfill & m=\frac{2}{v^2\tau }{\displaystyle {\int }_0^{t}I}(t),\hfill \end{array}$$

where I is the luminous intensity, P_0M is the power of a zero-magnitude meteor, M_A the range-corrected Gaia G-magnitude (normalized to 100 km), τ the luminous efficiency, and v the velocity of the meteor. P_0M and the term before the integral simply scale the mass, and as we are only interested in the slope of cumulative logarithms of mass, the scaling has no influence on the slope if we assume that all meteors have similar velocities. This is exactly the case, as all meteors are members of the same shower. Thus, we drop these terms and compute the dimensionless integrated luminous intensity I^* as a proxy for mass: $$I^{*}={\displaystyle {\int }_0^{t}1}{0}^{-0.4M_A(t)}.$$(18)

The two pairs of EMCCDs recorded a total of 68 double-station Draconids with complete light curves, ranging in peak magnitudes from + 1^M to + 6.6^M. We divided the data into two time bins: Bin (1) is a ~30 min bin from 00:03:33 UTC (beginning of observation) to 00:32:01 UTC that had a total of 30 Draconids and captures a part of the predicted peak; and Bin (2) is a ~60 min bin from 00:32:01 UTC to 01:28:58 UTC (end of observation) that had a total of 38 Draconids and captures the post-peak declining activity.

Fig. 13 Cumulative histogram of magnitudes using the Draconid data in Bin 1.

4.1 Bin 1

The first bin close to the peak, Bin 1, produced a population index of 1.98± 0.50 with an estimated completeness magnitude of + 4.30^M ± 0.97. Figure 13 shows that the observations roughly follow a power-law until the reference point at + 3.30^M ± 0.97, after which there is a rapid drop-off in the number of observed meteors. The model fit in this case does not seem to be robust, as indicated by the uncertainty, despite the large p-value of 0.843. The population index might have been underestimated due to the overabundance of meteors with peak magnitudes of around 4^M, possibly because of small-number statistics. Alternatively, the population index might have been increasing rapidly, as evidenced by the measured value in the second bin; thus the measured value is only the average of the population index function. Because this would break the underlying assumption of a fixed power-law, this population index, the associated effective meteor limiting magnitude of + 5.42^M ± 0.64, and the flux derived from it are to be taken with caution.

Figure 14 shows the cumulative histogram of dimensionless integrated luminous intensities and the mass index estimate of 1.74± 0.18 with the value of integrated intensity completeness being − 2.40 ± 0.39. The fit appears to be robust, and we note that the mass index estimate has a much lower uncertainty than the population index estimate. This value of the mass compares well with the value for the 2011 outburst at the peak (Koten et al. 2014).

4.2 Bin 2

For the data in the second bin, Bin 2, we measure a much larger and more uncertain population index of 3.36± 1.31 (Fig. 15), and a mass index of 2.32 ± 0.26 (Fig. 16). The large difference from the values estimated for Bin 1 indicates a fast change in the particle size distribution as the Earth moves through the stream. This rapid change is the probable cause of the increased uncertainty, and this 60-min bin likely encompasses a range of meteoroid size distributions. We note that the speed of the mass index increase is similar to what was observed for the 2011 Draconids, when the mass index changed from 1.84 to 2.30 in one hour (Koten et al. 2014).

Fig. 14 Cumulative histogram of integrated luminous intensities using the Draconid data in Bin 1.

Fig. 15 Cumulative histogram of magnitudes using the Draconid data in Bin 2.

Fig. 16 Cumulative histogram of integrated luminous intensities using the data in Bin 2.

5 Flux

Collecting areas for the two pairs of EMCCD cameras were calculated using the method of Campbell-Brown et al. (2016). A grid of points spaced by 4 km was placed at heights of between 92 and 98 km to capture the typical heights of Draconid meteors. At each height, the bias for each point in the grid was calculated, taking into account the radial speed of Draconids at the time of the observation, the sensitivities of the two cameras, and the range to each camera. The area of the 4 × 4 km region was scaled to account for meteors not seen due to lower sensitivity using the measured mass index of the stream. An average of the areas at the heights of interest was taken, though the collecting area of both systems was relatively constant at those heights.

The flux was calculated by dividing the number of two-station Draconids by the overlapping collecting area of eachcamera pair. The same two time intervals used in the mass index calculation were used for this calculation, and the flux of the first (half-hour) bin multiplied by two to get the number of meteoroids per hour. This calculation was performed for each pair of cameras separately, and then the numbers and collecting areas were added and the flux from both sets of cameras was calculated. The flux was then corrected using the s value from the system limiting magnitudes of +5.4 and +4.5 to a limiting magnitude of +6.5, which is standard when, for example, calculating ZHRs. All of these numbers, including the ZHR, are given in Table 1.

It is obvious that the collecting area of these systems is a strong function of the mass index: if the population index is higher (more faint meteors) the effective collecting area is smaller because many more meteors are expected to be missed in the less sensitive parts of the collecting area. We expect the flux in the second time bin to be lower than in the first because the peak occurred a short time before observations began, but the flux is nearly an order of magnitude larger, and the ZHR slightly higher. We note that both camera pairs were pointed close to the apparent radiant, which may skew flux calculations.

To get a sense of the uncertainty in the flux, which comes from the uncertainty in the mass index, for the first time bin we recomputed the collecting area with the mass index varying the uncertainty limits from 1.56 to 1.92. At the limiting magnitude of the system, this changed the flux from 0.047 meteoroids km⁻² h⁻¹ (at the lowest mass index) to 0.061 (at the highest); extrapolating to a limiting magnitude of +6.5 produced values from 0.081 to 0.15 meteoroids km⁻² h⁻¹. The uncertainty in the ZHR is significantly higher, varying from 389 to 1150 for the two different mass indices. There is some uncertainty in the limiting magnitude as well, which will contribute further uncertainty to these values. All of the ZHR values are higher than those observed visually by the IMO⁶, which were of order 100. Because of the rapidly changing mass index, the flux values here should be treated with caution.

6 Radiant distribution

Trajectories and orbits were computed with the Monte Carlo method of Vida et al. (2020a), and the initial velocities were estimated using the sliding fit proposed in this latter paper using 40% or more points from the beginning of the meteor. We did not apply the deceleration correction proposed by Vida et al. (2018) as accurate velocities were not the focus of this work, and this should not influence the radiant dispersion.

Figure 17 shows a probabilistic radiant map in Sun-centered geocentric ecliptic coordinates generated by combining trajectory solutions of all Monte Carlo runs (200 per meteor) using kernel density estimation and a Gaussian kernel. Figure 18 shows the nominal radiants with error bars in equatorial coordinates with color-coded solar longitude. Both plots include a comparison between the observed and simulated radiants as produced by Egal et al. (2018). In the simulation, the meteors were produced by meteoroids ejected in 1953, and we only selected simulations which hit the Earth in the range of solar longitude between 195 ° and 196 °.

The observed radiant is quite compact and the bulk of the shower falls within a circle of 0.6° in diameter (dashed circle in the plot), centered around a median shower radiant α_g = 261.997° and δ_g = 56.007°, or λ_g − λ_⊙ = 50.985° and β_g = 78.786°. Only radiants with high measurement uncertainties lie outside of this circle. The probabilistic radiant map shows an area of higher density at α_g = 262.140° and δ_g = 56.095°, or λ_g − λ_⊙ = 51.2° and β_g = 78.9°, which is marked with a smaller circle in Fig. 17, but we cannot exclude this being an artifact due to small-number statistics (only five concentrated radiants with small uncertainty). The observed radiants match the width of simulated radiants well, but there is an offset in solar longitude.

Following Kresák & Porubčan (1970), we compute the dispersion as the median angular offset from the mean radiant. Because the observation period was short, we did not perform the radiant drift correction. Figure 19 shows a histogram of angular offsets from the mean radiant. The median angular offset was 0.25 °.

In Fig. 18, the Egal et al. (2018) simulations show an obvious correlation of the radiant position with the solar longitude. Our observing period spanned the solar longitudes from 195.400696° to 195.456238°, but the simulated radiant positions matched the observations best for λ_⊙~ 195.8°, indicating that the position of simulated particles is shifted along the orbit when compared to observations. Due to the high precision of these observations, it may be possible to optimize the meteoroid ejection and integration parameters until a good match in both the position and time is achieved, improving future predictions.

Figure 20 shows histograms of simulated and measured geocentric velocities for the 2018 Draconids. The histograms were normalized so that their area is equal to 1. The mean and standard deviation were 20.05 ± 0.93 km s⁻¹ and 20.94 ± 0.08 km s⁻¹ for observedand simulated data, respectively. The observed geocentric velocities have a larger scatter than the simulations due to the limited velocity measurement precision. Furthermore, there is a systematic ~0.9 km s⁻¹ underestimation of the geocentric velocity, a combined effect of meteor deceleration prior to detection (Vida et al. 2018) and using the average velocity of the first 40% (or more) of the meteor as the initial velocity. In Table B.1 we give a list of all meteors and their geocentric radiants, orbits, magnitudes, and mass proxies. We note that due to the underestimation of the geocentric velocity from lack of deceleration correction, the semi-major axes listed in the table are also underestimated.

Fig. 17 Sun-centered geocentric ecliptic radiant distribution of the 2018 Draconids observed with the EMCCD systems. The normalized density is color coded.

Fig. 18 Geocentric radiant distribution of the 2018 Draconids observed with the EMCCD systems. The solar longitude is color coded, and the error bars show the one standard deviation uncertainty. The squares show the simulated radiants by Egal et al. (2018).

Fig. 19 Distribution of angular offsets from the mean radiant of the 2018 Draconids. The vertical line marks the median angular offset of 0.25°.

Fig. 20 Comparison of simulated (blue) and observed (red) geocentric velocities.

7 Conclusions

We analyzed 90 min of optical meteor data of the 2018 Draconid meteor shower outburst. The observations started shortly after the predicted peak at 00:00 UTC on October 9, 2018. We performed multi-station observations with the SOMN all-sky network and four sensitive EMCCD cameras. The all-sky system did not record any bright Draconid meteors, but 68 multi-station Draconids with complete light curves were recorded with the EMCCD systems. The astrometric and photometric calibration of the EMCCD systems has been discussed in detail, and we show that the EMCCD systems have a limiting stellar magnitude of + 10.0^M and an astrometric precision of 10 arcsec.

We developed a novel method of population and mass index computation which is based on maximum likelihood estimation. The methodfits a gamma function to the distribution of magnitudes or logarithms of mass to estimate the slope of the cumulative distribution (i.e., population or mass index) before the sensitivity effects become important. We validated this approach on past data obtained with the CAMO influx system. The method also enables robust error estimation.

We divided EMCCD data of the 2018 Draconid outburst into two time bins of approximately equal number of meteors: the first bin spanned the time from 00:00 UTC to 00:30 UTC, and the second bin spanned from 00:30 UTC to 01:30 UTC, on October 9, 2018. The mass index was 1.74 ± 0.18 in the first bin and2.32 ± 0.27 in the second bin, an increase of similar magnitude and rate of change to the 2011 Draconids (Koten et al. 2014).

Becauseof the rapid rate of change of the mass index of the shower and small-number statistics, the values derived for the mass index and flux are very uncertain, particularly for the first time bin. The fluxes in the two bins were 0.11 and 0.99 meteoroids km⁻² h⁻¹, respectively, corrected to a limiting meteor magnitude of +6.5. From visual observations, we expect the flux was actually falling from the first to the second time interval; the reversed trend in this data is most likely due to the assumption of a single power law in each time bin, which does not adequately describe the data. For this reason, we do not attempt to describe the flux profile of the shower.

The measured shower radiant diameter was ~0.6° which matches simulations by Egal et al. (2018) well, but there is an offset in solar longitude of 0.4°. The measured radiant dispersion is a good match to the 2011 outburst, both the high-precision radiant measurements by Borovička et al. (2014), and the meteoroid ejection simulations by Vida et al. (2020b).

Complete input data and trajectory solutions for EMCCD meteors

We provide the input data and detailed trajectory solutions for the Draconids observed with the EMCCD systems in electronic form at the CDS.

Acknowledgements

We thank Dr. Gunther Stober for helping with the equipment during the observing campaign. Furthermore, we thank Dario Zubović, Damir Šegon, Aleksandar Merlak, Damjan Nemarnik, Vladimir Jovanović, Peter Eschman, Dr. Dan Kinglesmith, Przemek Naganski, Dmitrii Rychkov, Stanislav Korotkiy, and Milan Kalina for supporting the RMS project. This work was supported in part by the NASA Meteoroid Environment Office under cooperative agreement 80NSSC18M0046. Funding from the Natural Sciences and Engineering Research Council of Canada (Grant no. RGPIN-2016-04433) and the Canada Research Chairs Program is gratefully acknowledged.

Appendix A Light curves of selected Draconids

Fig. A.1 Example light curves of EMCCD Draconids covering the range of observed peak magnitudes.

Appendix B Table of EMCCD Draconids

Appendix C Additional table

References

All Tables

All Figures

Fig. 1 Stack of raw frames showing a + 1M Draconid recorded with camera 01G on October 9, 2018 at 00:01:48 UTC. The meteor was first detected at an apparent magnitude of + 7.0M. V -band magnitudes of a selection of stars are shown as well. The vertical gradient of the image background is intrinsic to this specific camera and is compensated for by flat field correction during data reduction. The two thinner and fainter lines left of the meteor are satellites.

In the text

	Fig. 2 Absolute Gaia G-band magnitude light curve of the meteor in Fig. 1. The error bars are ~ 0.1 magnitude. Observations from camera 01G (Tavistock) are shown in blue, and observations from camera 02G (Elginfield) are shown in red.
In the text

	Fig. 3 Photometric calibration of the 01F camera. The uncalibrated magnitude is simply − 2.5 log10Is, where Is is the background-subtracted sum of pixel intensities of all star pixels (LSP in the figure). The sensor response was linear (Gamma = 1.0). The magnitudes of stars used for calibration ranged from + 2.3M to 9.9M.
In the text

	Fig. 4 Comparison of the cumulative density function for a power law and a gamma function.
In the text

Fig. 5 Distribution of peak R-magnitudes of sporadic meteors from the CAMO influx camera (histogram) and the gamma distribution fit (blue line is the main fit, black lines are all Monte Carlo runs). The jaggedness of the histogram is due to the limited precision of the magnitudes in the data file. The yellow shaded area around the inflection point represents the uncertainty in the inflection point location, and the green shaded area the uncertainty around the reference point.

In the text

	Fig. 6 Cumulative distribution of peak magnitudes of sporadic meteors captured by the CAMO influx camera. The red line indicates the slope of the fitted population index.
In the text

	Fig. 7 Distribution of masses of sporadic meteors from the CAMO influx camera (histogram) and the gamma distribution fit (blue line). Masses were computed using a fixed dimensionless luminous efficiency τ = 0.7%.
In the text

	Fig. 8 Cumulative distribution of masses of sporadic meteors from the CAMO influx camera. The red line indicates the slope of the fitted mass index, while the shaded black lines represent the Monte Carlo fits of the resampled set as described in the text.
In the text

	Fig. 9 Distribution of masses of 38 Geminid meteors from the CAMO influx camera and the gamma distribution fit. See Fig. 5 caption for the explanation of plot details.
In the text

	Fig. 10 Cumulative distribution of masses of 38 Geminid meteors from the CAMO influx camera. The line indicates the slope of the fitted mass index. See Fig. 5 caption for the explanation of plot details.
In the text

	Fig. 11 Annual variation in the population index for sporadic meteors using data from the CAMO influx camera.
In the text

	Fig. 12 Annual variation in the sporadic mass index (solid line) as derived using the CAMO influx camera. The dashed line is the mass index derived from the population index using the classical relation s = 1 + 2.3logr.
In the text

	Fig. 13 Cumulative histogram of magnitudes using the Draconid data in Bin 1.
In the text

	Fig. 14 Cumulative histogram of integrated luminous intensities using the Draconid data in Bin 1.
In the text

	Fig. 15 Cumulative histogram of magnitudes using the Draconid data in Bin 2.
In the text

	Fig. 16 Cumulative histogram of integrated luminous intensities using the data in Bin 2.
In the text

	Fig. 17 Sun-centered geocentric ecliptic radiant distribution of the 2018 Draconids observed with the EMCCD systems. The normalized density is color coded.
In the text

	Fig. 18 Geocentric radiant distribution of the 2018 Draconids observed with the EMCCD systems. The solar longitude is color coded, and the error bars show the one standard deviation uncertainty. The squares show the simulated radiants by Egal et al. (2018).
In the text

	Fig. 19 Distribution of angular offsets from the mean radiant of the 2018 Draconids. The vertical line marks the median angular offset of 0.25°.
In the text

	Fig. 20 Comparison of simulated (blue) and observed (red) geocentric velocities.
In the text

	Fig. A.1 Example light curves of EMCCD Draconids covering the range of observed peak magnitudes.
In the text