Meteor Shower Modeling: Past and Future Draconid Outbursts

TL;DR

This paper presents an updated numerical model of the Draconid meteoroid stream, accurately reproducing past outbursts and predicting future radio outbursts, including a previously unreported event in 1999.

Contribution

The study improves existing meteoroid stream models by incorporating comet dust production evolution and validates predictions with historical observations.

Findings

Accurately reproduces historical Draconid outbursts within half an hour.

Predicts three future radio outbursts in 2019, 2025, and 2029.

Identifies a previously unreported 1999 radio outburst confirmed by radar.

Abstract

This work presents numerical simulations of meteoroid streams released by comet 21P/Giacobini-Zinner over the period 1850-2030. The initial methodology, based on Vaubaillon et al. (2005), has been updated and modified to account for the evolution of the comet's dust production along its orbit. The peak time, intensity, and duration of the shower were assessed using simulated activity profiles that are calibrated to match observations of historic Draconid outbursts. The characteristics of all the main apparitions of the shower are reproduced, with a peak time accuracy of half an hour and an intensity estimate correct to within a factor of 2 (visual showers) or 3 (radio outbursts). Our model also revealed the existence of a previously unreported strong radio outburst on October 9 1999, that has since been confirmed by archival radar measurements. The first results of the model, presented in Egal et al. (2018), provided one of the best predictions of the recent 2018 outburst. Three future radio outbursts are predicted in the next decade, in 2019, 2025 and 2029. The strongest activity is expected in 2025 when the Earth encounters the young 2012 trail. Because of the dynamical uncertainties associated with comet 21P's orbital evolution between the 1959 and 1965 apparitions, observations of the 2019 radio outburst would be particularly helpful to improve the confidence of subsequent forecasts.

Tables

Table 1. Table 1: Comet and meteoroid characteristics considered by Model I. See the text for more details.

Comet parameter	Choice	Reference
Shape	Spherical
Diameter	2 km	[47]
Density	400 kg m-3
Albedo	0.05
Active surface	20 $\%$
Active below	3.7 au	[70]
Variation index ${\gamma }_1$, ${\gamma }_2$	2.1, 7.3	This work
$K_1$	$\sim$ 100	This work
$Af{\rho }_{max}$	$\sim$1600-1800 cm	This work
Particle parameter	Choice	Reference
Shape	Spherical
Size	${10}^{-4}$ to ${10}^{-1}$ m
Density	300 kg m-3	[8]
Ejection velocity model	“CNC” model	[15]
Ejection rate	each day
N${}_{\text{particles}}$ / apparition	$\sim$ 840 000

Table 2. Table 2: Comet and meteoroid characteristics used by the MSFC model (Model II).

Comet parameter	Choice	Reference
Shape	Spherical
Radius	1.7 km	See text
Active within	2.5 au
Particle parameter	Choice	Reference
Shape	Spherical
Size ($\beta$)	$\sim {10}^{-5}$ to ${10}^{-2}$
Density	1000 kg m-3
Ejection velocity model		[33]
Cap angle	30$°$
Ejection rate	every hour
N${}_{\text{particles}}$ / apparition	$\sim$ 600 000

Table 3. Table 3: Summary of the observed (left panel) and simulated (right panel) characteristics of the historic Draconid outbursts. *Peak time estimated from the stream median location; see Section 4.3 . a : reported full width at half maximum in 2018 (Molau, personal communication)

	Observations			Simulations
				Model I				MSFC model
Date	Time	Reported	ZHR	Time	Total	ZHR	Time*	Time	ZHR
(y/m/d)	(UT)	duration		(UT)	duration		(UT)	(UT)
1933/10/09	20h15	4h30	5400 to 30 000	20h13	5h	5830	20h08	20h14	strong
1946/10/10	3h40-50	3-4h	2000 to 10 000	3h38	4h30	7960	3h39	4h05	storm
1985/10/08	9h25-50	4h30	400 to 2200	9h42	8h30	370	10h43	10h33	outburst
1998/10/08	13h10	4h	700 to 1000	13h26	4h	900	13h22	13h19	strong
1999/10/09	11h15-45	5h	1250	11h13	6h45	1380	11h30	11h57	moderate
2005/10/08	16h05	$\ge$ 3h	150	16h06	7h	50	16h28	Oct. 9, 3h52	weak
2011/10/08	20h00-15	3-4h	300-400 to 560	20h06	4h	550	20h03	19h48	outburst
2012/10/08	16h40	2h	9000	16h20	10h30	3450	15h58	16h42	weak
2018/10/08	23h15-45	3h30a	100-150	22h51	8h	90	23h45	0h30	moderate

Table 4. Table 4: 2019 Draconid forecasts performed by independent modelers. Sources: (1) This work, (2) [ 64 ] , (3) Personal comment, (4) [ 56 ] , (5) [ 97 ]

Modeler	Trail	Date (y/m/d)	hour	$L_{\odot }$	ZHR	Comment
Egal (1)	Mult.	2019/10/08	12h	194.648	200 $\divideontimes$ 3	From the ZHR profile (main prediction)
	Mult.	2019/10/08	14h35	194.753	200 $\divideontimes$ 3	Maximum of the ZHR profile
	1959	2019/10/08	13h42	194.718	-	Based on stream median location
	1966	2019/10/08	12h01	194.648	-	Based on stream median location
MSFC (2)	Mult.	2019/10/08	15h06	194.775	$\sim$ to 2018	Primary peak
	Mult.	2019/10/08	4h20	194.332	-	Secondary peak
Vaubaillon (3)	Mult.	2019/10/08	14h41	193.522	4042	Old comet ephemeris
	Mult.	2019/10/07	11h01	192.811	363	Updated comet ephemeris (recommended)
Maslov (4)	1959	2019/10/08	14h44	194.759	5-10 (visual)	Likely higher activity in radio range
Ye (5)	1979			194.2 ($\in [194.0,194.5]$)		Lower rate than in 2018

Equations

$$\left\{\begin{aligned} N_{g}=&K_{2}\frac{J\Delta t}{2A(\Phi)}A_{1}(a_{1}^{\prime},a_{2}^{\prime})\color[rgb]{0,0,0}\cdot\color[rgb]{0,0,0}\left(K_{1}\Delta r_{h}\color[rgb]{1,1,1}\frac{1}{2}\color[rgb]{0,0,0}\right.\\ -&\left.\frac{\color[rgb]{0,0,0}Af\rho_{max}\color[rgb]{0,0,0}}{\gamma\ln(10)}\left[10^{-\gamma(r_{h}+\frac{\Delta r_{h}}{2}-q)}-10^{-\gamma(r_{h}-\frac{\Delta r_{h}}{2}-q)}\right]\right)\\ \gamma=&\gamma_{1}\text{ pre-perihelion and }\gamma=\gamma_{2}\text{ post-perihelion}\end{aligned}\right.$$

$$\mathcal{F}=\mathcal{F}_{s}\color[rgb]{0,0,0}\cdot V_{r}\color[rgb]{0,0,0}$$

$$\text{ZHR}=\frac{{\mathcal{A}}_{\text{s}}\mathcal{F}}{(13.1r-16.5)(r-1.3{)}^{0.748}}$$

$$\beta =\frac{F_{rad}}{F_{grav}}\simeq \frac{5.7\times 10^{-4}}{a\rho }$$

$$\text{IP}=\frac{R_{\oplus }+h_{atmos}}{D}$$

$$\text{Whipple 1951: }{V}_{ej}$$

$$\text{Jones 1995: }{V}_{ej}$$

$$\text{Modified: }{V}_{ej}$$

$$\text{Crifo 1997: }{V}_{ej}$$

$$\{\begin{array}{rl}& Af\rho (X)=K_1+Af\rho (X_{max})\ast 10^{-\gamma |X-X_{max}|}\\ & \gamma ={\gamma }_1\text{ if }X\ge X_{max}\text{ and }\gamma ={\gamma }_2\text{ if }X\le X_{max}\\ & X=(r_h-q)\frac{t(q)-t}{|t(q)-t|}\end{array}$$

$$Z_{{\text{H}}_2\text{O}}(r_h,\theta )=Q_{{\text{H}}_2\text{O}}(r_h)\frac{\cos \theta }{\pi R_c^2}$$

$$Z_g(r_h,\theta )=K(r_h)Z_{{\text{H}}_2\text{O}}(r_h,\theta )$$

$$h(a)=\frac{N}{a^u}$$

$$Z_g(r_h,\theta ,a)=K(r_h)Q_{{\text{H}}_2\text{O}}(r_h)\frac{\cos \theta }{\pi R_c^2}\frac{N}{a^u}$$

$$n_g(a,\theta ,d,r_h)=\frac{Z_g(r_h,\theta ,a)}{v_g(r_h,\theta ,a)}{\left(\frac{R_c}{d}\right)}^2$$

$$v_g(a,\theta ,r_h)$$

$$\text{where mbox }W$$

$$K(r_h)=\frac{JAf\rho (r_h)}{2N{Q}_{H_{2}O}(r_h)A(\Phi )}$$

$$Z_g(r_h,\theta ,a)=\frac{JAf\rho (r_h)}{2A(\Phi )}\frac{\cos \theta }{\pi R_c^2}\frac{1}{a^u}$$

$$N_g=$$

$$=$$

$${\int }_{a_1^{\prime }}^{a_2^{\prime }}{Z}_g(r_h,\theta ,a^{\prime })\sin \theta \text{mbox}d{a}^{\prime }d{r}_{h}d\varphi d\theta$$

$$\begin{array}{rl}{N}_g& =R_c^2\Delta t{\int }_{\theta }{\int }_{\varphi }{\int }_{r_h}{\int }_{a^{\prime }}\frac{JAf\rho (r_h)}{2A(\Phi )}\frac{\cos \theta }{\pi R_c^2}\\ & \times \frac{1}{a^{\prime u}}\sin \theta \text{mbox}d{a}^{\prime }d{r}_h{d}_{\varphi }d\theta \\ & =\frac{\Delta tJ}{2\pi A(\Phi )}{\int }_{\theta }\sin \theta \cos \theta d\theta {\int }_{\varphi }d\varphi \\ & {\int }_{r_h}Af\rho (r_h)d{r}_h{\int }_{a^{\prime }}\frac{1}{a^{\prime u}}d{a}^{\prime }\end{array}$$

$$\begin{array}{rl}{N}_g=& \frac{J\Delta t}{2A(\Phi )}{A}_1(a_1^{\prime },a_2^{\prime }){\int }_{r_h}Af\rho (r_h)d{r}_h\end{array}$$

$$\begin{split}I_{Af\rho}&=\int_{r_{h}-\frac{\Delta r_{h}}{2}}^{r_{h}+\frac{\Delta r_{h}}{2}}Af\rho(r_{h})dr_{h}\\ &=K_{1}\Delta r_{h}-\frac{\color[rgb]{0,0,0}Af\rho_{max}\color[rgb]{0,0,0}}{\gamma\ln(10)}\left[10^{-\gamma(r_{h}-q)}\right]_{r_{h}-\frac{\Delta r_{h}}{2}}^{r_{h}+\frac{\Delta r_{h}}{2}}\end{split}$$

$$\begin{split}N_{g}&=\frac{J\Delta t}{2A(\Phi)}A_{1}(a_{1}^{\prime},a_{2}^{\prime})\cdot\left(K_{1}\Delta r_{h}\color[rgb]{1,1,1}\frac{1}{2}\color[rgb]{0,0,0}\right.\\ &\left.-\frac{\color[rgb]{0,0,0}Af\rho_{max}\color[rgb]{0,0,0}}{\gamma\ln(10)}\left[10^{-\gamma(r_{h}-q)}\right]_{r_{h}-\frac{\Delta r_{h}}{2}}^{r_{h}+\frac{\Delta r_{h}}{2}}\right)\end{split}$$

$$\left\{\begin{aligned} N_{g}=&K_{2}\frac{J\Delta t}{2A(\Phi)}A_{1}(a_{1}^{\prime},a_{2}^{\prime})\cdot\left(K_{1}\Delta r_{h}\color[rgb]{1,1,1}\frac{1}{2}\color[rgb]{0,0,0}\right.\\ -&\left.\frac{\color[rgb]{0,0,0}Af\rho_{max}\color[rgb]{0,0,0}}{\gamma\ln(10)}\left[10^{-\gamma(r_{h}+\frac{\Delta r_{h}}{2}-q)}-10^{-\gamma(r_{h}-\frac{\Delta r_{h}}{2}-q)}\right]\right)\\ \gamma=&\gamma_{1}\text{ pre-perihelion and }\gamma=\gamma_{2}\text{ post-perihelion}\end{aligned}\right.$$

Full text

1 Introduction

The October Draconid shower (009 DRA) is an autumnal meteor shower known to episodically appear around the 9th of October since 1926. The Draconid shower, with apparitions irregular in time and intensity, has challenged the forecasts of modelers since its discovery. Its parent body is the Jupiter-family comet 21P/Giacobini-Zinner, discovered in 1900, which has an erratic and highly perturbed orbit [55]. The Draconid annual activity is usually barely perceptible (with a rate of a few visual meteors per hour), but the stream occasionally produces strong outbursts and storms (up to ten thousand meteors per hour) that are not directly correlated with the past geometrical configuration between the Earth and the parent comet [19]. In addition, some showers were observed by naked-eye witnesses or using video devices, while other outbursts (e.g. 2012) were caused by small meteoroids only detectable by radar instruments [97].

Previous attempts to predict upcoming Draconid activity have had mixed success. Of the nine historic outbursts of the shower (1926, 1933, 1946, 1952, 1985, 1998, 2005, 2011 and 2012), five were expected (1926, 1946, 1985, 1998 and 2011) and predictions of three events correctly estimated the peak time with an accuracy better than two hours. In 1946 the maximum activity occurred close to the descending node of the comet, simplifying timing estimates. [74] accurately predicted the shower return in 1998, caused by an encounter with the meteoroid stream ejected in 1926. More recently, a successful peak time prediction concerned the 2011 outburst. Numerical simulations of meteoroid streams ejected from 21P conducted by different authors estimated a maximum activity around 20h UT and a secondary peak around 17h on October 8 2011 [93, 90]; the existence of both peaks was indeed confirmed by independent observation campaigns [e.g. 58, 87, 36].

If different models provided an accurate prediction of the 2011 apparition date, they also led to the first reasonable estimate of the outburst’s strength. The activity of a shower, usually characterized by the zenithal hourly rate (ZHR) parameter (i.e. the number of meteors a single observer would see during an hour under ideal conditions), is a quantity that is hard to foresee. To our knowledge, no real intensity predictions were performed for Draconid outbursts prior to the 1998 apparition. The first attempt of [44] to estimate the strength of the shower outburst in 1998, expected to be lower than the one occuring in 1985, underestimated the 1998 shower’s intensity by a factor of ten.

In 2011, ZHR estimates of the main peak ranged from 40-50 [56] to 600 [93, 90], with intermediate activity predictions reaching storm level [ZHR of 7000 in 80]. Numerous observations of the 2011 Draconids estimated a ZHR between 300 [38, 58] and 400-460 [87, 36]. Among the wide range of ZHR predictions for 2011 published in the literature, some were in good agreement with the observations [93, 90]. This success could have settled the question of Draconid forecasts using numerical simulations. However, a completely unexpected radio storm reaching a ZHR of 9000$\pm$1000 meteors per hour was observed one year later [97], highlighting the need to improve the intensity predictions of meteor showers.

The motivation underlying the present work was therefore to be able to predict not only the date and time of appearance of a Draconid shower, but to also provide a quantitative estimate of its intensity with relative confidence. Since the numerical modeling of meteoroid streams ejected from the parent comet is the only efficient tool to predict such a shower, we performed simulations of the {21P, Draconid} complex using the model of [89]. For the first time, we simulated ZHR profiles for each year of predicted enhanced Draconid activity and compared them with the available observations. The simulated profiles were then calibrated against observations of the timing, activity profile and strength of historic Draconid outbursts to reinforce the reliability of future forecasts. This work presents the implementation and results of our meteoroid stream modeling, applied specifically to the Draconid shower. The structure of the paper is comprised of four main parts, some of which consist of multiple sections.

In Section 2, we detail the main characteristics of the historic Draconid outbursts observed between 1926 and 2012. Because the Draconids were widely analyzed by multiple observers using different instruments, the literature has a wide variety of information relating to the shower. For each event, we tried to highlight the conclusions shared by the highest number of authors or the best documented reports.

• 2.

Sections 3, 4 and 5, together with Appendices A, B & C, describe in detail the implementation and interpretation of the new meteoroid stream model used here to reproduce or predict past and future Draconid occurrences (Model I). In section 6, we offer a succinct summary of the independently developed NASA Meteoroid Environment Office’s MSFC meteoroid stream model, used to validate Model I.

• 3.

Section 7 presents the performance of Model I in reproducing the main visual and radio outbursts observed between 1933 and 2012, in terms of years of appearance and shower time, strength, and duration. A comparison between our simulations and observations of the previously unreported 1999 radio outburst, revealed by our model, is also presented. An initial comparison between the 2018 prediction performed in [19] and preliminary visual and radio observations of this recent outburst, shared by the International Meteor Organization (IMO), is given in section 8. The simulated activity profiles presented in these sections form the main validation of the approach employed by Model I.

• 4.

Finally, Section 9 details our predictions for three potential radio outbursts expected in the next decade (2019, 2025 and 2029). The reliability of this forecast, intrinsically correlated to Model I’s accuracy, is discussed in Section 10.

2 Draconid observations

The Draconids are known to have produced two storms observed optically in 1933 and 1946, as well as more moderate shower outbursts in 1926, 1998, and 2011. Other Draconid outbursts were detected mainly by radar techniques in 1952, 1985, and 2005, in addition to a radio storm in 2012. Very low activity from the shower was reported in 1972 instead of the intense outburst/storm that was originally predicted [61, 27].

In this work, we classify the Draconid apparitions into three categories: the outbursts which were visually observed or recorded using optical instruments (“visual outbursts”), the showers mainly observed by radio/radar devices (“radio outbursts”), and weak or poorly documented showers.

2.1 Established visual outbursts

1933

The 1933 storm occurred on October 9 around 20h15 UT [94], for a total duration of about 4h30 [R. Forbes-Bentley in 69]. ZHR estimates of the shower are very disparate and vary from around 5400 [94] to 10 000 [29] and even 30 000 [69, 14].

1946

The 1946 storm was by visual, photographic, and radar techniques mainly in Europe and North America. The shower peaked around 3h40-3h50 UT on October 10, and lasted 3 to 4 hours [53, 45]. Again, ZHR estimates in the literature vary between 2000 [28] and 6800 [45] or 10 000 [29].

1998

The 1998 outburst appeared on October 8, peaking around 13h10 UT [41, 3, 92] although [88] indicated a maximum around 13h35 UT. The total duration was about 4h [92], and the maximum ZHR reached 700 to 1000 meteors per hour [41, 3, 92].

2011

Following the concordant predictions of an outburst in 2011, the Draconids were observed on October 8 by many teams with many different observational techniques (radar, video, photography, visual). The main peak occurred around 20h-20h15 UT, for a total duration of about 3 to 4 hours and a maximum ZHR estimate varying from 300-400 to 560 [86, 38, 43, 63, 58, 87, 36].

2.2 Established radio outbursts

1985

The 1985 peak occurred on October 8 between 9h25 and 9h50 [12, 82], probably around 9h35 UT [52, 79]. The observed shower duration varies from 3h [52, 57] to 4h30 [79]. An equivalent ZHR is estimated to be from around 400-500 [52, 82] to possibly as high as 2200 [57]. In agreement with the radio measurements, visual observations carried out in Japan in 1985 confirmed a ZHR higher than 500 shortly before 10h UT on October 8 [40].

2005

The unexpected 2005 radio outburst was detected on October 8 around 16h05 UT [11]. The end of the shower was missed by radar, while the beginning was not observed using visual techniques [42]; producing only a lower limit of 3h for the shower duration. The equivalent ZHR was estimated to be around 150 [11].

2012

The 2012 storm happened on October 8 at 16h40, for a peak duration of about 2h and an equivalent ZHR of 9000 [97]. [22] measured a total duration of the shower of about 3 hours, with a peak time occurring between 16h20 and 17h40 UT.

2.3 Weak or controversial shower returns

1926

The meteor activity observed by visual observers in October 1926 allowed the shower to be linked to comet 21P/Giacobini-Zinner. From 36 meteors observed between 20h20 and 23h20 (G.M.T) on October 9 1926, a ZHR of about 20 meteors per hour was estimated [18]. A Draconid fireball, observed around 22h16 (G.M.T) that night, impressed the observers by leaving a persisting train during half an hour [18, 21].

1952/1953

Draconid activity was noticed by the Jodrell Bank radar on October 9 1952, with a maximum ZHR of 170-180 around 15h 40-50 min UT [16]. The shower duration was approximately 3 hours. The apparition of Draconid meteors the following year, in 1953, is controversial; no enhanced activity was detected by the Jodrell Bank radar at any time within 12 hours of the expected maximum time [32], while other authors suggest some weak activity [30].

1972

Due to the favorable geometrical configuration between the comet and the Earth, a potentially strong outburst/storm was expected on October 8 1972. However no significant activity was observed [62]. Radar observations revealed a diffuse component of the shower, with a weak maximum activity at October 8.2$\pm$0.3 [27]. Radio observations carried in Japan might however argue for a stronger radio activity, with a peak of 84 meteors detected in a 10 minute interval around 16h10 UT on October 8 1972 [54].

Having summarized the major observational features of the Draconids over the last century, our next goal is to reproduce the years of recorded activity, shower timing, strength and activity profile from a model. We attempt to do this employing a model with the fewest number of tunable parameters which also provides good fits to the most robustly determined characteristics of the shower.

3 Model I - Simulations

3.1 21P/Giacobini-Zinner

Comet 21P/Giacobini-Zinner was observed for the first time by M. Giacobini in 1900, and identified again by E. Zinner in 1913. Since its discovery, at least 16 apparitions of the comet have been observed and numerous orbital solutions produced. 21P is a typical Jupiter-family comet (JFC), with a period of approximately 6.5 years and a current perihelion distance of 1.03 au. The comet has suffered multiple close encounters with Jupiter through its history, resulting in our simulations in a global reduction of its semi-major axis and increase of its eccentricity in 400 years. The comet’s motion has also been affected by sudden and significant variations in the nongravitational forces (NGF) induced by its outgassing, causing 21P to be classified as an “erratic” comet [78]. A particular discontinuity in the transverse NGF coefficient between 1959 and 1965 dramatically modified 21P’s orbit [101]. The gas and dust production of the comet used to peak after perihelion before this epoch, while the maximum outgassing has occurred pre-perihelion since then [77, 7].

3.1.1 Ephemeris

Several studies have tried to reproduce the comet’s dynamical evolution during the 20th century, in particular a significant discontinuity in its orbital evolution observed between 1959-1965. Despite a moderately close approach with Jupiter in 1958, the gravitational influence of the giant planet did not seem to be responsible for the orbital variations observed around 1959 [98, 99].

[77, 46] were able to explain the NGF evolution of 21P by considering the precession of the comet’s spin axis. However, Sekanina’s model yielded unrealistic values of the comet’s oblateness and rotation period [100, 46], while the Królikowska analysis required many additional parameters derived from observations (e.g. relative times between the maximum activity and the perihelion passage for different apparitions) to correctly link all the comet’s apparitions. Another independent explanation for the erratic orbital evolution of 21P could be the activation of discrete source regions at the surface of the comet, at specific locations and for a certain duration [78].

In their current state, neither of these models permits the comet’s ephemeris to be reproduced without involving a detailed knowledge of its past apparitions; for our study, we therefore decided to rely directly on the observations. The comet’s motion at each apparition is integrated from an orbital solution provided by the JPL Small Body Data Center111https://ssd.jpl.nasa.gov/sbdb.cgi, with an external time step of 1 day and considering the asymmetric non-gravitational forces model of [101].

3.1.2 Integration period

Because of the 1959 orbital evolution discontinuity, we chose to build the ephemeris of each apparition of the comet from the closest available orbital solution, and not from the most recent measured orbit of 21P. However, older observations of 21P (prior to 1966) were determined with less accuracy, which may reduce the reliability of our integrations. In order to evaluate the impact of this approach on the comet ephemeris, we compared the trajectories integrated from 15 distinct orbital solutions of 21P obtained for the years 1900, 1913, 1926, 1933, 1940, 1946, 1959, 1966, 1972, 1979, 1985, 1992, 1998, 2006, and 2013.

Figure 1 illustrates the evolution with time of the minimum, maximum, and average similarity criterion $D_{SH}$ [84] for each pair of trajectories obtained (left panel). The middle panel presents the variations of the descending node’s heliocentric distance for each initial orbital solution considered. The right panel was obtained by considering the same initial orbit (2006 solution) but integrating the trajectory with different nongravitational force models (no NGF, symmetric outgassing, and asymmetric outgassing, respectively).

In Figure 1, we observe a slight increase of the $D_{SH}$ around 1960. The $D_{SH}$ evolution before this date coincides with the dispersion noticeable in the node distances plots. The first significant dispersion of the orbits occurs in 1898, as a consequence of a close encounter with Jupiter. Earlier close approaches with Jupiter in 1815 and 1732 also disperse the orbits. From this analysis, we conclude that the orbital solution selected for the comet ephemeris computation should not influence the prediction of the meteor showers caused by trails ejected after the comet’s discovery. For older streams, however, the initial integration epoch has a significant impact on the comet ephemeris, and hence on the meteoroid stream generated. Since the Draconids are often associated with young material [52, 96, 90], our integrations were limited to the 1850-2030 period. If any interpretation of the contribution of trails ejected prior to 1898 needs to be carefully conducted, their influence on our predictions is negligible. Particles ejected over the 1850-1898 period were only simulated for comparison of the nodal footprints of the simulated Draconids with previous works [e.g. with 90].

3.1.3 Physical properties

The lack of direct imaging of 21P’s nucleus forces us to estimate certain physical properties needed for our integrations. The radius of the nucleus of 21P is not well determined; estimates suggest it probably lies between 1 and 2 km [51, 83, 47, 70]. By default, we consider a nucleus density of 400 kg m*-3* and an albedo of 0.05 commonly assumed for JFCs [67, 48, 25]. The percentage of active surface, unknown for 21P, is fixed to 20%. A summary of the model’s parameters is presented in Table 1.

3.2 Meteoroid streams

3.2.1 Ejection

Particles are ejected at each simulation time step of the comet for heliocentric distances below 3.7 au [70]. Around 12.48 million dust particles were simulated over the period 1850-2020, covering the size bins $[10^{-4},10^{-3}]$ m ($10^{-9},10^{-6}$) kg: (160 000 particles per apparition), $[10^{-3},10^{-2}]$ m ($10^{-6},10^{-3}$) kg: (190 000 particles) and $[10^{-2},10^{-1}]$ m ($10^{-3},1$) kg: (130 000 particles). An additional sample of 120 000 particles in each size bin was ejected at each comet apparition between 1852 and 2025, and integrated until the year 2030 (9.72 million particles). In total, 22.2 million dust particles were simulated over the period 1850-2030.

The meteoroid density is assumed to be 300 kg m*-3*, as determined by Draconid meteor observations [8]. Simulated particles are isotropically ejected from the sunlit hemisphere of the comet, with velocities determined by the [15] model, which produces ejection speeds most in accord with recent in situ cometary measurements (c.f. A).

3.2.2 Integration

The stream integration is performed in Fortran 90 using a 15${}^{\text{th}}$ order RADAU integrator with a precision control parameter LL of 12 [20]. The external integration step is fixed to one day and the internal time steps are variable, allowing close encounters with planets to be dealt accurately. The gravitational attraction of the Sun, the Moon, and the eight planets of the solar system as well as general relativistic corrections are taken into account. Solar radiation pressure and Poynting-Robertson drag are also included. The Yarkovsky-Radzievskii effect was neglected since the size of our simulated particles doesn’t exceed 10 cm [91].

3.2.3 Impact selection

As a first selection, all particles that approach the Earth below a distance $\Delta X=V_{r}\Delta T$ are considered potential impactors, with $V_{r}$ the relative velocity between the planet and the particle and $\Delta T$ a time parameter depending on the shower duration. The length of a typical Draconid shower doesn’t exceed a few hours; we consider for this study a conservative time interval of $\Delta T=1$ day and a distance threshold $\Delta X$ of $1.15\times 10^{-2}$ au. Potential impacting particles are, in a second step, individually integrated with a time step of one minute, in order to precisely estimate the date and position of their closest approach with the Earth.

4 Model I - Interpretation

The primary goal of this work is to compare our simulated meteor showers to Draconid observations with the goal of matching the peak date, peak intensity, duration and overall activity profiles of past showers. Because of limited computational resources, the number of particles ejected in our model is not comparable to the true number of meteoroids effectively released by comet 21P/Giacobini-Zinner. It is therefore necessary to extrapolate from the comparatively finite number of simulated particles to the true number of particles in the stream.

Since the number of simulated particles that would physically impact the Earth is too small to derive a reliable flux profile, we consider each particle that passes within a sphere $\mathcal{S}$, centered on the Earth and with a given radius $R_{s}$, as contributing to the shower. The $R_{s}$ value must be smaller than the previous $\Delta X$ spatial criterion to exclude surrounding streams not producing any activity at the Earth. On the other hand, $R_{s}$ needs to be large enough to include a statistically useful sample of particles. The selection of the $R_{s}$ parameters for the Draconids simulations is discussed in Section 4.2.

Once particles are selected using the reduced distance criterion $R_{s}$, we assign to each impactor a weight $N_{g}$ which represents the scaling required to move from the small number of simulated meteoroids to the true number which would be released by the comet in the same ejection circumstances.

4.1 Weights

The weighting process implemented here relies mainly on [89], with some differences. First, the evolution of the dust production for 21P/Giacobini-Zinner with the heliocentric distance $r_{h}$ is extrapolated to each apparition from a photometric solution of the comet (cf. B). Second, the gas and dust production are not assumed to be proportional to each other; this decoupling allows the gas and dust production to evolve separately as a function of heliocentric distance [77]. The main hypotheses underlying the weighting solutions are:

The nucleus is spherical, homogeneous and composed of dust and water ice. 2. 2.

The dust production of the comet, represented by the $Af\rho$ parameter [2, cf. B], evolves with the heliocentric distance $r_{h}$. The gas-to-dust production ratio $K$ is not constant over $r_{h}$. 3. 3.

The water molecules and the dust particles are ejected from the sunlit hemisphere of the nucleus, with an intensity varying with the angle $\theta$ from the subsolar point. 4. 4.

The size distribution of the particles follows a power law of index $u$.

Details of the weighting computation are provided in C. For a given apparition of the comet with a perihelion distance $q$, the number of particles ejected in all directions in the size bin $[a_{1}^{\prime},a_{2}^{\prime}]$ during time $\Delta t$ and around $r_{h}$ is found to be:

[TABLE]

where:

$K_{2}$ is a tunable normalization coefficient constant over all the ejection conditions,

• -

$A(\Phi)$ is the nucleus albedo at phase angle $\Phi$,

• -

$Af\rho_{max}$ is the maximum $Af\rho$ obtained for a given apparition (ejection epoch),

• -

$r_{h}$ is the meteoroid heliocentric distance at the time of ejection (in au),

• -

$q$ is the perihelion distance (au),

• -

$A_{1}(a^{\prime}_{1},a^{\prime}_{2})$ is shorthand for $\int_{a^{\prime}_{1}}^{a^{\prime}_{2}}\frac{da}{a^{u}}$, with $a$ the particle radius (in m),

• -

$K_{1}$ is a constant governing the $Af\rho=f(r_{h})$ evolution (fit parameter),

• -

$\gamma$ is the $Af\rho$ variation index, with $\gamma=\gamma_{1}$ for pre-perihelion locations and $\gamma=\gamma_{2}$ for post-perihelion measurements, and

• -

$J$ is a function of the ejection velocity and the minimal and maximal particle radii that can be ejected from the comet (cf. C).

4.2 Flux and activity profile

The shower spatial number density $\mathcal{F}_{s}$ is estimated from the weighted number of particles entering the sphere $\mathcal{S}$ over time. In this work, the sphere radius is fixed to $R_{s}=V_{\oplus}\delta t$, with $V_{\oplus}$ the Earth’s heliocentric velocity and $\delta t$ a time parameter corresponding to the duration of the maximum activity of the shower. For the most intense Draconid outbursts (e.g. 1946), most of the activity was contained within a one-hour interval [16, 45]. For the main Draconid outbursts observed visually (1933, 1946, 1998 and 2011), we then choose $R_{s}=V_{\oplus}\times 1$ $h$ ($R_{s}<10$ Earth radii). When the total number of unweighted particles retained for the flux computation is lower than 15, $\delta t$ is gradually increased until either the number of particles exceeds 15 or $\delta t$ is 6 hours (assumed to be the maximal duration of the shower). The flux $\mathcal{F}$ of meteoroids is obtained by:

[TABLE]

with $V_{r}$ the relative velocity between the Earth and the meteoroid stream. The zenithal hourly rate of the shower is derived from $\mathcal{F}$ using [39]:

[TABLE]

where $r$ is the measured population index of the shower and $\mathcal{A_{\text{s}}}$ an observer’s surface area of the atmosphere at the ablation altitude ($\mathcal{A_{\text{s}}}\sim 37200\;\mathrm{km}^{2}$).

4.3 Peak time estimate

The predicted peak date of the shower is determined by computing a weighted average of the simulated activity profile. The presence of additional peaks, as well as the intensity and the duration of the shower, are also derived from the ZHR evolution. In the case of simulated outbursts with a low number of particles retained for the flux computation ($\delta t>1h$), the resolution of the activity profiles may compromise the validity of this approach. In such cases, the time of the peak maximum is taken to be the date of closest approach between the Earth and the median location of the meteoroid streams [89].

5 Model I - Calibration

Before predicting any future Draconid showers, we need to determine the unknown parameters ($K_{2},u$) of the weighting function in Eq. 1 which best reproduce the past observed outbursts. The results are particularly sensitive to the value of the size distribution index $u$ of the meteoroids at ejection, since this value influences the intensity and the shape of the activity profiles as well as the peak time estimate.

The calibration is first performed on the four strong visual showers of 1933, 1946, 1998, and 2011 and then refined using the radio outbursts. The simulated ZHR profiles, derived from the meteoroid fluxes using Equation 3 and a given population index $r$, are compared with the observations (see Figure 4 later). We have chosen in the calibration process not to allow a variable population index, but instead adopt a fixed value of $r=2.6$ based on measurements of the 2011 [86, 36] and potentially the 1933 shower [71].

Our philosophy is that since any model can be improved when more free parameters are introduced, we chose to limit the number of tunable parameters to the absolute minimum required. In doing this, we decrease the match between the simulated and observed intensities of a specific shower in fixing $r$. However, this allows us to determine the weighting solution that best reproduces all the outbursts with minimum adjustable parameters, which we feel improves the robustness of the model shower prediction.

The best agreement between the simulated and the observed activity profiles for the four visual showers was obtained for an ejecta size distribution index $u$ of 2.9. This value differs from our previous published estimate of 2.64 [19]; differences between the two works are due to several factors. First, in this work, we have increased the number of simulated particles by 75% and used a slightly different weighting solution which better reflects the observed $Af\rho$ evolution of comet 21P. Second, the size distribution index used here was determined considering two additional outbursts (1999 and 2018), not included in [19].

Our calibrated value of $u$ lies in the range of measurements performed for 67P [drifting from 2 beyond 2 au to 3.7 at perihelion, cf. 23]. The associated cumulative mass index ($s_{mc}=0.63$ for $u=2.9$) is somewhat lower than the value of 0.85 derived for comets 1P/Halley [59] and 81P/Wild 2 [24]. Unfortunately, no such in situ measurement have been performed for 21P. Indirect estimates of $u$, inferred by meteor observations, are hampered by the variability of the observed mass index changes with each apparition of the shower and even during an individual outburst [43]. In addition, because the orbital evolution of a meteoroid is size-dependent, a fundamental difference between the size distribution index at the time of ejection and observation is not unexpected. Since no in situ or indirect measurement of this parameter has been performed for 21P, the value of $u=2.9$ will be adopted through the rest of this study.

6 Model II - MSFC

For validation purposes, the Draconid streams simulated using Model I have been compared to an independent meteoroid modeling performed by the NASA Meteoroid Environment Office (MEO). The MEO’s MSFC Meteoroid Stream Model, similar to the model presented above, is detailed in [65] and [66]. The position and velocity of 21P/Giacobini-Zinner is taken from the 2006 orbital solution and ephemeris provided by JPL HORIZONS222https://ssd.jpl.nasa.gov/horizons.cgi. The comet radius is assumed to be 1.7 km, i.e. an average of [48, 13] and [68] measurements. 600 000 meteoroids are ejected at each apparition of the comet between 1594 and 2018, resulting in a total of about 35.4 million particles. Particles are released from the comet with a time step of 1 hour, using the ejection model of [33] with a spherical cap angle of ejection of 30$\degree$. The size of the particles are distributed uniformly over $\log\beta$, where $\beta$ ranges from $10^{-5}$ to $10^{-2}$. $\beta$ is defined as the ratio between the forces due to the radiation pressure ($F_{rad}$) and the Sun’s gravity ($F_{grav}$), and can be approximated by the following relation:

[TABLE]

with $a$ and $\rho$ respectively the radius (m) and density (kg m*-3*) of the particle. The selected $\beta$ range of $10^{-5}$ to $10^{-2}$ corresponds to masses between 1 $\mu$g and 1 kg for the selected density of 1000 kg m*-3*, and masses around ten times higher for the density of 300 kg m*-3* considered by Model I.

The meteoroid integration is performed exactly as described in Section 3.2.2. Particles are retained as potential impactors if they cross the Earth’s orbital plane within 0.01 au of the planet, and within $\pm$ 7 days of the expected shower peak. The time and intensity of the shower’s maximum activity are evaluated through computation of an impact parameter (IP), defined as

[TABLE]

where $R_{\oplus}$ is the Earth’s radius, $h_{atmos}$ the height of the atmosphere, and $D$ the Earth-particle distance at the time of nodal crossing. The impact parameter increases the weight of particles passing by the Earth more closely, and therefore more susceptible to contribute to a shower. A Lorentzian fit of the IP distribution as a function of the solar longitude allows the time of the shower peak to be estimated. The shower strength is usually estimated by scaling the IP distributions to historical observations.

Since less time and effort was allocated for the analysis of the MSFC simulations, the output of Model I and Model II cannot be compared for all the shower’s characteristics. Results of the MSFC model, successful in predicting previous Draconid outbursts (e.g. in 2011), are however presented in this work because they allow an assessment of the confidence level of the time and intensity predictions issued from Model I.

7 Outbursts post-prediction

The goal of this section is to investigate the agreement of our simulated activity profiles with historical observations of the Draconids in terms of peak time, peak intensity, and shower duration. In Section 7.1, the years of potential apparitions of the shower are appraised from the simulations. Details of each visual and radio outburst mentioned in Section 2 are detailed in Section 7.2.

7.1 Validation

7.1.1 Annual activity profile

Years of activity caused by particles of the three size bins simulated by model I are presented in Figure 2. In this diagram, all the particles meeting the large distance criterion $\Delta X=1.15\times 10^{-2}$ au for Model I and below 0.01 au for Model II are included. The annual variation of the number of impactors presented in Figure 2 is therefore not a good proxy for the shower’s strength, but instead highlights the periods of potential activity.

From our simulations, we predict intense Draconid activity in 1933, 1946, 1998, and 2011 for particles with a size larger than $10^{-3}$ m (i.e., those visible to the naked eye). We also predict activity in 1940, which to our knowledge was not reported by any observer, as well as some in 1959, 2018, and around 1926. When considering smaller particles in the size bin $[10^{-4},10^{-3}]$ m, we see additional significant activity in 1953, 1972, 1985, 1999, 2012, and 2019, as well as minor contributions over the period 1985-2030 (including the 2005 outburst).

Both Model I and Model II were successful in reproducing all the years for which a Draconid shower was observed, except the 1952 radio outburst detected only by the Jodrell Bank radar [32]. For Model I, all the visual showers are clearly distinguishable in the profile while the faintest radio events are more complicated to identify (e.g. 2005). Figure 2, shows that the models are good at reproducing the years of potential enhanced Draconid activity.

7.1.2 1999 outburst

Our model predicted strong activity rich in smaller particles on October 9 1999, which would be detectable using radio instruments. This model outcome was robust and persisted despite several attempts (not detailed here) to change the weighting parameters and distance criteria to remove this seemingly non-shower outburst year while still allowing fits to other years.

Few observations of Draconid activity in 1999 are recorded in the literature and no evidence for a strong radio shower was reported in the literature at the time. Some visual observations conducted from Europe on October 8 indicated low-level Draconid activity [49]. Visual observations performed in Japan on October 9 also indicated a weak Draconid shower, with “a broad peak from 10h00m to 13h00m UT” and meteor rates small compared with those in 1998 [ZHR = 20-30, cf. Iiyama in 75]. These observations are consistent with the simulations, which predict intense activity of particles detectable by radar, but too small to be seen optically.

In 1999, [9] conducted a Leonid observation campaign in the Canadian Arctic (at Canadian Forces Station Alert, Nunavut) involving an automated meteor radar. The radar, an early version of the Canadian Meteor Orbit radar (CMOR), consisted of the 29 and 38 MHz systems, but lacked the outlying receivers needed for orbital calculations. It was located at 82.455∘N, 62.497∘W. Only the 29 MHz system was in operation on October 9, 1999 as part of early calibration prior to the Leonids; its transmitter power was 3.46 kW. Data from this calibration interval had been collected but not examined in detail prior to the current study, at which time a strong outburst from the Draconids was found in the calibration data.

The radar echo data were processed to produce fluxes with the same single-station procedure as in [11], in which echoes occurring at 90∘ to the radiant are counted and background subtractions are applied. The collecting area of the radar was near maximum at the time of the observed shower, ranging from 500 to 750 km2. The Draconid radiant was close to its lowest elevation of 50∘. Approximately four hundred echoes occurred on the Draconid echo line in the six hours surrounding the peak, compared to 30 sporadics the previous day. Figure 3 presents the ZHR profile of the shower with time bins of 15 minutes. The peak activity occurred at 11h30$\pm$15m, for an approximate ZHR of 1250 [assuming a mass index of 1.82, c.f. 72]. Unfortunately, no other radar dedicated to meteor observation and recording in October 1999 has been identified.

Since our simulations allowed us to correctly reproduce the years of noticeable Draconid activity and provided a successful “prediction” of the 1999 radio outburst, we consider at this stage Model I to be validated. Next we investigate each simulated shower in more depth.

7.2 Activity profiles

ZHR profiles derived for each year of activity are presented in Figure 4. The first column of the figure groups the four main visual outbursts (1933, 1946, 1998, and 2011) and the second column the radio showers (1985, 1999, 2005, and 2012), whose characteristics are detailed below. The third column displays our simulations of the recent 2018 Draconids and of three potential future outbursts (in 2019, 2025 and 2029). These four cases are detailed in Sections 8 and 9.

When available, observed profiles (solid lines) have been superimposed on the simulation results (filled boxes). The maximum ZHR (number of meteors per hour) is compared here to the profile’s maximum (using a time bin of 1 hour for the model) while the peak time is estimated considering model time bins of 20 minutes. All the time estimates presented in this work have an accuracy of about half an hour. A summary of the main results shown in Figure 4 is provided in Table 3.

7.2.1 Established visual outbursts

1933

From the modelled activity profile, we predict a maximum on October 9, around 20h13 ($L_{\odot}=197.005\degree$), reaching a ZHR of 5830. The outburst is caused by particles ejected during the 1900 and 1907 apparitions, with a minor contribution from the 1886 and 1894 perihelion approaches. Most of the particles (71%) belong to the size bin $[10^{-3},10^{-2}]$ m, with a significant contribution of smaller particles (23% $\in[10^{-4},10^{-3}]$ m) and some large particles (5%). The simulated outburst duration is 5 hours, consistent with the observations.

1946

The simulated maximum occurs on October 10, around 3h38 ($L_{\odot}=196.986\degree$), reaching a ZHR of 7960. The outburst is caused by particles of several streams ejected between 1900 and 1940, with the largest contributions from the 1900 and 1907 trails. The vast majority of particles belong to the visible range, with a large number of particles having sizes between $10^{-2}$ and $10^{-1}$ m (65%). The total duration of the simulated outburst is 4h30, in good agreement with measurements. The predominance of the 1900 and 1907 trails influence on the 1933 and 1946 storms is consistent with the results of [74] and [90].

1998

The simulated activity peaks on October 8, around 13h26 ($L_{\odot}=195.087\degree$), reaching a ZHR of about 900. The outburst is caused by particles of the 1926 trail only, mainly with sizes in the ranges $[10^{-3},10^{-2}]$ m (67%) and $[10^{-4},10^{-3}]$ m (29%). The estimated shower duration is 4h. Contrary to [75], we don’t need to assume weak activity of the comet in 1926 to explain the intensity of the 1998 outburst compared to the 1933 and 1946 storms.

2011

The modelled maximum activity occurs on October 8, around 20h06 ($L_{\odot}=195.033\degree$), reaching a ZHR of 550. The outburst is mainly caused by particles from the 1900 and 1907 streams, potentially leading to two maxima around 19h30 and 20h-20h15 (cf. Section 2). However, the resolution of the simulated profile of Figure 4 does not permit reliable identification of the observed secondary peak. From our simulations, the outburst is mostly composed of particles of size $[10^{-3},10^{-2}]$ m (62%), with the remainder split evenly between smaller and larger particles. The estimated total duration of the shower is 4h, consistent with observations.

7.2.2 Established radio outbursts

1985

The simulated maximum occurs on October 8, around 9h42 ($L_{\odot}=195.260\degree$), reaching a ZHR of 370. The outburst is caused by the stream ejected in 1945, with a minor contribution from the 1933 apparition. Most of the particles contributing to the ZHR belong to the radio range ($\in[10^{-4},10^{-3}]$ m, 89%) with some candidates between $10^{-3}$ and $10^{-2}$ m (11%). No larger particles from the simulations were found to encounter the Earth. The total simulated duration (including all the activity detected and not only the central activiy plateau) is slightly shorter than 9h.

1999

The simulated activity peaks on October 9, around 11h13 ($L_{\odot}=195.728\degree$), reaching a ZHR of 1380. The outburst is exclusively caused by small particles ($\in[10^{-4},10^{-3}]$ m, 100%) ejected in 1959 and 1966. For this simulated profile, we encountered an unusual situation. A small number of particles ($<10$) belonging to the 1966 trail fulfilled all the ejection circumstances required to increase their weight (ejection at perihelion, low phase angle, very small size, etc.). This combination led to weights more than 10 times higher than any other of the 22.7 million particles simulated. Taken at face value this would produce a maximum ZHR of about 50 000 meteors per hour. Given the exceptional nature of the ejection and excessive weighting compared to all other simulated particles, they were removed from the analysis as being unphysical and from the profile of Figure 4. The total duration of the simulated outburst is about 6h45.

2005

The 2005 profile represents the worst match between our simulations and the observations. In exploring matches, we consistently found that weighting parameters adopted which better reproduce this outburst, led to a general model solution which was poorer for all other shower return fits. Allowing this poor fit between the simulations and observations for the 2005 return we feel is justified by the low intensity of the shower and the reduced confidence in the time estimates derived from the observations [11, 42].

From the coarse simulated profile, we estimate a maximum activity on October 8, around 16h06 ($L_{\odot}=195.401\degree$), reaching a ZHR of 50, composed only of particles in the radio range. Contrary to the simulations in [11], the best agreement with the observed time and intensity was found for particles ejected in 1953 and not 1946. Differences in the comet ephemeris might explain this divergence. For this poorly simulated profile, a maximum duration of 7h is expected.

2012

The simulated maximum occurs on October 8, around 16h20 ($L_{\odot}=195.610\degree$), reaching a ZHR of 3450. The outburst is only caused by small particles belonging to the 1966 stream. Similar to [97], we find that the nodal footprint of the 1966 trail does not cross the Earth’s trajectory (cf. Figure 5). This might account for our underestimate of the shower’s intensity; the 2012 outburst reached storm level with about 9000 meteors per hour. However, the particles reaching the planet’s vicinity allow us to derive an activity profile that is consistent in time with the observations, with a total duration slightly shorter than 11h.

7.2.3 Post-prediction summary

As a result of the calibration, each simulated profile of a visual outburst accurately matches the observations in terms of peak time ($<$ 30 min) and maximum intensity (within a factor of 2, consistent with the uncertainties related to the meteor observations). The duration of the simulated profiles slightly exceeds the real duration of the shower, but still allows the period of noticeably enhanced activity to be constrained.

By applying the same weighting parameters to all observed showers, the simulated profiles are less able to reproduce details of the Draconid radio outbursts. Because of the lower number of small model particles reaching the Earth ($N_{part}\simeq 15$), all the radio profiles presented in this section have irregular shapes, with gaps in solar longitude that make their comparison to radar measurements difficult. However, both peak time estimates provided in Table 3 are consistent with the observed maximum activity with an error of less than half an hour. The shower strength, which is much more variable, matched the observations within a factor of 3 in the worst cases such as 2005 and 2012. Despite the apparent discrepancies between the shower duration reported in the literature and our simulations (e.g. in 1985), a quick look at Figure 4 confirms that we correctly constrained the time window of most of the radio showers.

For the first time, to our knowledge, we have presented in this section simulated ZHR profiles that allow us to quantitatively estimate the peak time, intensity, and duration of the main Draconid outbursts with reasonable accuracy.

8 Draconids 2018

8.1 Predictions

From the simulation of 12 million particles using Model I and a slightly different weighting scheme, we performed a prediction of the 2018 Draconids in [19]. As found in several other studies, our work suggested that the Earth would cross the meteoroid stream through a gap between the 1946 and 1953 trails. The extremely small number of particles reaching the Earth’s vicinity and retained for the flux computation in this orbital configuration prevented us from confidently predicting a complete simulated profile. The peak time was therefore estimated from the date of the closest approach to the stream median location (cf. Section 4.3). This approach, which we have found to be more robust in cases with low numbers of model particles, led to a maximum activity estimate around 23h51 ($L_{\odot}=195.390\degree$) on October 8, 2018. The model maximum ZHR profile being about 80 (with a factor 3 uncertainty), we had predicted an activity not exceeding a few tens of meteors per hour. However, our predictions at the $L_{1}$ and $L_{2}$ Lagrange points suggested storm-level activity at these locations, leading to the temporary re-orientation of the Gaia satellite (personal communication from Serpell, 2018).

The updated simulated profile for 2018, calibrated using the larger sample of about 22 million particles, is presented in Figure 4. Compared to our previous work [19], we note a slight modification of the profile shape, but no significant variations of the maximum intensity (90 meteors per hour) or time (23h45 using the median location method, $L_{\odot}=195.348\degree$). In our simulations, all the particles reaching Earth belong to the 1953 trail, and have relatively small sizes ($[10^{-3},10^{-2}]$ m: 87%, $[10^{-4},10^{-3}]$ m: 12%,). Almost no large particles encountered the Earth.

8.2 Observations

At the time of writing, no definitive measurements of the Draconids 2018 have been published. However, as a preliminary comparison with our predictions, the ZHR profiles issued from the International Meteor Organization Visual Meteor Database (IMO/VMDB) and the IMO Video Meteor Network (IMO/VMN)333https://www.imo.net/draconids-outburst-on-oct-8-9/ accessed on 15/12/2018 were superimposed on the simulation results in Figure 6. The black solid curves represent the original measurements, and the light boxes the data re-binned in time steps consistent with the simulated profile (dark boxes). We observe that the length of the simulated profile matches well both sets of observations and that the peak time estimate is accurate to within a half an hour of uncertainty (depending on when the maximum activity really occurred). The simulated intensity underestimates the real strength of the shower, but remains correct within a factor of 2. From these preliminary observations, we conclude that despite the low number of particles involved in the shower simulation, our model led to one of the best predictions of the 2018 Draconids published [19].

9 Future Draconid outbursts

Finally, in an effort to predict Draconid activity over the next decade, we integrated a sample of particles generated by Model I until the year 2030. In total, 360 000 particles equally spread over our three size bins ($[10^{-4},10^{-3}]$ m, $[10^{-3},10^{-2}]$ m and $[10^{-2},10^{-1}]$ m) were ejected at each apparition of the comet from 1850 to 2030, resulting in a total of about 9.72 million simulated meteoroids. The model predicted overall peak intensity of the shower covering the period 2020-2030 is presented in Figure 2.

From this annual variation in the predicted shower strength, we estimate enhanced radio activity in 2019 and 2025, as well as a potential minor display in 2029. Comparing to the predictions of other modelers, we observe that enhanced activity in 2019 and 2025 is also predicted by [97] and [56]. We however do not share the forecast of a low intensity shower in 2021 mentioned in [97].

9.1 2019

On October 8 2019, the model predicted that the Earth will face a similar situation to that in 1999 and encounter the streams ejected by 21P in 1959 and 1966. Our estimated activity profile, presented in Figure 4, is exclusively composed of particles in the size bin $[10^{-4},10^{-3}]$ m. As with the radio outburst in 1999, a small number of unphysically heavy weighted particles belonging to the 1966 trail raised the ZHR maximum to a storm value of 5000 meteors per hour and have been removed from the analysis.

The low temporal resolution of the simulated profile, typical of our other simulated radio showers, decreases the precision of the time and intensity predictions compared to strong visual outbursts. We estimate, however, that the Earth might experience a radio outburst with a maximum ZHR of about 200 (within a factor of 3), and enhanced activity between 5h45 and 15h30 UT on October 8. The peak time is difficult to assess from the model profile. If the maximum ZHR is reached around 14h35 UT, ($L_{\odot}=194.753\degree$), the weighted average of the profile has a peak time of 12h UT ($L_{\odot}=194.648\degree$). When considering the date of the closest approach between the Earth and the median location of a specific stream, we estimate a maximum activity caused by the 1959 stream around 13h42 UT ($L_{\odot}=194.718\degree$) and a maximum caused by the 1966 stream at 12h01 UT. The total duration of the shower, among the longest for radio showers in our simulations, is 9h45.

Simulations performed by other authors differ from our predictions, either in time and intensity. Results from the MSFC model, differing from Model I mostly on the ejection circumstances and the weighting process, augurs for noticeable Draconid activity (nearly as strong as in 2018) culminating around 15h06 UT on October 8 ($L_{\odot}=194.775\degree$) consisting of mainly smaller particles. A secondary peak might be expected much earlier, around 4h20 the same day ($L_{\odot}=194.332\degree$). According to [56], the activity level from the 1959 trail in 2019 should be about half of that in 1999. Based on the 1999 ZHR estimated by [Iiyama in 75], Maslov expects “a small visual peak not higher than ZHR 5-10, on 2019 October 8 at 14h45m UT”. He however points out that “radio observations could show much higher activity” [56].

In [30], no enhanced Draconid activity was predicted by J. Vaubaillon in 2019. However, some revised predictions are presented in Table 4. The first one points toward a potential storm on October 8, with a peak time consistent with the MSFC and [56]’s results. His ZHR estimate of about 4000 meteors per hour is also in agreement with our simulations when we do not remove the heavily weighted particles. The second prediction, favored by that author because it uses the more accurate ephemeris of the comet provided by JPL, indicates a moderate outburst on October 7 around 11h UT (Vaubaillon, personal communication).

In [37], no specific activity is discernible for 2019 compared to historic Draconid outbursts. [97] attribute enhanced activity in 2019 to small particles belonging to the 1979 trail. The maximum activity is expected to be lower than in 2018 and to peak around $L_{\odot}=194.2\degree$ ($\in[194,195]\degree$).

A summary of the 2019 predictions is provided in Table 4. Different estimates of the shower’s characteristics are highly variable and no general agreement about the peak time and intensity is found. When considering the weighted average of our simulated activity profile, the peak time is driven by the influence of the 1966 trail and the resulting prediction is in advance of three hours compared to the MSFC or the [56] results. The maximum of the ZHR distribution is consistent with these two models, but was not used in our model’s validation. The expectations of the shower’s strength are even more uncertain, and vary from a low activity [56, 97, 37] to a radio outburst (this work, Model I and MSFC model).

9.2 2025

In our simulations, the Earth crosses a portion of the meteoroid stream composed of small particles at the 2025 Draconid return mainly ejected during the 2011-2012 apparition, with some contribution from the 2005 and potentially the 1999 trails. From the activity profile, we estimate a maximum activity around 16h18 UT on October 8 ($L_{\odot}=195.286$). The Earth is expected to approach the median location of the 2012 trail around 15h30 UT ($L_{\odot}=195.252\degree$) and the 2005 trail at 13h40 UT ($L_{\odot}=195.178\degree$). Our maximum intensity ($\divideontimes$ 3) reaches a ZHR of 950 meteors per hour.

According to [56], encounters with the 1907 to 1953 trails will lead to multiple submaxima at 05h01 UT (ZHR: 10-15), 07h25 UT (ZHR: 20-25), 09h06 UT (ZHR: 20-25), and 10h17 to 10h49 UT (ZHR: 50-60), with many fireballs expected. Consistent with our prediction, [56] also predicts an encounter with the 2012 trail on October 8 2025 around 15h14 UT. If the expected visual activity is low (10-40 in ZHR), the author anticipates a much higher radio activity going up to tens of thousand of meteors per hour. Because both forecasts point toward a very strong radio outburst/storm caused by an encounter with the young (and therefore probably dense) 2012 trail, special attention should be paid to this event.

9.3 2029

The 2029 ZHR profile and node plot (cf. Figures 4 and 5) are hard to interpret. Only some particles ($Npart<10$) ejected in 1940 approach the Earth in our simulations, leading to the low-resolution shape of Figure 4. For this specific shower, more simulations of the 1940 trail in particular are needed. Our first impressions are that we might detect some weak radio activity the night between October 7 and 8, which was not predicted by [97] or [56].

10 Discussion

Numerical models of meteor showers are forced to rely on multiple assumptions. Each attempt to better reproduce the comet’s behavior requires the incorporation of new parameters, which are usually not measured for the parent body in question. In this section, we address the influence that certain of our assumptions have on our forecasts of future Draconid outbursts.

10.1 Comet characteristics and ephemeris

As mentioned in Section 3.1, after one century of observations, comet 21P/Giacobini-Zinner remains a mysterious comet. Its physical characteristics are not well defined, and the accuracy of its ephemeris is poor prior to 1966 because of multiple close encounters with Jupiter and erratic variations in its nongravitational forces. A significant orbital discontinuity was observed in 1959-1965, changing the comet ejection patterns and orbital evolution. These considerations led us to adopt orbital solutions available for 21P based on observations at a particular apparition instead of computing its ephemeris backward in time from a precise and recent solution. This decision, forcing us to sometimes rely on less precise orbital measurements (but allowing us to take the 1959-1965 discontinuity into account), may explain the preponderance of different streams involved in our simulated 2005 and 2019 showers compared to other authors [11, 97]. As noticed in [90], the choice of the ephemeris solution has, therefore, an impact on the shower predictions.

By comparing the post-predictions performed using Model I and the MSFC model, which is similar in methodology but assumes different physical characteristics of the comet and ejection circumstances, we observe that both approaches provided consistent conclusions about the years of enhanced Draconid activity, the streams contributing to the shower, and the particles’ characteristics (size and ejection velocity). We also notice that the peak time estimates of the historic Draconid outbursts predicted by Model I (from the profiles and the stream median location) and by the MSFC model are similar. We therefore consider the exact choice of the comet’s physical parameters and the ejection circumstances to have only a moderate influence in the predictions of the existence and date of enhanced Draconid activity.

10.2 Weights and ZHR computation

The weighting solution for comet 21P, detailed in B and C, has an extremely significant impact on the shower’s predicted level of activity. In this work, we have updated the weighting solution of [89] considering pre- and post-perihelion dust measurements of comet 21P during the 2018 apparition. The heliocentric distance dependence of the amount of dust ejected by the comet was established from these observations and applied to all the previous apparitions considered in our model.

Despite the evident limitations of this approach, this solution offers a first approximation of the comet’s behavior throughout its period of activity and has the merit of being based on real measurements of the comet’s activity. Any further modification of 21P’s $Af\rho=f(r_{h})$ can also easily be implemented in the model to improve future shower predictions.

Among the unknown parameters of the weighting function 1, the most influential and uncertain one is the particle size distribution index at ejection, $u$. Coupled with the number of particles required to produce a ZHR profile ($>15$) and the maximum sphere size adopted for the flux computation ($6h\times V_{\oplus}$), the value of $u$ can substantially modify the shape, peak time, and intensity of the simulated profiles. It is therefore difficult to conclude that the profiles presented in Figure 4 provide an accurate prediction of future Draconid outbursts. However, by calibrating these parameters to strong and well-known visual outbursts (cf. Section 5), their impact on the final prediction is reduced.

To illustrate this, we compare the conclusions of the present work with the previous analysis published in [19]. We observe that with a different weighting solution, a different size distribution index, a higher number of particles required for the profile computation (15 instead of 10), and a larger sample of particles simulated, the accuracy of the time and intensity predictions over the period 1933-2018 are still comparable within uncertainties between the two works. The incorporation of the comet’s $Af\rho$ evolution into the model allowed [19] to obtain correct time and ZHR estimates for historic Draconid outbursts and to perform a high-fidelity prediction of the 2018 outburst (cf. Section 8). In this work, the incorporation within the weighting solution of new measurements of the comet’s activity (including post-perihelion observations) allowed an increase in the accuracy of our modeling. At the same time the calibration process maintained consistency with our previous work despite significant variations of our ZHR profile computation.

10.2.1 Potential 2019 outburst

As described in Section 9, a prediction for the 2019 Draconids was particularly difficult to establish and our results diverge from other studies. The low resolution of the simulated activity profile hampers a reliable peak time estimate, while the maximum intensity prediction is complicated by a few unphysically heavily-weighted 1966 particles that we chose to remove. Nevertheless, we predict a radio outburst on October 8 around 12h UT; most other authors expect a lower maximum activity later the same day. Radar observations of the 2019 Draconids will therefore be a good opportunity to both test and improve our model. Precise measurements would help to determine the veracity of our choice to limit the contribution of the 1966 stream, as well as the ephemeris computation implemented in our model.

The 2019 Draconid return will also offer the opportunity to study two of the most problematic trails in the Draconid complex. The 1959 and 1966 streams, both of which cross the Earth’s trajectory in 2019, previously encountered the planet in 1999 and caused a strong radio outburst that was unnoticed until this work. They also intercepted the Earth in 2012, producing a radio storm that was unpredicted by every shower modeler. The fact that 21P suffered an orbital discontinuity precisely during the period 1959-1965 might also point to the activation of local source regions over the comet’s surface [78] and a potential for much higher activity than that considered in this model. Special attention should therefore be paid to observations of the coming shower on October 8 2019.

11 Conclusion

In this work, we performed numerical simulations of comet 21P/Giacobini-Zinner and the Draconid meteoroid stream complex following the methodology of [89]. We generated a total of about 22 million dust particles, ejected by the comet since 1850. Particles were integrated until the year 2030, and candidates approaching the Earth were carefully selected and analyzed in order to identify past and future Draconid outbursts.

To our knowledge, this study is the first to produce ZHR profiles for each simulated Draconid shower and quantitatively estimate the peak time and intensity of all the reported historic outbursts. Our model provides a timing accuracy of order half an hour for the peak times and a ZHR estimate correct to within a factor of 2 (visual showers) or 3 (radio outbursts). These successful post-predictions are the result of an updated weighting scheme of the simulated meteoroids, allowing us to better assess the contribution of individual particles to the shower. The initial algorithm of [89] has been updated and modified to take into account the evolution of the parent comet dust production along its orbit, adjusted using recent measurements of 21P performed by the NASA Meteoroid Environment Office between May and December 2018. Thanks to this updated solution, we were able to uncover the existence of a previously unreported strong radio outburst in 1999 that we confirmed by archival radar measurements at the predicted date and time. A first set of simulations issued from this model also led to one of the best published predictions of the recent 2018 outburst [19].

In the next decade, we can expect up to three Draconid outbursts caused by radar-sized meteoroids, respectively in 2019, 2025 and potentially 2029. Observations of the 2019 outburst, which we find is among the most challenging of all returns to characterize from our simulations, should particularly help in improving our stream modeling for the Draconid returns of the 2020s.

From our model predictions we expect a strong radio outburst on October 8 2025 and a somewhat weaker shower in 2019. We encourage observers to be particularly vigilant in recording the Draconid shower in these years.

Acknowledgements

We are thankful to J. Vaubaillon for his support and advice regarding his model’s implementation. This work was supported in part by NASA Meteoroid Environment Office under cooperative agreement 80NSSC18M0046 and contract 80MSFC18C0011.