© A. Egal et al. 2020

1 Introduction

In recent years, thanks to modern numerical simulations of meteoroid streams, the timing of meteor outbursts has reached unprecedented accuracy (e.g., Asher 1999). In contrast, reliable estimates of meteor shower duration and intensity are less robust (Janches et al. 2020). This is due to both a lack of measurement of the parent body’s past activity and a lack of long-term, consistent records of the resulting meteor showers.

The modeling of comet 1P/Halley and Halleyid meteor showers, namely the η-Aquariids and Orionid, suffers from these observational limitations. Simulations of the meteoroids released by the comet a few millennia ago are limited by the increasing uncertainty of 1P/Halley’s past activity and dynamics. In addition, the lack of consistent observations of the Halleyids (and especially of the η-Aquariids) has hampered development and validation of new models of the shower (e.g., see Sekhar & Asher 2014).

In this work, we present the second part of our analysis of the Halleyid meteor showers. In Egal et al. (2020), we documented the long-term activity (around several decades) of the η-Aquariids and Orionids. The showers’ activity profiles, shapes, and year-to-year variations were defined through a compilation of visual meteor observations made since 1985, more recent measurements of the Canadian Meteor Orbit Radar (CMOR, cf. Brown et al. 2008, 2010), and the International Meteor Organization Video Meteor Network (IMO VMN; Molau & Rendtel 2009). That paper reported results of multi-instrumental measurements of the Halleyids over several decades using consistent analysis techniques. The long-term characterization of the Halleyids provides the observations needed to validate modeling of 1P/Halley’s meteoroid complex.

The purpose of the present work is to apply the new numerical model of meteoroid stream evolution developed in Egal et al. (2019) to the Halleyids. In particular, the long-term evolution of these streams can now be constrained by the 35 yr of observations reported in Egal et al. (2020). Through extensive numerical simulations and the application of an adapted weighting solution to the ejected particles (see, e.g., Egal et al. 2019), we constrained the model parameters using the η-Aquariid and Orionidintensity profiles and interannual variability since 1985. The simulation parameters are calibrated on the tens of annual apparitions of each shower to increase the model’s robustness. The dynamical evolution of 1P/Halley’s meteoroid streams are characterized and linked to the main characteristics of the present-day observations of the η-Aquariid and Orionidmeteors, including origin, age, and variability. After model calibration of the observed showers before 2020, we provide a forecast of the Halleyids activity up to 2050.

The paper is divided to four main parts as follows. Section 2 reviews the characteristics of 1P/Halley’s nucleus, dust production activity, and dynamics that will serve as starting conditions for our numerical model. Section 3 presents a general description of the η-Aquariid and Orionid meteor showers, summarizing the main conclusions of Egal et al. (2020). A short historical review of prior Halleyid models is provided in Sect. 4. Section 5 details our new numerical model of 1P/Halley’s meteoroid stream. This description is followed with a reconnaissance analysis of the meteoroid trails’ evolution (Sect. 6) and the predicted associated meteor activity on Earth (Sect. 7). Sections 8 and 9 present the core of our analysis. The simulated activity profiles are compared and calibrated with meteor observations, and the activity of the η-Aquariid and Orionidmeteor showers is forecast through 2050. The successes and limitations of our model are discussed in Sect. 10.

2 1P/Halley

1P/Halley is among the most famous of all comets. It is named after Sir Edmond Halley, who recognized the periodic nature of the body forthe first time in comet science and successfully predicted its return close to the Sun in 1758. 1P/Halley moves on a high-inclination retrograde orbit, with a period of revolution of about 70–80 yr. It is the only Halley-type comet in recorded history that can be observed twice in a human lifetime. The comet often appears bright in the sky, as seen from the Earth and depending on the viewing geometry, and is therefore commonly recorded in historical annals.

It has been suggested that the first observation of the comet may date to more than 2000 yr ago, based on Chinese records (Kiang 1972). In 1986, the return of 1P/Halley to perihelion generated unprecedented enthusiasm for cometary observations. The nucleus was approached by an international fleet of five spacecraft (Giotto, Vega 1, Vega 2, Suisei, and Sakigake), collectively termed the “Halley Armada.” In situ space missions were supported by additional space observations (performed by the ISEE-3, Pioneer 7, and Pioneer 12 spacecraft) as well as by ground-based telescopic observations. The observations of the comet and the associated meteoroid streams were coordinated, collected, and archived by the International Halley Watch (IHW) organization.

Most of our knowledge of 1P/Halley’s composition, morphology, and activity was obtained during the extensive observation campaigns organized during the 1986 apparition. In this section, we review the results from the IHW and contemporary studies relating to 1P/Halley that are pertinent for meteoroid stream modeling.

2.1 Nucleus

The spacecraft flybys of 1P/Halley revealed that the comet’s nucleus was asymmetric with an avocado-like shape having dimensions of 15 × 8 × 8 km (Whipple 1987) and a mean effective cross-sectional area of about 80–100 km² (Sagdeev et al. 1988; Keller et al. 1986). Estimates of the nucleus albedo were 0.04–0.05 (Whipple 1987; Keller et al. 1987; Hughes 1987a).

The density of the nucleus is not accurately known. By analyzing the influence of nongravitational forces on the motion of the comet, Rickman (1986) estimated the mass and the density of the nucleus. Using the nongravitational acceleration coefficients found by Yeomans (1977) and Kiang (1972), Rickman (1986) estimated a most probable nucleus bulk density of about 0.1 to 0.2 g cm⁻³, with an upper limit on the order of 0.3 g cm⁻³ when the asymmetric sublimation rate around perihelion was used. Applying the same technique, but adding data closer to perihelion, Sagdeev et al. (1988) estimated a density of 0.6 g cm⁻³. Work by Whipple (1987) suggested upper limits of 0.4 to 0.5 g cm⁻³.

2.2 Activity

2.2.1 Location

Ice sublimation and the associated dust ejection from 1P/Halley’s surface are generally confined to the sunlit hemisphere of the comet. Images taken by the Halley Multicolor Camera (HMC) on board the Giotto spacecraft revealed a concentration of the dust emission in a sector spanning an angle of 70 degrees of the solar direction (Keller et al. 1987). Thermal infrared observationsin March 1986 highlighted the existence of a sunward fan and the absence of a well-defined tail (Campins et al. 1987).

The Giotto camera recorded several active jets on the comet’s sunlit hemisphere, while the rest of the surface was largely inactive (Emerich et al. 1987; Sekanina 1987). The percentage of active area on the nucleus was estimated to be 10% (Keller et al. 1987)or 15% (Whipple 1987). Localized jet activity was found to be variable in time and location on the nucleus (Whipple 1987). The onset of activity was first observed at 6 AU, well beyond the limiting heliocentric distance for water ice sublimation (Wyckoff et al. 1985; Whipple 1987; Hughes 1987a).

2.2.2 Afρ

Several dust production rates were estimated during the Giotto and Vega 1-2 flybys of the nucleus (e.g., Edenhofer et al. 1987; Hanner et al. 1987; Mazets et al. 1987; Krasnopolsky et al. 1987). However, these estimates are strongly dependent on the spacecraft’s distance from the nucleus and do not necessarily reflect long-term variations in the cometary activity.

Cometary dust production is frequently characterized by the Afρ parameter defined by A’Hearn et al. (1984). The Afρ estimates depend on the comet’s Bond albedo, the phase angle, the filling factor of the dust grains within the field of view on the plane of the sky, and the aperture radius used for an observation. The Afρ measurements are independent of the aperture radius ρ when the comet’s radial brightness varies as ρ⁻¹ close to the nucleus. However, observations of 1P/Halley revealed a strong variation in Afρ as a function of aperture size, implying that the inner coma brightness of the comet decreases more steeply than ρ⁻¹ before perihelion (Schleicher et al. 1998).

Comparisons of different Afρ measurements of comet 1P/Halley therefore require a correction for the aperture radius. Examples of debiased Afρ estimates can be found in Fink (1994) and Schleicher et al. (1998). In a second step, raw measurements need to be corrected for phase angle effects.

The variation of 1P/Halley’s Afρ as a function of the phase angle θ has been determined by Schleicher et al. (1998) for θ between 0° and 60°. Schleicher et al. (1998)’s dust phase function of 1P/Halley became a reference relation to correct Afρ measurements of other comets at small phase angles (e.g., Blaauw et al. 2011). Current cometary dust phase corrections are usually performed using a composite dust phase function, combining the Schleicher et al. (1998) relation for θ ≤ 42°: $${\text{log}}_{10}Af\rho (\theta )=-0.01807\theta +0.000177{\theta }^2$$

and the Marcus (2007a,b) phase relation based on the Henyey-Greenstein function for θ > 42°.

The left-hand panel of Fig. 1 presents the original Afρ measurements of comet 1P/Halley, corrected for the aperture radius used, as presented in Fink (1994) and Schleicher et al. (1998). Phase-corrected estimates of Schleicher et al. (1998) are presented in the figure’s right-hand panel. Original estimates of Fink (1994)were decoupled from phase angle effects via the composite dust phase function detailed above and plotted along with the Schleicher et al. (1998) measurements in Fig. 1.

Following the procedure of Egal et al. (2019), we fit 1P/Halley’s Afρ profile from Fig. 1 with a double-exponential function of the form: $$\begin{array}{l}Af\rho (X)=K_1+Af\rho (X_{\text{max}})\times {10}^{-\gamma |X-X_{\text{max}}|},\\ \gamma ={\gamma }_1\text{ if }X\le X_{\text{max}}\text{ and }\gamma ={\gamma }_2\text{ if }X\ge X_{\text{max}},\\ X=(r_{\text{h}}-q)\frac{t-t(q)}{|t-t(q)|},\end{array}$$(1)

where r_h is the comet’s heliocentric distance at time t, q is the perihelion distance, and the γ_1,2, X_max, and K₁ parameters are to be determined by the fitting process. The best agreement with measurements of 1P’s 1985 apparition gives K₁ ~ 3500 cm, Afρ(X_max) ~ 105 285 cm, γ₁ = 2.97 (pre-perihelion branch), γ₂ = 0.99 (post-perihelion branch), and X_max = −0.04 AU, as illustrated by the black line in Fig. 1. The comet dust production is asymmetric, with a maximum around perihelion and a steeper pre-perihelion profile compared to post-perihelion.

Fig. 1 Raw (left panel) and phase-corrected (right) Afρ measurements of comet 1P/Halley during the 1986 apparition presented in Schleicher et al. (1998) and Fink (1994). The abscissa (X) represents the comet’s heliocentric distance to perihelion, arbitrarily negative for pre-perihelion measurements and positive for post-perihelion times ($X=(r_h-q)\frac{t-t(q)}{|t-t(q)|}$, with t(q) being the time at perihelion).

2.3 Dust

In situ measurements of the dust environment of comet 1P/Halley have provided insight into its meteoroids’ composition, ejection speed, size, and mass distributions. The Giotto spacecraft measured two distinct populations of extremely small dust grains, down to 10⁻¹⁶g in mass: one silicate-rich and the other referred to as CHON (for carbon, oxygen, hydrogen, and nitrogen) particles (Kissel et al. 1986a,b).

The dust’s bulk density was estimated to be 0.297 g cm⁻³ from Vega-2 measurements of the dust scattering (Krasnopolsky et al. 1987). From thermal infrared images, the estimated dust albedo decreased from 0.03 to 0.02 for phase angles between 20° and 50° (Sekanina 1987).

Measurements of the meteoroids’ mass distribution were performed by impact detectors on-board both the Vega and Giotto spacecraft at different locations in the coma (Mazets et al. 1986; Simpson et al. 1986; McDonnell et al. 1986; Zarnecki et al. 1986; Vaisberg et al. 1987). The mass index s was found to vary with the distance between the spacecraft and the nucleus in an irregular manner. In addition, a slope turnover in the mass distribution was identified for small masses (~ 10⁻⁸−10⁻¹³ kg), but this value showed intersensor differences as well as variations with nuclear distance (Hanner et al. 1987; McDonnell et al. 1987).

A summary of the mass distribution indices s measured by the different spacecraft is presented in McDonnell et al. (1987) and Hughes (1987b). Estimates of s (the differential mass distribution exponent in dN = Cm^−sdm, where dN isthe number of grains between m and m+dm) ranged between 1.49 and 2.03 for masses between 10⁻⁶ and 10⁻¹³ kg, with an average value of 1.9 to 2.1 (Hanner et al. 1987). Measurement of the meteoroids’ bulk mass density and the mass distribution index close to the nucleus are in general agreement with the values derived from meteor observations (see Sect. 3).

2.4 Ephemeris

A key parameter in modeling the Halleyid stream is the past orbital evolution of the parent comet, 1P/Halley. We expect that uncertainties in the osculating orbit with time will have a significant influence on our ability to reproduce the present day streams through modeling. As a result, we critically examined all available information on 1P/Halley’s previous orbital behavior with an eye toward identifying those elements of each stream that are most secure given the comet’s ephemeris uncertainty.

2.4.1 Since 1404 BCE

Observations of Halley’s comet could date back as far as 240 BCE in Chinese records (Kiang 1972) and 164 BCE in Babylonian records (Stephenson et al. 1985). Despite the existence of old observations that provide valuable constraints on the ephemeris of the comet, no single set of orbital elements successfully links all the recorded apparitions of 1P/Halley. Numerical integration of the comet’s motion back in time from recent observations starts diverging from ancient observations around 218–374 CE (Yeomans & Kiang 1981; Stephenson et al. 1985). It has been suggested that variations of the comet’s nongravitational acceleration might be responsible for the differences between its modeled and observed past motion (Sitarski & Ziolkowski 1987).

Different sets of osculating elements of the comet have been proposed in the literature (see for example Chang 1979; Brady 1982; Yeomans & Kiang 1981; Landgraf 1984, 1986; Sitarski & Ziolkowski 1987). To reproduce modern and ancient observations, Landgraf (1984) introduced a linear variation of the comet’s nongravitational parameters. Yeomans & Kiang (1981), in contrast, considered a new determination of the comet’s perihelion passage times in ancient Chinese observations to constrain the ephemeris. These authors assumed constant nongravitational coefficients with time but had to make subjective changes to the comet’s eccentricity during the Earth’s close approach of 837 CE to match older observations. Yeomans & Kiang (1981) integrated the motion of the comet back to 1404 BCE, when an additional close encounter with Earth and the lack of older observational constraints on the comet’s motion prevented extension of the analysis further back in time.

Despite the manual modification of the comet’s orbital elements in 837 CE by Yeomans & Kiang (1981), the resulting ephemeris provides the most accurate reproduction of ancient and modern observations of 1P/Halley (e.g., Stephenson et al. 1985; Sitarski & Ziolkowski 1987). Yeomans & Kiang (1981)’s solution for 1P/Halley’s ephemeris between 1404 BCE and 1910 CE is usually taken as the definitive reference for the comet’s orbital evolution for studies of its meteoroid streams (e.g., McIntosh & Jones 1988; Ryabova 2003; Sekhar & Asher 2014). The evolution of 1P/Halley’s nodal location since 1404 BCE, computed with Yeomans & Kiang (1981)’s ephemeris (cf. Fig. 1 in Egal et al. 2020), is presented in Appendix A. The appendix also reproduces the variations of the minimum orbit intersection distance (MOID) between the Earth and 1P/Halley around each node (cf. Fig. 2 in Egal et al. 2020), measuring the shortest Cartesian distance between thetwo osculating orbits when considering the pre-perihelion or the post-perihelion arcs of the comet’s orbit.

2.4.2 Long-term (pre-1404 BCE) evolution

A complete model of the Halleyids requires a thorough understanding of the comet’s orbital evolution prior to the 1404 BCE close encounter with Earth, as the streams are certainly much older. In the absence of earlier observations, the location of the comet becomes more unconstrained prior to 1404 BCE. However, the motion of 1P/Halley over longer periods of time has been investigated by several authors. As the orbital plane of the comet regresses backwards in time, the descending node of 1P/Halley approaches Jupiter’s orbit. The date of this unobserved, but probable, close approach with the giant planet can be used to constrain the dynamical lifetime of the comet.

Using this idea, Yeomans (1986) estimated the dynamical age of 1P/Halley to be at least 16 000 yr. Hajduk (1985) estimated a maximum dynamical age within the inner Solar System of 200 000 yr. The theoretical model of the Halleyids by McIntosh & Hajduk (1983) estimated the dynamical age of 1P/Halley at around 100 000 yr. Sekhar & Asher (2014)’s numerical integrations from Yeomans & Kiang (1981)’s 240 BCE orbital solution showed that close encounters with Jupiter significantly modified the comet orbit about 12 000 yr in the past.

We simulated the past orbital evolution of the comet in a similar way. We used as starting conditions a recent orbital solution provided by the Jet Propulsion Laboratory (JPL) and assumed constant nongravitational acceleration parameters with time. The comet was integrated from the JPL¹ solution until 15 000 BCE. A thousand clones of the comet’s orbit were created according to the multivariate normal distribution defined by the solution’s covariance matrix. The time evolution of the nominal (green), average (red), and standard deviation (blue) of the clones’ semi-major axes is presented in Fig. 2 (top panel). After a slight dispersion of the clone swarm around 800 BCE, the standard deviation of the semi-major axis increases almost linearly into the past. A close encounter with Jupiter around 10 500 BCE(black box) accelerated the rate of dispersion of the comet’s clones, significantly increasing the uncertainty in the comet’s ephemeris at earlier times. The orbital evolution of the other osculating elements are presented in Appendix A.

Our integration is in good agreement with Sekhar & Asher (2014)’s results but does not take into account the manual corrections of 1P/Halley’s ephemeris around 837 CE performed by Yeomans & Kiang (1981). To check the validity of our analysis, we also generated a second sample of clones from Yeomans & Kiang (1981)’s solution of the comet. The motion of 1P/Halley was integrated starting from 141 CE, at which time it matched the comet’s observed perihelion passage with an accuracy of 0.08 days. Since no measurement uncertainties are available for the Yeomans & Kiang (1981) solution, the comet’s clones were created using the standard deviation of the orbital elements determined from the first integration using the JPL solution in 141 CE (${\sigma }_{e, a (\text{AU}), q (\text{AU}), i {(}^{°}), \omega {(}^{°}), \Omega {(}^{°})}~[6\times {10}^{-6},2\times {10}^{-3},7\times {10}^{-5},6\times {10}^{-3},6\times {10}^{-2},6\times {10}^{-2}]$). The junction of the two samples of clones integrated over the periods 15 000 BCE- 141 CE and 141 CE - 1994 CE is presented in Fig. 2 (bottom panel).

The clones’ evolution is similar for the two initial orbital solutions back to 8000 BCE, at which point the standard deviation of the clones generated in 141 CE from the Yeomans & Kiang (1981) solution starts to increase much faster than that of the clones created in 1994 CE based on the JPL solution. The abrupt dispersion of the clones around 8000 BCE, 9600 BCE, and 10 200 BCE in the Yeomans & Kiang (1981) 141 CE starting condition integration are all attributable to close encounters of ~40% of the test-particles with Jupiter. Our simulations are consistent with the hypothesis of a close encounter between 1P/Halley and Jupiter about 8000 to 12 000 yr ago that significantly modified the orbital evolution of the comet. The dynamical lifetime of 1P/Halley is probably much lower than its physical lifetime, which is estimated to be fewer than 1000 revolutions (Whipple 1987).

Since the date of this close approach depends sensitively on the initial conditions used for the integration, we are not able to set a clear upper time limit for our meteoroid stream simulations. We suggest that if the exact location of 1P/Halley is uncertain prior to 1404 BCE, the comet orbit also becomes uncertain after 10 000 yr of backward integration. This affects our choice of timeframe in our meteoroid simulations; we will return to this topic in Sect. 10.

Fig. 2 Time-evolution of the nominal (green) and average (red) semi-major axis of 1000 clones of comet 1P/Halley generated in1994 CE (top) based on the JPL J863/77 solution and 141 CE (bottom) based on the solution of Yeomans & Kiang (1981). The blue curves represent the three standard deviation range about the average semi-major axis. See text for details.

2.4.3 Resonances

Mean motion resonances (MMRs) with Jupiter play a major role in modulating the level of activity of the η-Aquariid and Orionid showers. In particular, the 1:6 MMR with Jupiter is likely responsible for ancient and recent Halleyid outbursts (Sato & Watanabe 2007; Rendtel 2007, 2008; Sekhar & Asher 2014; Kinsman & Asher 2020). Sekhar & Asher (2014)’s analysis of the resonant behavior of the Orionid meteoroid stream highlighted the influence of the 1:6 and 2:13 MMR with Jupiter. The authors examined the evolution of the comet’s 1:6 and 2:13 resonant arguments over 6 000 yr and concluded that 1P/Halley was trapped in the 1:6 and 2:13 MMR with Jupiter from 1404 BCE to 690 BCE and 240 BCE to 1700 CE, respectively.

Figure 3 presents an extension of Sekhar & Asher (2014)’s computation of the resonant arguments of 1P/Halley between 12 000 BCE and 1986 CE. Here, the comet’s motion was integrated for 14 000 yr using the orbital solutions of Yeomans & Kiang (1981). For this ephemeris solution, though it is uncertain prior to 1404 BCE (see Sect. 2.4.1), we can still see that the comet has spent on the order of a few thousand years within the 1:6 MMR with Jupiter since 12 000 BCE. From our analysis in Fig. 3, we find that the times when 1P/Halley has been in the 1:6 and 2:13 resonances since 1404 BCE are the same as found by Sekhar & Asher (2014).

In addition to being influenced by MMRs, 1P/Halley is known to have undergone Kozai resonances during its long-termevolution (Kozai 1962, 1979; Emel’Yanenko 2001). The comet’s perihelion argument librates between 47° and 133°, while the longitude of the ascending node varies between −18° and 82°. Kozai librations of the comet’s orbit might be helpful in constraining the formation of the Halleyids meteoroid streams (McIntosh & Hajduk 1983).

Fig. 3 1:6 (top) and 2:13 (bottom) resonant arguments for 1P/Halley based on integrations of the nominal orbit using the orbital solutions of Yeomans & Kiang (1981). The comet is in the resonance during times when the resonant argument librates, as opposed to when it circulates. The time period when the orbit is less constrained, prior to 1404 BCE, is shown in gray.

3 Halleyids

Comet 1P/Halley produces two annual meteor showers on Earth: the η-Aquariids in early May and the Orionids in late October. The η-Aquariids occur at thedescending node of the comet, while the Orionids occur at the ascending node. Modern observations of the η-Aquariids and Orionids date back to the late 19th century (Egal et al. 2020). However, analysis of ancient observations of meteor showers suggests that both showers have been active for at least a millennium, and likely longer (Ahn 2003, 2004; Imoto & Hasegawa 1958; Zhuang 1977).

The quantitative characteristics of the Halleyids meteor showers since 1985 are described in detail in Egal et al. (2020). That work documented the showers’ average activity, duration, profile shape, annual variations, and individual apparitions based on visual, video, and radar data collected since 1985, 2011, and 2002, respectively. The purpose of Egal et al. (2020) was to provide a consistent and multi-instrumental monitoring of the two meteor showers over several decades in order to constrain models of 1P/Halley’s meteoroid streams. That study is the main source we use to validate our numerical model of the Halleyids in the present work.

The next subsections briefly review the general characteristics of the Halleyid showers as determined in Egal et al. (2020). The results presented here, along with additional figures presented in that paper (annual variations and intensity profiles), are used throughout the rest of this study to calibrate our model with meteor observations. As usual, the activity of the showers are characterized by the zenithal hourly rate (ZHR), which represents the number of meteors a single observer would observe per hour under standard reference conditions.

3.1 η-Aquariids

The η-Aquariids are usually active between 35° and 60° solar longitude (SL), with the highest meteor rates recorded between 44° and 50°. The average maximum ZHR lies around 65 to 70 meteors per hour, but in the interval since 1985 the shower displayed two notable outbursts (ZHR > 100), in 2004 and 2013. The average profile of the η-Aquariids and Orionids was found to be variable from year to year. Over the period 2002 to 2019, the η-Aquariids average activity profile as a function of time is asymmetric, with a rise in intensity steeper than the subsequent decline.

The mass index s of the shower is also variable each year as a function of time during the shower’s period of activity. Lower mass indices are generally measured close to the main peak of activity, while higher values are found at the beginning and end of the shower. Estimates of s derived from visual observations range between 1.72 and 1.94 (Rendtel 1997; Dubietis 2003), and from 1.75 to 1.95 in radar records (Blaauw et al. 2011; Campbell-Brown & Brown 2015). The average mass index of the shower lies around 1.87 to 1.9 (Dubietis 2003; Egal et al. 2020).

3.2 Orionids

The Orionids are generally observed between 195° and 220° SL, and present a broad maximum between 206° and 211°. The average profile of the Orionids is more symmetric than the η-Aquariids around the maximum of activity.

Maximum meteor rates range between 20 and 40 meteors per hour, with occasional outbursts reaching two to four times these values (e.g., in 2006 and 2007). Recent CMOR measurements reinforce the hypothesis of a ~12 yr cycle in the Orionids activity. However, additional and longer term observations are required to definitively confirm this periodicity.

As with the η-Aquariids, the Orionids mass index varies with the solar longitude and by year. Estimates of s were found to vary between 1.46 or less (in 2006) and 1.96 in most years, with an average value of 1.87 (Dubietis 2003; Rendtel 2008). Estimates of the mass index from radar observations range between 1.65 and 1.93 from CMOR measurements (Blaauw et al. 2011), and around 1.95 in MAARSY data (Schult et al. 2018).

The comparison of visual, video, and radar observations since 2002 highlights the existence of multiple peaks in the Halleyid activity profiles, which vary by location and intensity for each apparition of the showers. The consistency between measurements of different systems supports the idea that there is no significant size sorting of the particles within the stream. However, the changing activity profiles indicate the existence of a complex fine structure within the meteoroid stream, as found in previous work (Egal et al. 2020).

4 Previous models

Many characteristics of the η-Aquariid and the Orionid meteor showers have challenged researchers for decades. In particular, the similar activity level and duration of the two showers (with such different nodal distances of the comet at the present day) cannot be explained with a toroidal meteoroid stream model. If the meteoroid streams were simply centered on the present orbit of the comet, the shower intensity and duration should be much higher for the η-Aquariids than for the Orionids. The toroidal model also cannot explain the gradual displacement of maximum activity observed for some Orionid apparitions (see review in Egal et al. 2020).

4.1 The shell or ribbon-like model

Perhaps the most successful model of the Halleyids is the theoretical “shell” model developed by McIntosh & Hajduk (1983). This model assumes that the evolution of meteoroids results in belts or “ribbon-like” structures intersecting Earth’s orbit at the time of the η-Aquariid and Orionid meteor showers.

The formation and the thickening of each meteoroid belt is explained by gravitational perturbations induced by Jupiter. In Yeomans & Kiang (1981)’s ephemeris, the variations of 1P/Halley’s longitude of the ascending node ΔΩ induced by Jovian perturbations range between 0.2° and 2.2°. If we assume that meteoroids are spread all along the comet’s orbit, then about 1/6 (~ P_Halley∕P_Jupiter) of the meteoroid stream is subject to the full range of ΔΩ perturbations induced by Jupiter at each revolution of the planet. After a few revolutions, the meteoroid stream will be distributed along a belt characterized by different Ω values. In other words, dust particles tend to diffuse to where the comet was a long time ago, along the present comet orbit, as well as on orbits that will be likely reached by the comet in the future. The thickness of each individual belt is expected to be relatively uniform along the intersection of the meteoroid’s nodes with Earth’s orbit (McIntosh & Hajduk 1983).

The shell model offers a natural explanation for the similar duration and intensity of the contemporary Halleyids showers. The libration of the argument of perihelion of the orbit suggested by Kozai (1979) creates a superposition of the meteoroid belts over time, also explaining the presence of multiple peaks of activity in the η-Aquariid and Orionid annual profiles. To reproduce the five Orionid activity peaks reported by Hajduk (1970), this model would require five libration cycles of the comet representing a period of about 100 000 yr (McIntosh & Hajduk 1983).

The theoretical shell model was investigated via numerical simulations by McIntosh & Jones (1988). Several sets of about 500 particles were ejected from the comet in 1404 BCE using Yeomans & Kiang (1981)’s orbital solution and integrated until 1986. The authors observed a quick elongation of the particles’ nodal locations, as well as the formation of a fine structure within the stream. This work resulted in global confirmation of the ribbon-like structure suggested by McIntosh & Hajduk (1983), but the similar duration and intensity of the meteor showers were not explained from the numerical simulations. In addition, this work suggested that particles ejected after 1404 BCE and with masses larger than 1 g would not intersect Earth’s orbit today. Therefore, the model suggests that such large meteoroids, recently ejected from 1P/Halley, would not contribute to the present day Orionid meteor shower (McIntosh & Jones 1988).

The results of McIntosh & Jones (1988) show that the fine structure of the stream can be created much faster than the 100 Kyr suggested by McIntosh & Hajduk (1983), possibly in as little as a few thousand years. The authors suggested that the meteoroid stream observed today at Earth developed after a close approach with Jupiter about 20 000 yr ago, but that no conclusive evidence of this close encounter exists because of the uncertain comet ephemeris prior to 1404 BCE. In a further analysis of the Orionids using radar, TV, and visual observations, Jones et al. (1989) estimated a total mass of the Halleyids meteoroid stream of 1.3 × 10¹³ kg, resulting in an age of 2500 to 62 000 yr. The author’s preferred age estimate lies around 23 000 yr.

4.2 Previous dust trail simulations

4.2.1 General characteristics

Using a similar approach to McIntosh & Jones (1988), Ryabova (2003) integrated the motion of discrete simulated meteoroid trails ejected in 1404 BCE, 141 CE, 837 CE, and 1910 CE. Each trail was composed of 5000 particles, ejected with a tenth of the velocity predicted by Whipple (1951)’s method and were within a 70° sunward cone. Particles of masses 1 and 0.001 g were simulated, with a density of 0.25 g cm⁻³. The orbital solution of Yeomans & Kiang (1981) was used for the comet ephemeris.

The simulations predicted a meteoroid spatial density structure closer to layers than filaments. From the trails’ nodal locations, Ryabova (2003) suggested that no simulated meteoroid of mass >0.001 g would contribute to the present day Orionid meteor shower. Particles producing the contemporary η-Aquariids were found to have been ejected prior to 837 CE at least, and only trails older than 1404 BCE could explain the observed duration of the shower. In these simulations, the nodal footprint of particles of 1g in mass are shifted compared to particles of 0.001 g. The mass distribution of the two showers was therefore expected to be dissimilar.

4.2.2 Resonant trails

More recently, classical dust trail simulations of 1P/Halley’s meteoroid streams were undertaken with the goal of identifying the years when past and future outbursts of the η-Aquariid and the Orionid meteor showers may occur. After the observation of the unexpectedly strong Orionid activity in 2006, the role of MMRs with both Jupiter and Saturn in producing Halleyid outbursts was investigated by several authors.

McNaught & Asher (1999), Sato (2003), and Sato & Watanabe (2007), following the methodology described in Kondrateva & Reznikov (1985), simulated dust trails released from the comet at each perihelion passage since 1404 BCE. They identified the 1:6 MMR to be responsible for the 2006 Orionid outburst. From similar calculations, Sato & Watanabe (2014) predicted an enhanced η-Aquariids activity on May 6, 2013, due to the encounter with meteoroids ejected in 910 BCE and 1197 BCE. The outburst was recorded by radio, visual, and video techniques as predicted (e.g., Molau et al. 2013; Cooper 2013; Steyaert 2014; Campbell-Brown & Brown 2015), confirming the validity of the numerical model.

Adopting the same approach, Kinsman & Asher (2017) find a correlation between probable past outbursts of the η-Aquariids in CE 400-600 and Mayan records of meteor activity. The model applied to the Orionids revealed the existence of probable significant outbursts in 417 CE and 585 CE, due to 1:6 and 1:7 MMRs with Jupiter (Kinsman & Asher 2020).

The role of MMRs with Jupiter on the occurrence of Orionid outbursts was investigated more precisely by Sekhar & Asher (2014). The authors find the 2:13 and especially the 1:6 resonances to be responsible for enhanced meteor activity on Earth over several consecutive years. From the analysis of these resonant regions and specific dust trail calculations similar to Sato & Watanabe (2007), Sekhar & Asher (2014) explained the existence and date of outbursts observed in 1436–1440, 1933–1938, and 2006–2010 as being due to the encounter by Earth of meteoroids trapped in the 1:6 MMR. The simulations showed that enhanced activity in 1916 and 1993 was caused by 2:13 MMR meteoroids, which they also predicted would cause a future Orionid outburst in 2070.

Sekhar & Asher (2014)’s simulations highlighted the importance of resonant structures in the meteoroid stream to explain the majority of the Orionids outbursts. Unfortunately, the analysis was not extended to the η-Aquariids due to the paucity of observations of the shower available at that time. With the publication of multiple decades of η-Aquariid observations (e.g., in Campbell-Brown & Brown 2015; Egal et al. 2020), the time is now ripe for the development, calibration, and validation of a more comprehensive model of 1P/Halley’s meteoroid stream.

4.3 The need for a new model

Several broad characteristics of the η-Aquariid and Orionid meteor showers have been successfully described by different studies detailed previously. McIntosh & Hajduk (1983)’s shell model, for example, offers a theoretical explanation for the similar intensity and duration of the two showers, as well as a general explanation for the existence of fine structure within the meteoroid streams. Past numerical simulations of 1P/Halley’s meteoroid stream generally agree with the stream evolution described by the shell model. In addition, the more recent recognition of the importance of resonant structures within the stream provides a natural explanation for outbursts and serves as an efficient tool to predict the occurrence of future Orionid and η-Aquariid outbursts.

However, several aspects of the streams remain unexplained. In particular, models developed so far do not provide quantitative estimates of the showers’ intensity variations, nor are they able to reproduce their comparable duration and strength or details of the radiant structure of the streams. While the rare major outbursts of the Halleyids are understood as being associated with MMRs, the mechanism behind the observed 11–12 yr modulation of annual shower activity previously reported for both streams is unknown.

Moreover, the range of ages of meteoroids comprising the stream for a given year’s apparition is only crudely constrained, as are variations in either the age or mass of particles as a function of time during any single return. Finally, the mass distribution as a function of particle size has not been predicted or explained by published models to date.

The broad aim of this project is to develop a new comprehensive numerical model of 1P/Halley’s meteoroid streams. Custom-made simulations of meteoroids ejected from the comet are calibrated using the 35 yr of Halleyid observations compiled in Egal et al. (2020), with the same model parameters applied to both the η-Aquariid and Orionid meteor showers. In particular, this model aims to reproduce the intensity profiles and annual ZHR variations of both showers in an automated way in order to forecast Halleyid activity until 2050. We also wish to estimate each shower’s age, radiant structure, and the presence (if any) of mass sorting within the streams. Finally, we seek to produce a model mass distribution for both streams over a large range of sizes appropriate for existing observations.

5 Simulations

Meteoroid stream simulations performed in this work follow the methodology described in Egal et al. (2019). Using an ephemeris of the comet’s past orbital behavior as its basis, millions of test particles are ejected from the comet nucleus and numerically integrated forward in time. The characteristics of Earth-impacting particles, as well as the dynamical evolution of the meteoroid stream, are recorded and analyzed. In order to calibrate the model with meteor observations, our primary simulation outputs are: the date and duration of potential meteor activity on Earth; the characteristics (age, size, etc.) of the trails involved; the stream’sstructure close to Earth; the radiants of Earth-intercepting meteors; and the approximate shower intensity.

The purpose of this section is to summarize the main model parameters chosen for our Halleyid meteoroid simulations. Additional information about the model is available in Egal et al. (2019).

5.1 Nucleus

In our model, we took 1P/Halley to have a spherical nucleus of 11 km diameter (Lamy et al. 2004), corresponding to an effective cross-section of about 95 km². We assumed a nucleus density of 0.3 g cm⁻³, which matches that of Rickman (1986). The geometric albedo of the nucleus was taken as 0.04 (cf. Sect. 2.1).

The motion of the comet was integrated between 1404 BCE and 2100 CE with an external time step of one day as described in Sect. 5.3. As in Egal et al. (2019), the ephemeris of each apparition of the comet was computed from the closest available orbital solution. Apparitions prior to 1986 CE were built from the set of osculating elements provided by Yeomans & Kiang (1981), converted into ecliptic J2000 coordinates. The comet coordinates since 1986 CE were determined from the JPL J863/77 solution². The model parameters are summarized in Table 1.

5.2 Stream simulation

As described in Egal et al. (2019), particles are ejected from the comet at each time step where the heliocentric distance is below 6 AU (Whipple 1987). Simulated meteoroids are released each day from the sunlit hemisphere of the nucleus, following the ejection velocity model of Crifo & Rodionov (1997). Crifo & Rodionov (1997)’s model is restricted to the sublimation of frozen water, responsible for ~90% of the volatiles flowing out of a comet below 3 AU (Combi et al. 2004). Our simulations extend meteoroid ejection beyond this distance in order to mimic the activity of 1P/Halley due to the sublimation of other volatiles farther from the Sun (Mazzotta et al. 2007; Womack et al. 2017). Particles are taken to have a constant density of 0.3 g cm⁻³ (Krasnopolsky et al. 1987).

Around 95 000 dust particles were ejected for each apparition of the comet between 1404 BCE and 2100 CE, for a total of about 4.48 million particles simulated. The particles were equally divided among the following three size/mass/magnitude bins: [10⁻⁴, 10⁻³] m, [10⁻⁹, 10⁻⁶] kg, [+10 to +3] mag; [10⁻³, 10⁻²] m, [10⁻⁶, 10⁻³] kg, [+3 to −4] mag; and [10⁻², 10⁻¹] m, [10⁻³, 1] kg, [− 4 to −10] mag. These bins approximately correspond to meteors detectable by radar (bin 1), visual and video means (bin 2), and as fireballs (bin 3).

Despite the unreliability of the comet’s location prior to 1404 BCE (Yeomans & Kiang 1981), additional simulations were conducted to assess the contributions of older trails, as some parts of the stream are almost certainly populated by ejecta older than 1404 BCE. 1P/Halley’s motion was integrated back to 7000 BCE using the oldest orbital solution of Yeomans & Kiang (1981). For these older ejections, we released a few thousand particles at each apparition of the comet at heliocentric distances inside 6 AU, covering the period [6982 BCE, 1404 BCE]. Between 18 000 and 6000 meteoroids were ejected at each return of the comet, for a total of about 1 million simulated particles. The number of particles ejected decreased arbitrarily as a function of the ejection epoch (cf. Table 2) to reduce the integration time of the oldest streams.

These older ejected particles were integrated forward until 2050 CE to investigate the influence of older trails on the Orionids and η-Aquariids. Because of the more uncertain ephemeris of the comet in this time period compared to the era post-1404 BCE, results of these simulations should be viewed with caution. The number of particles released at each comet apparition is summarized in Table 2.

Despite the large number of particles generated in all our simulations, the number of model meteoroids is negligible compared to the real number of meteoroids released by comet 1P/Halley at each perihelion return. To bridge the gap between our simulations and the real meteor shower activity, each simulated particle was assigned a weight representingthe actual number of meteoroids released by the comet under similar ejection circumstances. This weight depends upon, among other things, the particle’s radius and the comet’s dust production at these sizes at the heliocentric distance of ejection. In this work, we use the weighting scheme described in Egal et al. (2019). Details of the calibration of the weighting factors incorporating observations will be presented in Sect. 8.1.1.

5.3 Integration and data analysis

The comet and meteoroids are integrated with a 15th order RADAU algorithm (Everhart 1985), with a precision control parameter LL of 12 and an external time step of one day. The integration was performed considering the gravitational influence of the Sun, the Moon, and the eight planets of the Solar System. General relativistic corrections and nongravitational forces are also taken into account. The meteoroids are integrated as test particles in the stream, under the influence of solar radiation pressure and Poynting-Robertson drag. The Yarkovsky-Radzievskii effect is neglected since the size of our simulated particles does not exceed 10 cm (Vokrouhlický & Farinella 2000).

The integrations are performed in a Sun-centered, ecliptic J2000 coordinate system. However, the physical motion of the comet (and the particles) is only well described with heliocentric elements when the body is close to perihelion (typically within Jupiter’s orbit). Indeed, the periodic motion of the Sun around the Solar System’s barycenter induces oscillations of the comet’s heliocentric elements when computed far from perihelion. These oscillations are an artifact of the non-inertial heliocentric frame and are only apparent in the osculating elements. They can be removed by estimating the body’s orbital elements in the (inertial) barycentric frame. For this reason, the results presented in this work (particles orbits, nodes location) are determined and presented in the Solar System’s barycentric frame for convenience.

In order to avoid confusion, we employ some specific terminology in our simulation analysis. The terms “impactors” and “impacting particles” are applied to simulated particles approaching Earth within the distance ΔX and time ΔT defined in Sect. 7. Physically, these represent test particles whose actual positions in space best represent meteoroids that may actually collide with the Earth. The term “nodal crossing location” will refer to the particle’s position when physically crossing the ecliptic plane. “Nodes,” “node location,” or “nodal footprint” refer to the ecliptic position of the particle’s ascending and descending nodes as determined by their osculating orbit at the epoch of interest, even when the particle itself is far from the ecliptic plane.

6 Stream evolution

Before moving to an in-depth investigation of the simulation results, we present in this section a reconnaissance analysis of the meteoroid stream’s evolution. For this broad overview of general stream characteristics, we ejected 1000 particles at each apparition of the comet between 1404 BCE and 2050 CE as described in Sect. 5. Each particle’s position, orbit, and node location is recorded at a regular time step of 10 yr during the whole integration period, and every year between 1980 CE and 2030 CE. Three independent sets of simulations were conducted for particles in our three bins, representing sizes of 0.1–1, 1–10, and 10–100 mm.

6.1 Current structure

Figure 4 presents the current location of the simulated meteoroid stream, including all meteoroids that have nodal crossings between 1990 CE and 2020 CE, to enhance its visibility using this reconnaissance survey. The theoretical shell model presented inMcIntosh & Hajduk (1983) and McIntosh & Jones (1988) is reproduced at the top of Fig. 4 for comparison.

The simulated meteoroid stream intersects the Earth’s orbit at two locations and reproduces the ribbon-like shape predicted by McIntosh & Hajduk (1983). When looking at the edges of the meteoroid belt, the fine structure of the stream is clearly noticeable. As a first approximation, the gross behavior of the meteoroid stream is in good agreement with the McIntosh & Hajduk (1983) and McIntosh & Jones (1988) shell model.

Figure 5 shows the nodal footprint of all the simulated particles in 2020 CE, at the location of the η-Aquariid (top panel) and the Orionid showers (bottom panel). The stream nodes intersect Earth’s orbit between 40° SL and 54.5° SL at the time of the η-Aquariids, and between194° SL and 211° SL at the time of the Orionids. Trails ejected between 1404 BCE and 600 CE reproduce much of the current η-Aquariid activity. The present Orionids, observed between 195° SL and 220° SL are, in contrast, not reproduced as well. The Orionids seem to be composed of a larger percentage of older particles, ejected before 600 BCE (see Fig. 5). Thus we conclude that simulations of meteoroid trails ejected since 1404 BCE can reproduce the showers’ peak date but probably underestimate their observed duration.

In Fig. 5, we see that particles ejected between 1404 BCE and 600 BCE contribute to the present Orionid meteor shower. Simulations of McIntosh & Jones (1988) and Ryabova (2003) implied that meteoroids ejected after 1404 BCE and of masses higher than 1 and 0.001 g, respectively, have not approached Earth’s orbit in recent times. To compare our model with these results, we present in Fig. 6 the node locations of particles with masses higher than 1 g (filled circles) and between 0.001 and 1 g (empty circles). As expected, the nodal footprint is reduced when removing the smallest particles, and the width of oursimulated stream does not reproduce the shower’s observed duration. However, and in contrast with Ryabova (2003), we find thatparticles belonging to these mass ranges may still approach Earth’s orbit at the current epoch. Differences between our model and these earlier efforts may be due to different assumed ejection conditions (location, direction, and velocity).

From Figs. 5 and 6, we see that, in the case of the η-Aquariids, the simulations predict the onset of the shower at low SLs with a mixture of trails ejected at different epochs. However, by the end of the shower activity period, it consists primarily of the oldest trails. The simulated Orionids display a weaker gradient in ejection age as the Earth crosses the meteoroid stream, but the existence of such a gradient is hard to assess because of the overlapping locations of the trail’s nodes.

Fig. 4 Top panel: schematic representation of the Halleyid shell model as presented in McIntosh & Hajduk (1983) and McIntosh & Jones (1988). Bottom panel: 3D location of the simulated particles between 1990 CE and 2020 CE, color coded as a function of the particles’ ejection year.

Fig. 5 Node locations in 2020 CE for all particles as a function of their ejection date. The daily positions of the Earth on its orbit are indicated by blue circles. Numbers along the orbit correspond to the SL with a time step of 5 days. This distribution represents a crude estimate of the overall duration of the shower and peak activity location expected at the return in 2020 CE.

Fig. 6 Node locations in 2020 CE for particles of masses 10−6−10−3 kg (empty circles) and 10−3−1 kg (filled dots) as a function of the year of ejection.

6.2 Time evolution

To illustrate the evolution through time of the simulated streams, we independently analyzed the location of the nodes of each ejected trail as a function of time. As an example, we show the evolution of the meteoroid stream ejected around 1404 BCE until 2050 CE in Fig. 7. In this figure, we compile the locations of the particles’ descending (top panel) and ascending (bottom panel) nodes at different epochs from their ejection time (in green) to 2050 CE (in orange).

As noted in McIntosh & Jones (1988), the nodes of the ejected trail (in green) quickly spread, exceeding 2 AU in extent within 3000 yr. The median location of the nodes from the 1404 BCE ejection at each epoch is indicated with black dots in Fig. 7. The median node location of the 1404 BCE trail has changed considerably since the particles were ejected. The median descending node (top plot) is located well outside the Earth’s orbit in 1404 BCE and migrates inside the planet’s orbit after a few thousand years, intersecting the Earth’s orbit and contributing to η-Aquariid activity from roughly 0–1000 CE.

In contrast, the median ascending node of the 1404 BCE trail in Fig. 7 (bottom) starts much closer to the Earth’s orbit at the time of ejection but drifts outward with time to a present day distance of about 1 AU from the Earth. Here we see that the trail would significantly contribute to Orionid activity circa 1404–1000 BCE but not as much today. The evolution of the median location of every meteoroid trail ejected since 1404 BCE is presented in Appendix B. Each of the other simulated trails behave very much like the 1404 BCE trail, with little difference as a function of the particle sizes.

Figure 8 presents a zoom around Earth’s orbit of the median node locations for trails composed of 0.1 to 1 mm particles. The median node trajectories are in agreement with the 2020 CE nodal footprints observed in Fig. 5, confirming the picture that some particles remain on orbits that the comet was on a long time ago. When observing the intersection of the median ascending nodes with Earth’s orbit, the gradient in age of the Orionids is not clear. Instead, we observe a separation between trails ejected after or prior to 1000 BCE. The existence of such an age gradient needs to be re-examined, including the simulation of older meteoroid trail ejections from 1P/Halley.

Fig. 7 Evolution as a function of time (color coded as a function of the epoch) of the descending nodal footprint representing the η Aquariids (top) and ascending nodal footprint representing the Orionids (bottom) of a simulated trail ejected in 1404 BCE. The particle sizes are between 0.1 and 1 mm. The Earth’s orbit in 2020 CE, computed in the barycentric frame, is represented with a blue line.

6.3 Time to disperse around the orbit

Up to this point, we have investigated the evolution of the meteoroid streams by looking at the particles’ nodes at different epochs, which, because of their long orbital periods, are often somewhat different from the locations at which they actually cross the ecliptic plane (their “nodal-crossing positions” as defined in Sect. 5.3). In the previoussection, we found that trails ejected prior to 600 CE and 600 BCE could contribute to the present η-Aquariid and Orionid activity, respectively. However, this lower age-limit of the stream is valid only if these trails have had time to stretch along the whole orbit before the present epoch, which we verify below.

For several ejection epochs and particle sizes, we estimated the time T_S required for a newly ejected trail to entirely spread along the orbit. After ejection, we identified the average orbital period (μ) and standard deviation (σ) of the simulated particles and considered two particles of periods, P₁ = μ + 5σ and P₂ = μ − 5σ. After n revolutions, particle 1 lags particle 2 byΔT = n(P₁ − P₂). We therefore approximated T_S to the time required to satisfy ΔT = P₁ (i.e., T_S = P₁P₂∕(P₁ − P₂)). In this formulation, the influence of radiation and gravitational perturbations, which could accelerate the particles’ spread, is not taken into account.

For several trails ejected between 1404 BCE and 2050 CE, we estimated an average T_S time of 310 yr for 0.1 to 1 mm particles, 640 yr for 1 to 10 mm particles, and 1400 yr for particles of 10 to 100 mm in size. The dependence of T_S on meteoroid size is not surprising. Indeed, the combined influence of the radial radiation pressure and the meteoroids velocity relative to the nucleus (higher for small particles) tend to disperse the meteoroids immediately after ejection; the large meteoroids remain closeto the nucleus for a longer period of time, while small particles extend ahead and behind the comet more quickly, a general result of the modification of the orbits of meteoroids ejected from parent comets (Pecina & Simek 1997).

In our simulations, every trail ejected prior to 600 CE (for each size bin) has had time to spread along the whole orbit. The implication of these stretching time calculations is that it is reasonable to investigate the general evolution of the stream through the position of the particle nodes instead of the nodal crossing locations for trails ejected before 600 CE. However, this approach might not be valid for analyzing the fine structure of the stream, which is conditioned by planetary perturbations.

Fig. 8 Evolution of the median location of the descending (top panel) and ascending (bottom panel) nodes for each trail ejected between 1404 BCE and 1986 CE, color coded as a function of the epoch. Empty circles refer to trails ejected prior to 1000 BCE, filled circles to trails ejected between 1000 BCE and 1 CE, and crosses and the plus symbol to trails ejected after 1 CE. Symbols mark the median location with a time step of 70 yr ending in 2020 CE. Most trails have median locations outside the plot by 2020 CE. Particle sizes are 0.1–1 mm.

7 Extracting particles crossing their nodes and/or impacting Earth

The goal of our simulations is to convert our individual ejected test particles into equivalent shower activity at Earth. To accomplish this conversion, we must decide which particles represent meteoroids that can intersect the Earth. Since many of our test particles are from very old ejections, we do not expect the exact location of such particles to be highly accurate except in a statistical sense. In this case, the osculating elements of the particle may be a less accurate but more robust measure of the intersection of trails of a certain age with the Earth’s orbit. Using the nodal location based on the osculating orbit as described earlier can therefore provide a good global view of the age and extent of the overall trail, in particular since we have good number statistics.

In contrast, to achieve the highest fidelity in terms of actual particle numbers encountering the Earth in a given year, we must choose a volume box around the Earth in which to declare a particle an actual impactor. This population will have far fewer numbers but reproduce better year-to-year variations in activity as this approach captures features such as meteoroids residing in MMRs. We next describe the details of how we declare a simulated particle an impactor.

7.1 Impacting particles

During the integration, particles that approach Earth within a distance ΔX = V_rΔT are considered impactors. Here, V_r is the relative velocity between the planet and the particle and ΔT is a time parameter that depends on the accuracy of the shower prediction. For this study, ΔT was fixed at 12 h, leading to a distance criterion of 1.88 × 10⁻² AU for impactor selection.

7.1.1 Statistics

Overall, about 4% of the particles were flagged as impacting Earth and removed from the integration. Most of the impactors struck Earth when the MOID between the comet and the planet was at its lowest range (<0.02 AU), that is to say around 800 BCE for the Orionids and 500 CE for the η-Aquariids. Only 222 of the ~4.5 million simulated particles were impactors over the period 1980 to 2019. Among them, 62 particles approached Earth at the time of the Orionid meteor shower and 160 at the time of the η-Aquariids. Most years, fewer than five particles met the criteria of an impactor. With such low numbers, no comparison of the simulations with observed meteor rates can be attempted directly; other methods to increase our statistics must be considered.

When arbitrarily increasing the selection criteria ΔX to 0.05 AU (an increase of a factor of two), as commonly adopted for meteor shower forecasting (cf. review in Egal 2020), 1268 particles are retained as impactors over the period 1980 to 2050. The unweighted number of particles selected as a function of time for both showers is presented in Fig. 9. Despite the increase in ΔX, most of the annual apparitions of η-Aquariids and Orionids in this period contain fewer than ten impacting particles. Due to the small number of particles retained and the lack of realistic weights applied to the simulations (cf. Sect. 8.1.1), again no reliable intensity prediction can be generated from Fig. 9. So we again conclude that our statistics must be increased by other means. But first we will validate our simulations by comparing these initial impactors with observations.

7.1.2 Validation

Figure 9 shows impactors between 1980 and 2020. To the zeroth order, the larger the number of model impactors in a given year, the higher we expect η-Aquariid and Orionid activity to be. Almost all the particles impacting Earth in a given year were found to be ejected during a single apparition of the comet. The ejection year of the impactors since 1980 is shown in Fig. 10. For each year (abscissa), the ejection dates of particles impacting Earth are marked with circles for the η-Aquariids and squares for the Orionids. The color of the symbol represents the number of simulated impactors for that year (low activity in black, strong activity in yellow – see sidebar). The red and blue rectangles illustrate time periods when comet 1P/Halley was trapped in the 1:6 and 2:13 MMR with Jupiter.

Our model reproduces the basic features of the streams. In particular, no model particles ejected after 600 CE and 800 BCE become η-Aquariids or Orionids (respectively) on Earth at the current time, a finding we expect because of the evolution of 1P/Halley’s node away from the Earth after these times. In addition, the years with higher numbers of model impactors (i.e., 2013–2014 for the η Aquariids and 2006–2007 for the Orionids) correspond to ejecta potentially trapped in the 1:6 resonance, as predicted by Sato & Watanabe (2007, 2014), Rendtel (2008), and Sekhar & Asher (2014).

These first comparisons are encouraging and suggest our basic model captures the important details of the stream. However, no quantitative comparison with meteor observations can be performed without weighting the simulated meteoroids and increasing the number of particles selected. Section 8 addresses this issue.

7.2 Annual nodal crossings distribution

One approach to increase particle statistics in our analysis is to increase the temporal criterion ΔT. The temporal ΔT selection amounts to replacing each simulated particle p₀ by a swarm of particles on the same orbit, centered around p₀ and crossing the ecliptic plane at p₀’s nodal location during a period of 2ΔT. Small ΔT values (of a few hours or a few days) are required to reproduce short-lived activity variations of a simulated meteor shower. On the other hand, high ΔT parameters (of several months or even years) are useful for the investigation of the general evolution of a meteor shower’s activity, incorporating good statistics but at the expense of temporal resolution.

In this section, we consider a large ΔT criterion of half of a year, meaning that we retain all the test particles crossing the ecliptic plane at less than 0.05 AU from Earth’s orbit over an entire year (i.e., ±6 months from when the Earth intersects the shower node). With this relaxed selection threshold, our number statistics improve dramatically.

We find that a few hundred to a few thousand particles belonging to the η-Aquariids are retained every year (at least 200). For the Orionids, at least a few tens of (about 40) and sometimes a few hundred particles per year encounter the Earth over the period 1980 to 2050. In some specific years (e.g., 2002, 2003, 2004, 2014, and 2015), no particles ejected after 1404 BCE approach Earth at the time of the Orionids.

Figure 11 shows the stacked nodal crossing locations of particles traversing the ecliptic plane in 2020. An animation of the nodal crossing locations of our simulated particles ejected between 1980 and 2050 is available as supplementary material. The particles are color coded as a function of the ejection epoch. Meteoroids ejected when the comet was trapped in the 1:6 resonance are highlighted with red circles.

These plots show an overview of the years of potentially enhanced activity together with the age of the trails involved. From the raw number of nodes in the simulation located around Earth’s orbit, enhanced Orionid rates are expected in 2006 and 2007, along with relatively lower rates in 1997–1998 and 2010. The annual number of particles reaching the descending node suggests increased activity for the η-Aquariids in 1994, 2003–2004, 2012–2013, 2023–2024, and 2045–2046.

These probable outburst years, based on the nodal crossing locations, are in good agreement with years showing higher numbers of impactors (cf. Sect. 7). From the animation, we also see that the density of the nodal footprint around Earth is sufficient to produce meaningful activity profiles of the η-Aquariid meteor shower every year. With the selected ΔX and ΔT values of 0.05 AU and±6 months, several Orionid activity profiles around the outburst years have sufficient particle numbers for good profiles, but further expansion ofnumber statistics using different methods are required to produce reliable predictions of Orionid profiles and peak intensity in other years.

Fig. 9 Number of particles identified as impactors between 1980 and 2050. Striped and dashed boxes refer to η-Aquariids and Orionids ejected after 1404 BCE while the solid-colored boxes represent particles ejected before 1404 BCE. Vertical arrows show the years of recent η-Aquariid (red) outbursts in 2004 and 2013 and Orionid (blue) outbursts in 2006-2007.

Fig. 10 Particle ejection year from 1P/Halley as a function of impacting year at the Earth for the η-Aquariids (circles) and the Orionids (squares). Symbols are color coded as a function of the number of impactors (see sidebar). When several trails from different ejection epochs are involved in producing impactors (e.g., in 2012 and 2013), connected symbols frame the ejection period. The rectangles illustrate the periods when the ejection epoch occurred with 1P/Halley trapped in the 1:6 (red) and 2:13 (blue) MMR with Jupiter. Meteoroids ejected during these intervals may remain in the MMRs and form comparatively denser “clumps” in the stream. Vertical lines show the years of recent η-Aquariid (red) outbursts in 2004 and 2013 and Orionid (blue) outbursts in 2006-2007.

Fig. 11 Example of the distribution of nodal locations of particles when crossing the ecliptic plane in 2020 CE as a function of their time of ejection from 1P/Halley. The regions where meteoroids have their ascending node (Orionids) and their descending node (η-Aquariids) are shown in magnified form in the middle and right panels of the figure. The position of Earth along its orbit is represented by blue dots and dates (gray arrows indicating month and dots the day of each month) or SLs (middle and right panels). Particles are color coded as a function of the ejection epoch. Meteoroids ejected when the comet was trapped in the 1:6 resonance are highlighted with red circles. Trails ejected between 3000 BCE and 1404 BCE are presented in gray as these are associated with the time frame when 1P/Halley has a less certain ephemeris. The animated figure is provided online.

8 Postdictions

The final step in our modeling is to calibrate our model and determine particle weighting for the Halleyid meteoroid streams through direct comparison with meteor observations. We first computed synthetic intensity profiles (model shower activity as afunction of time around the maximum of each shower, each year) from our simulations. This provides a quick analysis of the meteor showers’ shape, duration, average activity, and year-to-year variations to compare to the measurements provided in Egal et al. (2020). In addition, the simulated radiants and particles’ mass and size distributions at Earth can be compared to meteor observations, all of which both constrain and validate the simulations.

8.1 Simulated activity profiles

8.1.1 Computation

With the ΔX and ΔT values selected in the previous sections, tens to thousands of particles crossing the ecliptic plane every year are included for activity profile computations. Each particle is assigned a SL, corresponding to the date when the Earth is the closest to the particle’s node. Within a given SL bin, the total number of particles is converted into an equivalent flux $\mathcal{F}$, which is then transformed into a ZHR by Koschack & Rendtel (1990)’s relation: $$\text{ZHR}=\frac{{\mathcal{A}}_{\text{s}}\mathcal{F}}{(13.1r-16.5){(r-1.3)}^{0.748}},$$(2)

where ${\mathcal{A}}_{\text{s}}$ is the typical surface area for meteor detection by a visual observer in the atmosphere at the ablation altitude (${\mathcal{A}}_{\text{s}}~37200{\text{km}}^2$), and r is the measured differential population index. We assumed a common population index of 2.5 for the two showers. This lies in between the values of 2.46 and 2.59 generally measured for the η-Aquariids and the Orionids (Egal 2020).

For the ZHR model estimates to be meaningful we first require a realistic, physical simulated flux. This is accomplished by weighting each particle as discussed, in Sect. 5.2, through calibration against measurements.

The initial weighting scheme applied to our simulation follows the methodology detailed in Egal et al. (2019). In this approach, each particle weight depends on: the initial number of particles of a given radius ejected at a given epoch; 1P/Halley’s dust production at the heliocentric distance of ejection, measured via Afρ; and the differential size frequency distribution of the particle radii after ejection, characterized by the differential size index u.

In practice, the final weight is normalized by the number of particles that share the same ejection characteristics (epoch of ejection and size) and is increased for low heliocentric distances as the comet dust production peaks around the perihelion. The weighting varies significantly as a function of the differential size index u at ejection that is used. Additional details regarding the weighting scheme are available in Egal et al. (2019).

Given our spatial and temporal accuracy, the aforementioned weighting scheme proves to be effective when the raw number of simulated particles in a shower profile exceeds a few hundred (i.e., in the case of the η-Aquariid analysis). For cases with fewer simulated particles, additional empirical weights are found to improve the quality of our postdictions. These additional weights are somewhat ad hoc and were found semi-empirically but are physically motivated. They are necessary for obtaining sufficient statistics. The three complementary weights applied to our simulations are as follows.

First, there is a coefficient C₁ proportional to the inverse distance of the particle’s node to the Earth’s orbit. The somewhat large ΔX criterion of 0.05 AU, which increases the number of simulated particles selected, is corrected by this coefficient, which gives more importance to nodes close to Earth’s orbit. This weight is an approximation of the impact parameter described in Moser & Cooke (2008).

Second, there is a coefficient C₂ proportional to 1∕n, with n being the number of revolutions completed by the trail since its ejection. This parameter, similar to the fM coefficient of Asher (1999), gives more weight to young trails that are known to be more dense than older streams. The inclusion of this corrective factor has been considered in past meteoroid stream simulations, for example by Watanabe & Sato (2008) to model the Draconid meteor shower.

Third, there is a coefficient function of the time δT_p between the time a particle crosses the ecliptic plane and the date of the prediction (by definition, δT_p < ΔT). As noted in the previous section, a ΔT of 0.5 yr leads to no simulated particles available for certain Orionid apparitions. so this parameter is increased to 1.5 yr in some cases (amounting to concatenation of three consecutive years of simulations at the epoch being investigated). The importance of these additional particles is, however, decreased by a weight of 1∕γ^δT_p, with γ a positive integer to be determined in our calibration fits.

The inclusion of these additional weights has only a small effect on years with high predicted activity, but it considerably improved our results for years with only a few simulated particles approaching Earth’s orbit. Of these three complementary weights, the last correction is arguably the least justified. However, we find that the inclusion of these additional particles and their relative weighting when determined empirically does improve model agreement with measurement of several apparitions of the Orionids. We also note that stacking several years of nodal crossing is a common practice in meteor shower simulations; this method has been used by several authors to study the general structure of streams when limited by small-number statistics (e.g., Jenniskens & Vaubaillon 2008; Šegon et al. 2017).

In summary, the weighting scheme applied to our raw model simulations has three tunable parameters, namely the particle size distribution index at ejection, u, the contribution of particles crossing the ecliptic plane at different dates, γ, and a final normalization coefficient, K₂, used to scale all the ZHR predictions to the observed activity level of the meteor showers as in Egal et al. (2019). Since our simulations aim to model the Halleyid meteoroid streams self-consistently, we applied the exact same weighting scheme to the simulated particles comprising both the η-Aquariid and the Orionid showers. In this way we can test the robustness of the model fits by limiting the weighting coefficients to be the same for thetwo showers, as is expected to be the case on purely physical grounds.

8.1.2 Calibration

The best weighting solution is determined by calibrating the activity profiles on past observations of the meteor showers. At present, the longest and most consistent continuous monitoring of the Halleyids is provided by CMOR, operational since 2002 (Egal et al. 2020). Our simulated activity profiles are calibrated to CMOR observations at 29 MHz, first for the η-Aquariid meteor shower (which has more particles) and then refined with the Orionid simulations.

The determination of the weights was performed in a semi-automated way. We first estimated the value of the particles’ size distribution index u using a particle swarm optimization (PSO) algorithm (Kennedy & Eberhart 1995). The PSO was run with two different cost functions, related respectively to the profiles’ main peak ZHR and shape. The first cost function compared the computed and observed maximum ZHR reached for all the apparitions between 2002 and 2019. The second cost function reflected the difference between the observed and the simulated activity profiles when the main peaks are scaled to the observed maximum at each apparition of the shower.

With the PSO, we find the u value that best reproduces the shower activity profiles’ shape to be 3.4 (corresponding to a mass index of 1.80), in good agreement with Giotto in situ measurements close to the nucleus (see Sect. 2.3). On the other hand, we find that the interannual variations of the main peak intensity are better modeled with a size index of 3.9 (i.e., a mass index of 1.97). The best u is also found to vary both as a function of the number of apparitions considered in the calibration and the specific years investigated. Since one of the main goals of this work is to predict future η-Aquariid and Orionid outbursts, we adopted a size index u of 3.9 (which provided a better match to the peak intensities) for the rest of this study.

The application of the three additional weights (C₁, C₂, and γ) was performed manually in a second, three-pronged step. (1) The weight modulation as a function of the particle’s nodal distance or ejection epoch was selected so as to improve specific Orionid apparitions. (2) Each particle was divided by the approximate number of revolutions it had performed since its ejection n. This operation is found to improve the average level of the η-Aquariid maximum activity between 1985 and 2020. (3) The γ coefficient mainly affects the year-to-year variations of the main peak ZHR, which are reduced as γ decreases. We find the best value to lie around γ = 3, which does not affect the η-Aquariid activity variations but improves the statistics of lower-activity Orionid apparitions.

8.2 Results

In the following subsections, we investigate the agreement of our simulated activity profiles with video, visual, and radar observations of the η-Aquariid and Orionid meteor showers presented in Egal et al. (2020). In this analysis, we consider every simulated trail ejected since 3000 BCE as described in Sect. 5.2. The inclusion of trails ejected before the 1404 BCE apparition, which helps improve our results in some cases, is discussed in Sect. 10.

8.2.1 Average activity profile

The average activity profile of the simulated η-Aquariid and Orionidapparitions between 2002 and 2019 is presented in Fig. 12. The average profiles as measured with the CMOR 29 MHz and 38 MHz systems, with the IMO VMN, and with visual observationsare plotted for comparison. The simulated activity profiles were scaled to the observations of the η-Aquariid shower, and the same scaling factor was applied to the modeled Orionid profile.

The average η-Aquariid profiles generated from the simulations are in good agreement with the observations. The modeled activity peaks around 45.5° SL, with the rise in activity more sudden than the post-maximum decrease in intensity in accordance with measurements. The simulation’s post-maximum activity matches that of the observed profiles well, while the pre-maximum shows some divergence at SLs below 42°. The simulated profile underestimates the total shower duration (of around 25 days) by about 5 days, a potential indicator that the shower contains material even older than the oldest included in our simulations (3000 BCE).

With the same scaling factors, the average Orionid activity deduced from our simulations peaks around 25 meteors per hour, that is to say about a factor of 2.5 times less than the η-Aquariids. Though our simulations do not exactly reproduce the observed activity ratio between the two showers (of about 2), we obtain, for the first time, comparable meteor rates with a self-consistent stream model. The simulated profile peaks around the 208° to 209° SL at the observed location of the first peak of activity but underestimates the plateau and falling portion of the activity curve. The total duration of the shower is also underestimated by about 12 days. These mismatches between the model and shower are an even stronger indication that a more significant fraction of the meteoroids comprising the Orionids are from ejections earlier than 3000 BCE compared to the η-Aquariids.

In summary, the good match between our simulated profiles and the observations suggests that most of the η-Aquariid activity can be reproduced with trails ejected since 3000 BCE. On the other hand, as expected from the discussion in previous sections, the Orionid average profile cannot be modeled completely with the same trails and requires the addition of older material. Trails ejected after 3000 BCE contribute to the core of the Orionid activity profile, but the width of the profile and the second maximum of activity reported around 210° SL are probably caused by material ejected before 3000 BCE.

Fig. 12 Average activity profile of the η-Aquariid (top) and Orionid (bottom) meteor showers, as measured with the CMOR 29 MHz system (blue) and the 38 MHz system (green) between 2002 and 2019, the VMN (orange) between 2011 and 2019, and visual observations from the IMO between 1988 and 2007 (Rendtel & Arlt 2008) or between 2002 and 2019 (VMDB in black). Video and radar observationswere scaled to visual measurements as described in Egal et al. (2020). The best-weighted simulated average profile usingone degree SL binning is drawn in red. We note the factor of two difference in the ordinate scale between the top and bottom plots.

8.2.2 Average mass index

As we generate absolute fluxes from our model scheme, we are also able, for the first time for any shower model to our knowledge, to generate predictions of the synthetic mass distribution of a shower at the Earth. The cumulative mass distribution of the simulated showers, computed based on weightings of all the particles encountering Earth between 2002 and 2019, is presented in Fig. 13. Only particles contributing to the main peak of activity of the η-Aquariids (44°–46° SL) and Orionids (208°–2010° SL in the simulations) were included. The simulated mass distributions display breaks and changes in the slope, though an overall rough power-law trend is evident. This is particularly true for the Orionids because of the smaller number of simulated particles in that sample.

The differential mass index s of each shower is estimated over the mass ranges [10⁻⁴, 10⁻⁶] kg, [10⁻⁶, 10⁻⁷] kg, and [10⁻⁷, 10⁻⁸] kg, corresponding to visual (M ~ [−2, +3]), video (M ~ [+3, +6]), and radar (M ~ [+6, +8]) magnitudes, respectively. Average s values estimated from the model span from 1.78 to 1.98 for the η-Aquariids and from 1.44 to 2.09 for the Orionids, in relatively good agreement with observations of these showers (cf. Table 3).

Estimates in the visual range, of, respectively, 1.84 and 1.92, are particularly consistent with the average value of 1.87 reported by Dubietis (2003) and Rendtel (1997). A key prediction of the modeled mass distribution is that while the overall distribution is roughly a power law, the specific power law value that is obtained depends very much on the mass interval used for measurements. In particular, in the radar and video mass ranges, our modeling suggests significant “kinks” in the distribution. This naturally can lead to quite different mass index values depending on the exact lower and upper bound to the mass range being measured and may account for some of the variance in mass index for the showers reported in the literature.

8.2.3 Annual variation

The ZHR of the main peak of each simulated η-Aquariid and Orionid apparition between 1985 and 2020 is presented in Fig. 14. The showers’ individual activity profiles between 2002 and 2019 are provided in Appendix C. As noted in Sect. 8.2.1, our model slightly underestimates the observed average activity level of the Orionids, potentially a consequence of it being an older shower compared to the η-Aquariids. In Fig. 14 and Appendix C, the overall modeled and observed profiles have been rescaled so they reach the same maximum level. Because of this correction, the relative intensity of the two simulated showers, which from Sect. 8.2.1 is about 2.5, is not respected anymore. However, we find this scaling necessary to compare the predicted and observed activity profile shapes and maximum rates with clarity; it is likely due to the model not including older ejections which proportionately make up more of the Orionids.

The yearly variations of the simulated η-Aquariid maximum ZHRs in Fig. 14 are in generally good agreement with the observations. The selected weighting scheme reproduces the observed 2004 and 2013 outbursts, as well as the overall ZHR variations observed since 2000. On the other hand, this model underestimates the shower activity prior to 1998, between 2006 and 2010, and in 2017. In our simulations, these apparitions are more dependent on the presence of trails ejected before 1404 BCE, which, we recall, are in an era where the ephemeris for 1P/Halley is less certain.

We can improve our agreement with the maximum η-Aquariid ZHRs observed prior to 1995 by increasing the weights of older trails in our simulations. However, a consequence of this improvement is that our models’ 2003 and 2014 ZHRs are raised by a factor of about 1.5. Since the purpose of this work is to predict future Halleyid outbursts and their peak activity, we prefer to underestimate the usual intensity level of the shower than to lose reliability around the years of expected enhanced activity. This situation demonstrates that our weighting scheme is subjective and not unique, and that it is conditioned by which observations are used and most emphasized in the calibration.

For the Orionids, despite the incompleteness of the activity profiles above 210° SL, our simulations reproduce the variations of the Orionid main peak intensities prior to 1996, around the 2006 and 2007 outbursts, and since 2010. The ZHR values at other times may be underestimated, with the exception of the 1998 apparition, which is significantly overestimated. In the IMO Visual Meteor DataBase (VMDB), we find the observed maximum Orionid rate in 1997 reached 27 meteors per hour³, suggesting that our model also overestimates the shower activity for that particular year. As was the case for the η-Aquariids, modifications of the weighting scheme (and especially the corrective factor depending on the particles nodal distance) would allow us to remove the excess intensity in 1997 and 1998. However, such modifications also worsen the fit for other apparitions of the two showers and were not retained here.

For most apparitions since 1985, our chosen model reproduces the reported maximum ZHR to within 30% or less.This translates to an absolute error of about 20–30 meteors per hour for each shower, depending on the meteor observation considered (this is indicated by the red region of Fig. 14). We therefore consider this model successful in reproducing the global yearly activity variations of the η-Aquariids and Orionids.

Fig. 13 Modeled mass distribution of the η-Aquariids (top) and Orionids (bottom) based on simulated meteoroids encountering the Earth between 2002–2019 near the time of the main peak of activity. The logarithm of the cumulative number of simulated meteors (weighted) is drawn as a function of the logarithm of the particles mass (in kg). Estimates of the differential mass index s in the visual(red), video (green), and radar (blue) mass ranges are provided for comparison. See text for details.

8.2.4 Individual apparitions

The strengths and limitations of our model are well described by the comparison of the simulated average profiles and year-to-year ZHR variations with meteor observations discussed previously. However, individual apparitions of the η-Aquariids and Orionids are also of interest, and some of them will be briefly detailed in this section.

As expected from the average activity profile, our simulated apparitions of the η-Aquariids (see Appendix C) are generally most consistent with observations. Except for 2019, the simulated maximum of activity matches the reported main peak, though with some dependence on the observations considered. In several years the model even reproduces the shape of the observed profile in some detail (e.g., in 2012, 2015, 2016, 2017, etc.).

The model is less successful in reproducing the Orionid profiles, especially for SLs above 210°. We interpret this as a likely consequence of it being more composed from ejections predating our oldest modeled ejections. For this reason, the ~12 yr periodicity of the Orionid main peak location (Egal et al. 2020) cannot be reliably investigated with our simulations. The pre-maximum branch of activity is in general good agreement with the observations, with some notable exceptions (e.g., in 2002–2003, 2008, and 2014), suggesting it is composed of relatively younger material. However, we note that the observed profiles from the various detection networks often differ as well, so varying observational sensitivities may also account for some of the differences.

The characteristics of the modeled outbursts (η-Aquariids in 2004 and 2013, Orionids in 2006 and 2007) are detailed in Appendix D. The activity profiles of the showers are presented along with the trails involved and the cumulative mass distributions. In these four cases, the simulated profile reproduces the intensity, duration, and shape of the observed activity well. The mass distributions are variable and s indices also vary with the range of masses considered (cf. figures in Appendix D).

In agreement with Rendtel (2007), Sato & Watanabe (2007), and Sekhar & Asher (2014), we find several trails ejected while 1P/Halley was in the 1:6 resonance to be responsible for the 2006 and 2007 Orionid outbursts. In our model, trails ejected in 1196 BCE and 1264 BCE are particularly strong contributors to these years. We also find multiple trails ejected after 3000 BCE to have moderate contributions to the 2004 and 2013 η-Aquariid outbursts. The main contributor to the 2004 outburst is the 314 BCE trail and, to a lesser extent, the 910 BCE and 1964 BCE streams. In our model, the 1128 BCE trail is the principal source of the 2013 outburst, together with particles ejected in 2731 BCE, 1196 BCE, 835 BCE, and 762 BCE (among others). This result slightly diverges from that of Sato & Watanabe (2014), who identified the 1196 BCE and 910 BCE trails as the main cause of this outburst.

Fig. 14 Annual variations of the η-Aquariid (top) and Orionid (bottom) peak ZHRs for the years 1985–2019. The maximum rates recorded by CMOR 29 MHz (blue dots) and 38 MHz (green diamonds) are plotted along with the results of the VMN system (orange triangles) and visual observations (black dots), as in Egal et al. (2020). The envelope of estimated simulated peak variations, with an arbitrary uncertainty of ±20 meteors per hour that allows the framing of most Orionid yearly maximum ZHR rates, are presented in red. The phase of the Moon for each year at the time of the shower peak is also shown, as this significantly affects the quality and quantity of visual and video data.

Fig. 15 Radiants of the η-Aquariids (top) and Orionids (bottom) measured with CAMO (red circles) and the GMN (black triangles) video networks in 2019. The average shower radiant location in one degree SL bins determined from the wavelet analysis of CMOR data is identified with black crosses. For GMN and CAMO, the formal radiant uncertainty was found using the method of Vida et al. (2019b). Simulated radiants computed for particles ejected before and after 1450 BCE are plotted with cyan and dark blue circles, respectively.

8.3 Radiant

Another way to validate our model with real observations is through comparing the meteor radiants. Measurements of the 2019 η-Aquariids (top panel of Fig. 15) and Orionids (bottom panel) recorded by video cameras of the Canadian Automated Meteor Observatory(CAMO, see Vida 2020) and the Global Meteor Network (GMN, cf. Vida et al. 2018, 2019a) are presented in Fig. 15. The figure includes the general radiant location computed from CMOR using stacked data from 2002–2016, as well as the modeled 2019 radiants of particles ejected from 1P/Halley before (cyan) or after (dark blue) 1450 BCE.

For both showers, the simulated radiants of trails ejected after 1450 BCE are located close to the observed radiants. The area delineated by the radiant dispersion of all the modeled η-Aquariids covers the area of average CMOR radiants for longitudes below 295°, but poorly reproduces radiants located above this longitude. The simulated Orionid dispersion reproduces the radiant structure around CAMO measurements well, but the simulated meteors are shifted or absent for higher and lower latitudes. A complete modeling of the radiant structure probably requires the inclusion of meteoroid trails ejected before 3000 BCE, particularly for the parts of the shower away from the core. For example, the simulation of older trails might be necessary to explain the formation of the observed Orionid radiant tail reported by Jenniskens et al. (2016).

9 Forecast

The comparison of our new Halleyid model with existing meteor observations in the previous section provides us with confidence that the model can reproduce the overall profile of the showers and their interannual peak activity. With the inclusion of trails ejected since 3000 BCE, our model satisfactorily reproduces the average activity profile, mass distribution, radiant dispersion, and yearly intensity variations of the η-Aquariid meteor shower. The same model applied to the Orionids leads to a 30% error of the shower’s average activity level and underestimates its duration, but reproduces the first observed peak of activity and the year-to-year maximum ZHR variation. Our simulations indicate that the Orionids are probably composed of much older material, as highlighted in McIntosh & Jones (1988) or Ryabova (2003), particularly during the descending branch of the shower.

Our simulations do not reproduce the totality of the observed activity profiles or radiant structure, perhaps not surprising given that the comet’s ephemeris and dust production rates have been uncertain and/or completely unknown for much of its history. Despite these limitations, the general trends of the main peak ZHR annual variations are reasonably well reflected in our model. The models’ success in the postdiction of characteristics of the η-Aquariids and Orionids since 1985 gives us enough confidence to forecast the future activity of these showers. In the following sections, we present the predictions of our model for Halleyid activity in the period 2020 to 2050.

9.1 Annual variation

The main peak ZHR of each η-Aquariid and Orionid apparition between 2020 and 2050 based on our calibrated model is shown in Fig. 16. Consistent with Sekhar & Asher (2014), we expect no strong Orionid activity over the analyzed period. The predicted meteor rates vary between 10 and 40 meteors per hour, corresponding to the usual range of ZHRs reported for the shower (Egal et al. 2020). In contrast, Fig. 16 shows enhanced η-Aquariid activity in 2023–2024 and 2045–2046, reaching peak rates similar to those observed in 2004 and 2013. These will be examined below.

9.2 Outburst characteristics

Details of the predicted future η-Aquariid outbursts are presented in Fig. 17. The 2045 and 2046 predictions are similar (in profile shape, trails involved, and mass distribution) and therefore only the 2046 outburst is detailed in the figure. The four showers are predicted to display noticeable activity between 40° and 55° in SL, with a maximum around 44.1° ± 0.5° (2023), 45.2° ± 0.5° (2024), and 44.8° ± 0.5° (2045 & 2046).

Most of the expected activity in 2023 is produced by particles ejected from 1P/Halley around 390 BCE, with the presence of potentially resonant particles ejected in 829 BCE. Additional trails (e.g., in 2215 BCE or 218 CE) have only a moderate influence on the profile’s shape. The main peak of activity in 2024 arises from the 985 BCE trail, with smaller contributions of particles ejected during the 835 BCE, 1058 BCE, and 314 BCE apparitions. The combined contributions of several trails ejected while the 1P/Halley was trapped in the 1:6 MMR with Jupiter increases the expectation of enhanced η-Aquariid rates in 2024. The 2045 and 2046 enhancements are caused by more recent material, ejected around 218 CE. In contrast with the η-Aquariid outburst in 2013, these two years will mainly show enhancement due to material ejected outside the 1:6 MMR, making them more similar to the situation in 2004.

As discussed earlier, the modeled mass distribution of the predicted outbursts does not strictly follow a power law, showing slope variations for different mass ranges. Observations of the mass index s therefore need to be computed with caution, carefully taking into account the mass range involved. Predicted s values span from 1.85 to 1.96 in 2023, 1.59 to 2.09 in 2024, 1.72 to 2.1 in 2045, and 1.73 to 1.96 in 2046 over typical radar, video, and visual mass intervals.

Fig. 16 Predicted variations of the η-Aquariid (top) and Orionid (bottom) main peak ZHRs between 2020 and 2050.

10 Discussion

10.1 Limitations

Most of the limitations of the proposed model stem from the modifications implemented to increase the number of particles retained. The larger ΔX and ΔT thresholds selected can lead to increasing prediction errors, artificially increasing ZHRs in years around an outburst year or confusing the identification of the dominate contributing trail for a particular apparition. Future meteor observations will help to update the weighting constraints developed in this work, which in turn should lead to refinement in modeled future activity of the Halleyids.

The inclusion of trails ejected between 3000 BCE and 1404 BCE in our analysis also introduces some uncertainty. Section 2.4.2 suggests that the comet’s orbit is known with some accuracy over the past several millennia, even if its exact location is not precisely determined before 1404 BCE (Yeomans & Kiang 1981). To investigate the influence of 1P/Halley’s ephemeris prior to 1404 BCE in our simulations, we ejected 5000 particles of 1–10 mm in size from different clones of the comet around 3000 BCE. Five clones of the comet were selected from the analysis of Sect. 2.4.2 to cover the full range of possible semi-major axes (from 16.6 to 18.8 AU). The 25 000 simulated meteoroids were ejected as described in Sect. 5 and integrated until 2020 CE, when the particles’ nodal location and average activity profiles were compared.

This analysis revealed that the choice of the comet ephemeris around 3000 BCE has a moderate impact on the showers’ average activity profile (mainly on the modeled total duration), with a potential variation of 20 to 40% of the meteor flux in periods away from the main peak of activity. The selection of different comet clones prior to 1404 BCE leads to only a small variation of the modeled peak ZHR for a given year and does not modify the characteristics of the predicted future η-Aquariid outbursts. This reflects the fact that the majority of the core of both streams is composed of material ejected more recently than 1404 BCE, though this is more true of the η-Aquariids than the Orionids. In addition, we find that the differences in the modeled radiants using different clones prior to 1404 BCE still produce radiants which lie within the observed radiant dispersion.

We therefore consider the results presented in Sects. 8 and 9 as valid despite the inclusion of trails ejected from the nominal comet orbit between 3000 BCE and 1404 BCE, where the comet ephemeris ismore uncertain. Despite the model’s limitations, our simulations reproduce the main characteristics of the η-Aquariid apparitions since 1985, as well as the general trends in ZHR variations measured for the Orionids, including the years of outburst activity, the mass distribution of the showers and the location of the shower peaks. Our model shower durations are smaller than is observed, consistent with the notation that much of the material at the extreme start and end of the activity periods are older than 3000 BCE.

These successes increase our confidence that our model is able to accurately predict activity of the two showers through 2050, and allows for the examination of their potential behavior over a longer time period. In particular, we can apply the model to investigate the long-standing question as to the cyclic variation of the η-Aquariid and Orionidactivity, which we discuss below.

10.2 Long-term periodicity

The existence of a periodicity in either the Orionid peak activity or the timing of maximum (or both) has been proposed by several authors in the past (cf. Egal 2020). In particular, a cyclic variation of the shower’s maximum rates of about 12 yr was suggested by McIntosh & Hajduk (1983) and also by Dubietis (2003). The analysis of CMOR measurements between 2002 and 2019 by Egal (2020) found a possible periodicity in the SL of the Orionid peak timings with a period of 11.88 yr. However, no strong conclusion about the shower’s cyclic variations in either peak time or activity could be established without additional and longer-term meteor observations. Egal et al. (2020) find no clear periodicity in the annual activity level of the η-Aquariids from analysis of the CMOR data. With our new model reproducing calibrated fluxes as a function of time for the Halleyid streams, we are now in a position to examine this question in more detail.

The annual variation in the observed (black line) and modeled (dark blue line) η-Aquariid and Orionidpeak ZHRs between 1985 and 2050 is given in Fig. 18. A sinusoidal signal was fit to both the predicted and observed ZHR estimates using the PSO method. The amplitude of the sine wave was allowed to vary to improve the fit with the modeled ZHR. This best fit solution to the model for each shower is presented in red in Fig. 18.

The best PSO solution for the peak Orionid activity has a period of 11.83 yr, close to the value of 11.88 found in Egal (2020) and Jupiter’s orbital period of 11.86 yr. Though the agreement between the PSO fit and observations showing outburst activity is not perfect, the sine function fit generally reproduces the observed and modeled intensity variations of the shower.

To investigate the underlying physical cause of this periodic behavior, we explored the contribution to each annual return from particles ejected when 1P/Halley was trapped in the 1:6 MMR with Jupiter. The result is illustrated with the blue boxes in Fig. 18. There is a clear correlation between the activity level of the shower and the number of particles reaching the Earth ina given year that were ejected when 1P/Halley was in the 1:6 MMR. This suggests that the 1:6 MMR induces a ~12-yr cyclic variation of Orionid activity.

The annual η-Aquariid activity is better described with a sinusoid having a period of 10.72 yr. The sinusoidal signal fit is in good agreement with the modeled peak ZHR since 2005, but does not reproduce earlier measurements of the shower well. When examining the possible contributionof 1:6 resonant particles, we note that these particles are involved in some observed and predicted outbursts (like in 2013 and 2024) but not all of them (e.g., in 2004 or 2045–2046). The existence of a periodic variation of the η-Aquariid annual activity is suggested by our analysis, but it is not as convincing as the Orionid fits. This periodicity needs to be confirmed by future radar and optical observations of the shower.

Fig. 17 Details of model predictions for future η-Aquariid outbursts. From top to bottom: predicted model ZHR profile, relative contribution of simulated trails as a function of SL for a given year and the model cumulative mass distribution. The middle trail plots show each ejection relative contribution by circle size (small to large) and color following the scale bar, with black showing trails with the lowest contribution and yellow trails the highest contribution. The colored horizontal periods show the intervals when 1P/Halley was in the 1:6 MMR in red and in the 2:13 MMR in blue. The mass distribution plots follow the conventions from Fig. 13, with estimates of the differential mass indices s in the visual(red), video (green), and radar (blue) mass ranges provided for comparison.

Fig. 18 Annual variation of the observed (in black) and modeled (dark blue line) Orionid (top) and η-Aquariid (bottom) peak ZHRs between 1985 and 2050. This figure represents a combination of Figs. 14 and 16. A sinusoidal curve of variable amplitude was fit to the modeled and observed meteor rates. This fit result is presented in red in the figure. The best match with the model yields a period of 11.83 yr in the case of the Orionids and 10.72 yr for the η-Aquariids. The relativecontribution of particles intercepted by the Earth in a given year, which were ejected while 1P/Halley was trapped in the 1:6 MMR with Jupiter, are shown with blue boxes.

11 Conclusions

We developed a new model of the meteoroid streams released by comet 1P/Halley and used it to characterize the associated η-Aquariid and Orionid meteor showers. Through simulation of more than 5 million particles and the application of a custom-made weighting solution, we modeled the showers’ activity profiles, radiants, mass distributions, and maximum ZHR variations between 1980 and 2050. The simulation results were calibrated using the peak activity and profiles of 35 yr of Halleyid observations, self-consistently applying a unique weighting scheme to the two showers.

From the model we find that the η-Aquariid characteristics (activity profile, radiants, and annual variations) can be largely reproduced by material ejected between 3000 BCE and 600 CE. In particular, the core of the stream is almost entirely made up of ejecta from this interval. The wings of the activity are less well fit and suggest that even older meteoroids are present in the stream, particularly near the start and end of the activity period.

The model fit to the Orionid activity profile, in contrast, leads to an underestimate of the activity. This indicates that the Orionids are older than the η-Aquariids on average. The model is generally able to reproduce the first (ascending) part of the broad Orionid maximum but fails to reproduce the reported second activity peak around 210° SL and the descending branch of activity. Based on our model, the core and early portion of the Orionid activity are primarily populated by material ejected between 3000 BCE and 610 BCE. Despite this limitation, the simulations reproduce the general trends of the maximum interannual ZHR variations, providing us with some confidence in our activity forecasts.

From our calibrated model we produced estimates of the mass distribution of meteoroids in the core of both streams. While the overall distribution roughly follows a power law, our modeling predicts significant kinks in the distribution, notably in the mass range 10⁻⁵ to 10⁻⁷ kg. This may explain the variability in measurements of the mass index reported by different instruments for the Halleyid streams, as the values will depend sensitively on the upper and lower mass ranges used for a given fit.

Our simulations suggest that the age of the present day η-Aquariids slightly exceeds 5000 yr, while the Orionids are older. Future simulation efforts could better constrain the upper range to the Orionid age by simulating meteoroids ejected from several clones of 1P/Halley’s nominal orbit before 3000 BCE. This also implies that the minimal physical age for 1P/Halley is 5000 yr.

Our model provides predictions for the peak ZHRs expected from the Halleyids in the period 2020 to 2050. These predictions suggest no significant Orionid outbursts in this interval and four potential η-Aquariid outbursts in 2023, 2024, 2045, and 2046. The expected maximum ZHRs for the η - Aquariid outbursts vary from 120 to 160 meteors per hour, with a 30% confidence on the predicted rates. Our model suggests that most of the η -Aquariid activity expected in 2024 originates from particles ejected while comet 1P/Halley was trapped in the 1:6 MMR resonance, similar to several previously observed η-Aquariid and Orionidoutbursts.

Among the MMRs with Jupiter and Saturn that can shape the long-term evolution of the Halleyid meteoroid stream, the 1:6 MMR with Jupiter has the strongest effect on the stream’s structure. The 1:6 MMR is a significant, but not the only, source of Halleyid outbursts. This is consistent with the results of other authors (cf. Sato & Watanabe 2007, 2014; Rendtel 2008; Sekhar & Asher 2014) and also suggests that the resonance induces a cyclic variation of the Orionids’ average ZHR levels with a period of about 11.8 yr. A potential periodicity of 10.7 yr for the η-Aquariids needs to be confirmed by future observations of the shower. Visual, video, and radar observations of each η-Aquariid and Orionidapparition are strongly encouraged to provide additional constraints on numerical models.

Movie

Movie of Fig. 11 Access here

Acknowledgements

This work was supported in part by NASA Meteoroid Environment Office under cooperative agreement 80NSSC18M0046 and by the Natural Sciences and Engineering Research Council of Canada (Grants no. RGPIN-2016-04433 & RGPIN-2018-05659), and by the Canada Research Chairs Program. We are thankful to Aswin Sekhar for his careful review that helped improve the manuscript.

Appendix A 1P/Halley orbital motion

Fig. A.1 Left: ecliptic location of 1P/Halley’s ascending and descending node as a function of time, in the barycentric frame (ecliptic J2000). The Earth’s orbit is represented in blue. Right: time variations of the MOID (in AU) of comet 1P/Halley and Earth close to the comet’s ascending and descending node.

Fig. A.2 Time-evolution of the nominal (green) and average (red) orbital elements of 1000 clones of comet 1P/Halley generated in 1994 CE and based on the JPL J863/77 solution. The blue curves represent the one standard deviation range about the average semi-major axis.

Appendix B Average precession

Fig. B.1 Median location of the trails’ descending node (left panels) and ascending node (right panels) as a function of time for particles of sizes 0.1–1mm (top panels) and 10–100 mm (bottom panels). Each line represents a meteoroid trail ejected during a single apparition of the comet, between 1404 BCE and 1910 CE.

Appendix C Simulated activity profiles

Fig. C.1 Activity profiles of the η-Aquariids between 2002 and 2020, without the outbursts of 2004 and 2013. Outburst years are detailed in Appendix D.

Fig. C.2 Activity profiles of the Orionids between 2002 and 2019, without the outbursts of 2006 and 2007. Outburst years are detailed in Appendix D.

Appendix D Past Halleyid outburst characteristics

Fig. D.1 From left to right: simulated (red line) and measured ZHR profiles, relative contribution of simulated trails, and cumulative mass distribution. The trails’ relative contribution is indicated by the circles size (small to large) and color (black: low significance, to yellow: highest contribution). The red band in the middle panel indicates the times when 1P/Halley was in the 1:6 resonance and blue the 2:13, as discussed in Sect. 2.4.3.

References

All Tables

All Figures

Fig. 1 Raw (left panel) and phase-corrected (right) Afρ measurements of comet 1P/Halley during the 1986 apparition presented in Schleicher et al. (1998) and Fink (1994). The abscissa (X) represents the comet’s heliocentric distance to perihelion, arbitrarily negative for pre-perihelion measurements and positive for post-perihelion times ($X=(r_h-q)\frac{t-t(q)}{|t-t(q)|}$, with t(q) being the time at perihelion).

In the text

	Fig. 2 Time-evolution of the nominal (green) and average (red) semi-major axis of 1000 clones of comet 1P/Halley generated in1994 CE (top) based on the JPL J863/77 solution and 141 CE (bottom) based on the solution of Yeomans & Kiang (1981). The blue curves represent the three standard deviation range about the average semi-major axis. See text for details.
In the text

	Fig. 3 1:6 (top) and 2:13 (bottom) resonant arguments for 1P/Halley based on integrations of the nominal orbit using the orbital solutions of Yeomans & Kiang (1981). The comet is in the resonance during times when the resonant argument librates, as opposed to when it circulates. The time period when the orbit is less constrained, prior to 1404 BCE, is shown in gray.
In the text

	Fig. 4 Top panel: schematic representation of the Halleyid shell model as presented in McIntosh & Hajduk (1983) and McIntosh & Jones (1988). Bottom panel: 3D location of the simulated particles between 1990 CE and 2020 CE, color coded as a function of the particles’ ejection year.
In the text

	Fig. 5 Node locations in 2020 CE for all particles as a function of their ejection date. The daily positions of the Earth on its orbit are indicated by blue circles. Numbers along the orbit correspond to the SL with a time step of 5 days. This distribution represents a crude estimate of the overall duration of the shower and peak activity location expected at the return in 2020 CE.
In the text

	Fig. 6 Node locations in 2020 CE for particles of masses 10−6−10−3 kg (empty circles) and 10−3−1 kg (filled dots) as a function of the year of ejection.
In the text

	Fig. 7 Evolution as a function of time (color coded as a function of the epoch) of the descending nodal footprint representing the η Aquariids (top) and ascending nodal footprint representing the Orionids (bottom) of a simulated trail ejected in 1404 BCE. The particle sizes are between 0.1 and 1 mm. The Earth’s orbit in 2020 CE, computed in the barycentric frame, is represented with a blue line.
In the text

Fig. 8 Evolution of the median location of the descending (top panel) and ascending (bottom panel) nodes for each trail ejected between 1404 BCE and 1986 CE, color coded as a function of the epoch. Empty circles refer to trails ejected prior to 1000 BCE, filled circles to trails ejected between 1000 BCE and 1 CE, and crosses and the plus symbol to trails ejected after 1 CE. Symbols mark the median location with a time step of 70 yr ending in 2020 CE. Most trails have median locations outside the plot by 2020 CE. Particle sizes are 0.1–1 mm.

In the text

	Fig. 9 Number of particles identified as impactors between 1980 and 2050. Striped and dashed boxes refer to η-Aquariids and Orionids ejected after 1404 BCE while the solid-colored boxes represent particles ejected before 1404 BCE. Vertical arrows show the years of recent η-Aquariid (red) outbursts in 2004 and 2013 and Orionid (blue) outbursts in 2006-2007.
In the text

Fig. 10 Particle ejection year from 1P/Halley as a function of impacting year at the Earth for the η-Aquariids (circles) and the Orionids (squares). Symbols are color coded as a function of the number of impactors (see sidebar). When several trails from different ejection epochs are involved in producing impactors (e.g., in 2012 and 2013), connected symbols frame the ejection period. The rectangles illustrate the periods when the ejection epoch occurred with 1P/Halley trapped in the 1:6 (red) and 2:13 (blue) MMR with Jupiter. Meteoroids ejected during these intervals may remain in the MMRs and form comparatively denser “clumps” in the stream. Vertical lines show the years of recent η-Aquariid (red) outbursts in 2004 and 2013 and Orionid (blue) outbursts in 2006-2007.

In the text

Fig. 11 Example of the distribution of nodal locations of particles when crossing the ecliptic plane in 2020 CE as a function of their time of ejection from 1P/Halley. The regions where meteoroids have their ascending node (Orionids) and their descending node (η-Aquariids) are shown in magnified form in the middle and right panels of the figure. The position of Earth along its orbit is represented by blue dots and dates (gray arrows indicating month and dots the day of each month) or SLs (middle and right panels). Particles are color coded as a function of the ejection epoch. Meteoroids ejected when the comet was trapped in the 1:6 resonance are highlighted with red circles. Trails ejected between 3000 BCE and 1404 BCE are presented in gray as these are associated with the time frame when 1P/Halley has a less certain ephemeris. The animated figure is provided online.

In the text

Fig. 12 Average activity profile of the η-Aquariid (top) and Orionid (bottom) meteor showers, as measured with the CMOR 29 MHz system (blue) and the 38 MHz system (green) between 2002 and 2019, the VMN (orange) between 2011 and 2019, and visual observations from the IMO between 1988 and 2007 (Rendtel & Arlt 2008) or between 2002 and 2019 (VMDB in black). Video and radar observationswere scaled to visual measurements as described in Egal et al. (2020). The best-weighted simulated average profile usingone degree SL binning is drawn in red. We note the factor of two difference in the ordinate scale between the top and bottom plots.

In the text

Fig. 13 Modeled mass distribution of the η-Aquariids (top) and Orionids (bottom) based on simulated meteoroids encountering the Earth between 2002–2019 near the time of the main peak of activity. The logarithm of the cumulative number of simulated meteors (weighted) is drawn as a function of the logarithm of the particles mass (in kg). Estimates of the differential mass index s in the visual(red), video (green), and radar (blue) mass ranges are provided for comparison. See text for details.

In the text

Fig. 14 Annual variations of the η-Aquariid (top) and Orionid (bottom) peak ZHRs for the years 1985–2019. The maximum rates recorded by CMOR 29 MHz (blue dots) and 38 MHz (green diamonds) are plotted along with the results of the VMN system (orange triangles) and visual observations (black dots), as in Egal et al. (2020). The envelope of estimated simulated peak variations, with an arbitrary uncertainty of ±20 meteors per hour that allows the framing of most Orionid yearly maximum ZHR rates, are presented in red. The phase of the Moon for each year at the time of the shower peak is also shown, as this significantly affects the quality and quantity of visual and video data.

In the text

Fig. 15 Radiants of the η-Aquariids (top) and Orionids (bottom) measured with CAMO (red circles) and the GMN (black triangles) video networks in 2019. The average shower radiant location in one degree SL bins determined from the wavelet analysis of CMOR data is identified with black crosses. For GMN and CAMO, the formal radiant uncertainty was found using the method of Vida et al. (2019b). Simulated radiants computed for particles ejected before and after 1450 BCE are plotted with cyan and dark blue circles, respectively.

In the text

	Fig. 16 Predicted variations of the η-Aquariid (top) and Orionid (bottom) main peak ZHRs between 2020 and 2050.
In the text

Fig. 17 Details of model predictions for future η-Aquariid outbursts. From top to bottom: predicted model ZHR profile, relative contribution of simulated trails as a function of SL for a given year and the model cumulative mass distribution. The middle trail plots show each ejection relative contribution by circle size (small to large) and color following the scale bar, with black showing trails with the lowest contribution and yellow trails the highest contribution. The colored horizontal periods show the intervals when 1P/Halley was in the 1:6 MMR in red and in the 2:13 MMR in blue. The mass distribution plots follow the conventions from Fig. 13, with estimates of the differential mass indices s in the visual(red), video (green), and radar (blue) mass ranges provided for comparison.

In the text

Fig. 18 Annual variation of the observed (in black) and modeled (dark blue line) Orionid (top) and η-Aquariid (bottom) peak ZHRs between 1985 and 2050. This figure represents a combination of Figs. 14 and 16. A sinusoidal curve of variable amplitude was fit to the modeled and observed meteor rates. This fit result is presented in red in the figure. The best match with the model yields a period of 11.83 yr in the case of the Orionids and 10.72 yr for the η-Aquariids. The relativecontribution of particles intercepted by the Earth in a given year, which were ejected while 1P/Halley was trapped in the 1:6 MMR with Jupiter, are shown with blue boxes.

In the text

	Fig. A.1 Left: ecliptic location of 1P/Halley’s ascending and descending node as a function of time, in the barycentric frame (ecliptic J2000). The Earth’s orbit is represented in blue. Right: time variations of the MOID (in AU) of comet 1P/Halley and Earth close to the comet’s ascending and descending node.
In the text

	Fig. A.2 Time-evolution of the nominal (green) and average (red) orbital elements of 1000 clones of comet 1P/Halley generated in 1994 CE and based on the JPL J863/77 solution. The blue curves represent the one standard deviation range about the average semi-major axis.
In the text

	Fig. B.1 Median location of the trails’ descending node (left panels) and ascending node (right panels) as a function of time for particles of sizes 0.1–1mm (top panels) and 10–100 mm (bottom panels). Each line represents a meteoroid trail ejected during a single apparition of the comet, between 1404 BCE and 1910 CE.
In the text

	Fig. C.1 Activity profiles of the η-Aquariids between 2002 and 2020, without the outbursts of 2004 and 2013. Outburst years are detailed in Appendix D.
In the text

	Fig. C.2 Activity profiles of the Orionids between 2002 and 2019, without the outbursts of 2006 and 2007. Outburst years are detailed in Appendix D.
In the text

Fig. D.1 From left to right: simulated (red line) and measured ZHR profiles, relative contribution of simulated trails, and cumulative mass distribution. The trails’ relative contribution is indicated by the circles size (small to large) and color (black: low significance, to yellow: highest contribution). The red band in the middle panel indicates the times when 1P/Halley was in the 1:6 resonance and blue the 2:13, as discussed in Sect. 2.4.3.

In the text