Modeling the past and future activity of the Halleyids meteor showers

We present a new numerical model of the eta-Aquariid and Orionid meteor showers. The model investigates the origin, variability and age of the eta-Aquariid and Orionid apparitions from 1985 to the present day, in order to forecast their activity over the next several decades. Through the numerical integration of millions of simulated meteoroids and a custom-made particle weighting scheme, we model the characteristics of every eta-Aquariid and Orionid apparition between 1985 and 2050. The modeled showers are calibrated using 35 years of meteor observations including the showers activity profiles and interannual variability. Our model reproduces the general characteristics of the present-day eta-Aquariids, and part of the Orionid activity. Simulations suggest that the age of the eta-Aquariids somewhat exceeds 5000 years, while a greater fraction of the Orionids are composed of older material. The 1:6 mean-motion resonance with Jupiter plays a major role in generating some (but not all) Halleyid stream outbursts. We find consistent evidence for a periodicity of 11.8 years in both the observations and modeled maximum meteor rates for the Orionids. A weaker evidence of a 10.7 year period in the peak activity for the eta-Aquariids needs to be investigated with future meteor observations. The extension of our model to future years predicts no significant Orionid outburst through 2050 and four significant eta-Aquariid outbursts in 2023, 2024, 2045 and 2046.


Introduction
In recent years, thanks to modern numerical simulations of meteoroid streams, the timing of meteor outbursts have reached unprecedented accuracy (e.g., Asher 1999). In contrast, reliable estimates of meteor showers 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 lack of long-term, consistent records of the resulting meteor showers.
The modeling of comet 1P/Halley and Halleyids' meteor showers, i.e. the η-Aquariids and Orionids, 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 and Asher 2014).
In this work, we present the second part of our analysis of the Halleyids meteor showers. In Egal et al. (2020), we have documented the long-term (∼several decades) activity of the η-Aquariids and Orionids. The showers' activity profiles, shapes, and year to year variations were defined through a compilation of visual meteor observations since 1985, more recent measurements of the Canadian Meteor Orbit Radar (CMOR, cf. Brown et al. 2008Brown et al. , 2010 and the IMO Video Meteor Network (VMN) (Molau and Rendtel 2009). That paper reported results of multi-instrumental measurements of the Halleyids over seve-mail: aegal@uwo.ca eral decades using consistent analysis techniques. This longterm characterization of the Halleyids provides the needed observations to validate modeling of 1P/Halley's meteoroid complex.
The purpose of the present work is to apply a 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 years 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 for e.g., Egal et al. 2019), we constrain the model parameters using the η-Aquariids and Orionids intensity profiles and interannual variability since 1985. The simulation parameters are calibrated on these 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 Orionid meteors including origin, age and variability. After model calibration of the observed showers before 2020, we provide a forecast of the Halleyids activity until 2050.
The paper is divided to four main parts as follows: A) 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. B) 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 Halleyids models is provided in Section 4.
Article number, page 1 of 30 arXiv:2008.03589v1 [astro-ph.EP] 8 Aug 2020 A&A proofs: manuscript no. Egal2020c C) 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 (Section 6) and the predicted associated meteor activity on Earth (Section 7). D) 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 Orionid meteor showers is forecast through 2050. The successes and limitations of our model are discussed in Section 10.

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 for the 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 years. 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. The first observation of the comet has been suggested to date to more than 2 000 years 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 were 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 pertinent for meteoroid stream modeling.

Nucleus
The spacecraft flybys of 1P/Halley revealed that the comet's nucleus was asymmetric with an avocado-like shape having dimensions of 15x8x8 km (Whipple 1987) and a mean effective cross-sectional area of about 80-100 km 2 (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 1987b).
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 non-gravitational 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 −3 , with an upper limit of order 0.3 g cm −3 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 −3 . Work by Whipple (1987) suggested upper limits of 0.4 to 0.5 g cm −3 .

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 observations in 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 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 1987b).

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). 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. Afρ measurements are independent of the aperture radius ρ when the comet's radial brightness varies as ρ −1 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 ρ −1 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: 1. Schleicher et al. (1998)'s relation for θ ≤ 42°: log 10 A f ρ(θ) = −0.01807θ + 0.000177θ 2 2. Marcus (2007a,b) phase relation based on the Henyey-Greenstein function for θ > 42°T he left panel of Figure 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) Schleicher et al. (1998) andFink (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) t−t(q) |t−t(q)| , with t(q) the time at perihelion).
estimates of Fink (1994) were decoupled from phase angle effects via the composite dust phase function detailed above, and plotted along with Schleicher et al. (1998)'s measurements in Figure 1. Following the procedure of Egal et al. (2019), we fit 1P/Halley's Afρ profile of Figure 1 with a double-exponential function of the form: where r h is the comet's heliocentric distance at time t, q the perihelion distance, and γ 1,2 , X max and K 1 parameters which are to be determined by the fitting process. The best agreement with measurements of 1P's 1985 apparition gives K 1 ∼ 3500 cm, A f ρ(X max ) ∼ 105 285 cm, γ 1 = 2.97 (pre-perihelion branch), γ 2 = 0.99 (post-perihelion branch) and X max = −0.04 AU as illustrated by the black line in Figure 1. The comet dust production is asymmetric, with a maximum around perihelion and a steeper pre-perihelion profile compared to post-perihelion.

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 −16 g in mass: one silicate-rich and the other referred as CHON (Carbon, Oxygen, Hydrogen and Nitrogen) particles (Kissel et al. 1986a,b). The dust's bulk density was estimated to be 0.297 g cm −3 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 −8 − 10 −13 kg) but this value showed intersensor differences as well as variations with nuclear distance 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 (1987a). Estimates of s (the differential mass distribution exponent in dN=Cm −s dm, where dN is the number of grains between m and m+dm) ranged between 1.49 and 2.03 for masses between 10 −6 and 10 −13 kg, with an average value of 1.9 to 2.1 . 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 Section 3).

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

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 which 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 and Kiang 1981;Stephenson et al. 1985). It has been suggested that variations of the comet's non-gravitational acceleration might be responsible for the differences between its modeled and observed past motion (Sitarski and Ziolkowski 1987).
Different sets of osculating elements of the comet have been proposed in the literature (see for example Chang 1979;Brady 1982;Yeomans and Kiang 1981;Landgraf 1984Landgraf , 1986Sitarski and Ziolkowski 1987). To reproduce modern and ancient observations, Landgraf (1984) introduced a linear variation of the comet's nongravitational parameters. Yeomans and 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 in the comet's eccentricity during the Earth's close approach of 837 CE to match older observations. Yeomans and 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 and 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 and Ziolkowski 1987). Yeomans and Kiang (1981)'s solution for 1P/Halley's ephemeris between 1404 BCE and 1910 CE is usually taken as the definitive reference of the comet's orbital evolution for studies of its meteoroid streams (e.g., McIntosh and Jones 1988;Ryabova 2003;Sekhar and Asher 2014). The evolution of 1P/Halley's nodal location since 1404 BCE, computed with Yeomans and Kiang (1981)'s ephemeris (cf. Figure 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. Figure 2 in Egal et al. 2020), measuring the shortest Cartesian distance between the two osculating orbits when considering the pre-perihelion or the post-perihelion arcs of the comet's orbit.

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 years. Hajduk (1985) estimated a maximum dynamical age within the inner solar system of 200 000 years. The theoretical model of the Halleyids by McIntosh and Hajduk (1983) estimated the dynamical age of 1P/Halley at around 100 000 years. Sekhar and Asher (2014)'s numerical integrations from Yeomans and Kiang (1981)'s 240 BCE orbital solution showed that close encounters with Jupiter significantly modified the comet orbit about 12 000 years in the past.
We have simulated the past orbital evolution of the comet in a similar way. We use as starting conditions a recent orbital solution provided by the JPL, and assume constant nongravita-  Yeomans and Kiang (1981). The blue curves represent the three standard deviations range about the average semi-major axis. See text for details. tional acceleration parameters with time. The comet was integrated from the JPL 1 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 axis is presented in Figure 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 and Asher (2014)'s results, but does not take into account the manual corrections of 1P/Halley's ephemeris performed by Yeomans and Kiang (1981) around 837 CE. To check the validity of our analysis, we have also generated a second sample of clones from Yeomans and 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 and 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 (σ e, a (AU), q (AU), i (°), ω (°), Ω (°) ∼ [6 × 10 −6 , 2 × 10 −3 , 7 × 10 −5 , 6 × 10 −3 , 6 × 10 −2 , 6 × 10 −2 ]). The junction of the two samples of clones integrated over the periods 15000 BCE -141 CE and 141 CE -1994 CE, is presented in Figure 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 and 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 and Kiang (1981) 141 CE starting conditions 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 of 1P/Halley with Jupiter about 8000 to 12 000 years 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 years of backward integration. This affects our choice of timeframe chosen in our meteoroid simulations; we will return to this topic in Section 10.

Resonances
Mean motion resonances (MMR) 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 Orionid outbursts (Sato and Watanabe 2007;Rendtel 2007Rendtel , 2008Sekhar and Asher 2014;Kinsman and Asher 2020). Sekhar and 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 years, 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 and 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 years using the orbital solutions of Yeomans and Kiang (1981). For this ephemeris solution, though it is uncertain prior to 1404 BCE (see Section 2.4.1), we can still see that the comet spent of order a few thousand years within the 1:6 MMR with Jupiter since 12 000 BCE. From our analysis in Figure 3, we find that the times when 1P/Halley is in the 1:6 and 2:13 resonances since 1404 BCE are the same as found by Sekhar and Asher (2014).  Yeomans and 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 grey.
In addition to being influenced by mean motion resonances, 1P/Halley is known to have undergone Kozai resonances during its long term evolution (Kozai 1962(Kozai , 1979Emel'Yanenko 2001). The comet's perihelion argument librates between 47°and 133°, while the longitude of the ascending node varies between -18°a nd 82°. Kozai librations of the comet's orbit might be helpful in constraining the formation of the Halleyids meteoroid streams (McIntosh and Hajduk 1983).

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 the descending 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 . 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(Ahn , 2004Imoto and 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.

η-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 zenithal hourly rate 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 of 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 and Brown 2015). The average mass index of the shower lies around 1.87 to 1.9 (Dubietis 2003;Egal et al. 2020).

Orionids
The Orionids are generally observed between 195°and 220°S L, 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 to 40 meteors per hour, with occasional outbursts reaching 2 to 4 times these values (e.g., in 2006 and 2007). Recent CMOR measurements reinforce the hypothesis of a ∼12 years 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 confirm the existence of multiple peaks of activity during a given year's apparition in both showers, variable in location and intensity. 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 ).

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) can not 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 can not explain the gradual displacement of maximum activity observed for some Orionid apparitions (see review in Egal et al. 2020).

The shell or ribbon-like model
Perhaps the most successful model of the Halleyids is the theoretical 'shell' model developed by McIntosh and 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 and 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 and 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 Orionids activity peaks reported by Hajduk (1970), this model would require five libration cycles of the comet representing a period of about 100 000 years (McIntosh and Hajduk 1983). The theoretical shell model was investigated via numerical simulations by McIntosh and Jones (1988). Several sets of about 500 particles were ejected from the comet in 1404 BCE using Yeomans and 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 fine structure within the stream. This work resulted in global confirmation of the ribbon-like structure suggested by McIntosh and 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 and Jones 1988).
The results of McIntosh and Jones (1988) show that the finestructure of the stream can be created much faster than the 100 Kyr suggested by McIntosh and Hajduk (1983), possibly in as little as a few thousands of years. The authors suggested that the meteoroid stream observed today at Earth developed after a close approach with Jupiter about 20 000 years 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 13 kg, resulting in age of 2 500 to 62 000 years. The author's preferred age estimate lies around 23 000 years.

General characteristics
Using a similar approach to McIntosh and 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 g and 0.001 g were simulated, with a density of 0.25 g/cm −3 . The orbital solution of Yeomans and 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.001g 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.001g. The mass distribution of the two showers was therefore expected to be dissimilar.

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 mean motion resonances with both Jupiter and Saturn in producing Halleyids outbursts was investigated by several authors.
Following the methodology described in Kondrateva and Reznikov (1985), McNaught and Asher (1999) and Sato (2003), Sato and Watanabe (2007) 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 and 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 and Brown 2015), confirming the validity of the numerical model.
Adopting the same approach, Kinsman and Asher (2017) found a correlation between probable past outbursts of the η-Aquariids in CE 400-600 and Maya 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 mean motion resonances with Jupiter (Kinsman and Asher 2020).
The role of mean motion resonances with Jupiter on the occurrence of Orionid outbursts was investigated more precisely by Sekhar and Asher (2014). The authors found 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 and Watanabe (2007), Sekhar and 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 to cause a future Orionid outburst in 2070.
Sekhar and 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 η-Aquariids observation (e.g., in Campbell- Brown and Brown 2015;Egal et al. 2020), the time is now ripe for development, calibration and validation of a more comprehensive model of 1P/Halley's meteoroid stream.

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 and Hajduk (1983)'s shell model, for example, offers a theoretical explanation of 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 provide a natural explanation for outbursts and serve 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 year 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 is variations in either the age or mass of particles as a function of time during any one 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 years of Halleyids 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 Halleyids 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 to existing observations.

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 behaviour as its basis, millions of test particles are ejected from the comet nucleus and numerically integrated forward in time. The characteristics of Earthimpacting 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: 1. the date and duration of potential meteor activity on Earth 2. the characteristics (age, size, etc.) of the trails involved 3. the stream's structure close to Earth 4. the radiants of Earth-intercepting meteors 5. the approximate shower intensity The purpose of this section is to summarize the main model parameters chosen for our Halleyids meteoroid simulations. For additional information about the model, the reader is referred to Egal et al. (2019).

Nucleus
In our model we take 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 2 . We assume a nucleus density of 0.3 g cm −3 , which matches that of Rickman (1986). The geometric albedo of the nucleus is taken as 0.04 (cf. Section 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 Section 5.3. As in Egal et al. (2019), the ephemeris of each apparition of the comet is computed from the closest available orbital solution. Apparitions prior to 1986 CE are built from the set of osculating elements provided by Yeomans and Kiang (1981), converted into ecliptic J2000 coordinates. The comet coordinates since 1986 CE are determined from the JPL J863/77 solution 2 . The model parameters are summarized in Table 1.

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 velocities model of Crifo and Rodionov (1997). Crifo and Rodionov (1997)'s model is restricted to the sublimation of These bins approximately correspond to meteors detectable by radar (bin 1), visual/video means (bin 2) and as fireballs (bin 3).
Despite the unreliability of the comet's location prior to 1404 BCE (Yeomans and 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 and 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 is assigned a weight representing the 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 Section 8.1.1.

Integration and data analysis
The comet and meteoroids are integrated with a 15 th order RADAU algorithm (Everhart 1985), with a precision control parameter LL of 12 and an external time step of one day. The integration is performed considering the gravitational influence of the Sun, the Moon, and the eight planets of the solar system. General relativistic corrections and non-gravitational 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ý and 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 well described with heliocentric elements only 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 induce oscillations of the comet's heliocentric elements when computed far from perihelion. These oscillations are an artifact of the non-inertial heliocentric frame and 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 simulations analysis. The terms "impactors" and "impacting particles" are applied to simulated particles approaching Earth within the distance ∆X and time ∆T defined in Section 7. Physically, these represent test particles whose actual position in space best represent meteoroids which 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.

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 one thousand particles at each apparition of the comet between 1404 BCE and 2050 CE as described in Section 5. Each particle's position, orbit and node location is recorded at a regular time step of 10 years 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 mm, 1-10 mm and 10-100 mm.
6.1. Current structure Figure 4 presents the current location of the simulated meteoroid stream, including all meteoroids having nodal crossings between 1990 CE and 2020 CE to enhance its visibility using this reconnaissance survey. The theoretical shell model presented in McIntosh and Hajduk (1983) and McIntosh and Jones (1988) is reproduced at the top of Figure 4 for comparison.
The simulated meteoroid stream intersects the Earth's orbit at two locations, and reproduces the ribbon-like shape predicted by McIntosh and 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 behaviour of the meteoroid stream is in good agreement with the McIntosh and Hajduk (1983) and McIntosh and 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 between 194°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 as well reproduced. The Orionids seem to be composed of a larger percentage of older particles, ejected before 600 BCE (see Figure 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 Figure 5, we see that particles ejected between 1404 BCE and 600 BCE contribute to the present Orionid meteor shower. Simulations of McIntosh and Jones (1988) and Ryabova (2003) implied that meteoroids ejected after 1404 BCE and of mass higher than 1 g and 0.001 g respectively do not approach Earth's orbit in recent times. To compare our model with these results, we present in Figure 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 our simulated stream does not reproduce the shower's observed duration. However, and in contrast with Ryabova (2003), we find that particles 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 Figures 5 and 6, we see that in the case of the η-Aquariids, the simulations predict the onset of the shower at low solar longitudes 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 600 -610 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 solar longitude 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. meteoroid stream, but the existence of such a gradient is hard to assess because of the overlapping locations of the trail's nodes.

Time evolution
To illustrate the evolution through time of the simulated streams, we independently analyze 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 Figure 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 and Jones (1988), the nodes of the ejected trail (in green) quickly spread, exceeding 2 AU in extent within 3000 years. The median location of the nodes from the 1404 BCE ejection at each epoch is indicated with black dots in Figure 7. The median node location of the 1404 BCE trail has clearly 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 Eta Aquariid activity roughly from 0 -1000 CE.
In contrast, the median ascending node of the 1404 BCE trail in Figure 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 behaves 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 mm to 1 mm particles. The median node trajectories are in agreement with the 2020 CE nodal footprints observed in Figure 5, confirming the picture that some particles remain on orbits which 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.

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 Fig. 7. The evolution as a function of time (color coded by ejection year) of the descending nodal footprint representing the Eta 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 mm and 1 mm. The Earth's orbit in 2020 CE, computed in the barycentric frame, is represented with a blue line. somewhat different from the locations at which they actually cross the ecliptic plane (their "nodal-crossing positions" as defined in Section 5.3). In the previous section, we found that trails ejected prior to 600 CE and 600 BCE could contribute to the present η-Aquariids and Orionids 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 estimate the time T S required for a newly ejected trail to entirely spread along the orbit. After ejection, we identify the average orbital period (µ) and standard deviation (σ) of the simulated particles and consider two particles of periods P 1 = µ + 5σ and P 2 = µ − 5σ. After n revolutions, particle 1 lags particle 2 by ∆T = n(P 1 − P 2 ). We therefore approximate T S to the time required to satisfy ∆T = P 1 , i.e. T S = P 1 P 2 /(P 1 − P 2 ). In this formulation the influence of radiation and gravitational perturbations, that could accelerate the particles' spread, is not taken into account.
A&A proofs: manuscript no. Egal2020c For several trails ejected between 1404 BCE and 2050 CE, we estimate an average T S time of 310 years for 0.1 to 1 mm particles, 640 years for 1 to 10 mm particles and 1400 years for particles of 10 to 100 mm in size. The dependence of T S with 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 close to 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 and 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 nodalcrossing locations for trails ejected before 600 CE. However this approach might not be valid to analyze the fine structure of the stream, which is conditioned by planetary perturbations.

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 which can intersect the Earth. Since many of our test particles are from very old ejections, we don't 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 better reproduce year to year variations in activity as this approach captures features such as meteoroids residing in mean motion resonances. We next describe the details of how we declare a simulated particle an impactor.

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 a time parameter depending on the accuracy of the shower prediction. For this study, ∆T was fixed at 12h, leading to a distance criterion of 1.88 × 10 −2 AU for impactor selection.

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), i.e. around 800 BCE for the Orionids and 500 CE for the η-Aquariids. Only 222 of the ∼4.5 million particles simulated 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 5 particles met the criteria of a 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 (a factor of two increase), 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 Figure 9. Despite the increase in ∆X, most of the η-Aquariids and Orionids annual apparitions in this period contain less than 10 impacting particles. Due to the small number of particles retained and the lack of realistic weights applied to the simulations (cf. Section 8.1.1), again no reliable intensity prediction can be generated from Figure 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. Figure 9 shows impactors between 1980 and 2020. To zeroth order, the larger the number of model impactors in a given year, the higher we expect η-Aquariid and Orionid activity. 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 Figure 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.

Validation
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 and Watanabe (2007), Rendtel (2008), Sato and Watanabe (2014) and Sekhar and 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.

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 0 by a swarm of particles on the same orbit, centred around p 0 , and crossing the ecliptic plane at p 0 '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 less than 0.05 AU from Earth's orbit over an entire year (i.e., ± 6 months either side of the time when the Earth intersects the shower node). With this relaxed selection threshold, our number statistics improve dramatically.
We find 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 (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 at http://astro.uwo.ca/~wiegert/Halleyids/ Halleyids_animated_nodes.mp4. 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 sug-gests 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. Section 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 val-ues 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 of number statistics using different methods are required to produce reliable predictions of Orionid profiles and peak intensity in other years.

Postdictions
The final step in our modeling is to calibrate our model and determine particle weighting for the Halleyids meteoroid streams through direct comparison with meteor observations. We first compute synthetic intensity profiles (model shower activity as a function 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.

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 solar longitude, corresponding to the date when the Earth is the closest to the particle's node. Within a given solar longitude bin, the total number of particles is converted into an equivalent flux F , which is then transformed into a ZHR by Koschack and Rendtel (1990)'s relation: where A s is the typical surface area for meteor detection by a visual observer in the atmosphere at the ablation altitude (A s ∼ 37200 km 2 ) and r is the measured differential population index. We assume 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 Section 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: 1. The initial number of particles of a given radius ejected at a given epoch, 2. 1P/Halley's dust production at the heliocentric distance of ejection, measured via Afρ, 3. The differential size frequency distribution of the particles radii after ejection, characterized by the differential size index u.
In practice, the final weight is normalized by the number of particles sharing 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 which is used. For additional details about the weighting scheme, the reader is referred to 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 η-Aquariids analysis). For cases with fewer simulated particles, additional empirical weights were 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 to obtain sufficient statistics. The three complementary weights applied to our simulations are: 1. A coefficient C 1 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 that gives more importance to nodes close to Earth's orbit. This weight is an approximation of the impact parameter described in Moser and Cooke (2008). 2. A coefficient C 2 proportional to 1/n, with n the number of revolutions completed by the trail since its ejection. This parameter, similar to the f M 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 and Sato (2008) to model the Draconid meteor shower. 3. A coefficient which is a 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 years leads to no simulated particles available for certain Orionids apparitions. so this parameter is increased to 1.5 years in some case (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 γ δ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 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 found 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 and 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 2 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 Halleyids meteoroid streams self-consistently, we apply the exactly 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 the two showers as expected to be the case on purely physical grounds.

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 the Canadian Meteor Orbit Radar (CMOR), operational since 2002 ). 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 Orionids simulations.
The determination of the weights was performed in a semiautomated way. We first estimated the value of the particles' size distribution index u using a Particle Swarm Optimization (PSO) algorithm (Kennedy and 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 found 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 Section 2.3). On the other hand, we found 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 was 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 η-Aquariids and Orionids 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 1 , C 2 , and γ) was performed manually in a second 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 is divided by the approximate number of revolutions it has performed since its ejection n. This operation was found to improve the average level of the η-Aquariids 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.

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 Section 5.2. The inclusion of trails ejected before the 1404 BCE apparition, that helps improve our results in some cases, is discussed in Section 10.

Average activity profile
The average activity profile of the simulated η-Aquariid and Orionid apparitions between 2002 and 2019 is presented in Figure 12. The average profiles as measured with the CMOR 29 MHz and 38 MHz systems, with the IMO Video Meteor Network (VMN) and visual observations are plotted for comparison. The simulated activity profiles were scaled to the observations of the η-Aquariids shower, and the same scaling factor was applied to the modeled Orionids profile.  Egal et al. (2020). The best-weighted simulated average profile using one degree solar longitude binning is drawn in red. Note the factor of two difference in the ordinate scale between the top and bottom plots.
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 of 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 premaximum shows some divergence at solar longitudes 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, i.e. about a factor 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 suggest that most of the η-Aquariids activity can be reproduced with trails ejected since 3000 BCE. On the other hand, as expected from the discussion in previous sections, the Orionids average profile can not be modeled completely with the same trails, and require 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.

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 Figure 13. Only particles contributing to the main peak of activity of the η-Aquariids (44°-46°SL) and Orionids (208°-2010°S L 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.
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 which 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.

Annual variation
The ZHR of the main peak of each simulated η-Aquariid and Orionid apparition between 1985 and 2020 is presented in Figure 14. The showers' individual activity profiles between 2002 and 2019 are provided in Appendix C. As noted in Section 8.2.1, Fig. 13. The 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. 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 Figure 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 section 8.2.1 is about 2.5 is not respected anymore. However, we find this scaling necessary to compare the predicted and observed activity profiles shapes and maximum rates with clarity and is likely due to the model not including older ejections which make up proportionately more of the Orionids.  Egal et al. (2020). The envelope of estimated simulated peak variations, with an arbitrary uncertainty of ±20 meteors per hour allowing to frame most Orionids 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.   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 to raise our models' 2003 and 2014 ZHR by a factor of about 1.5. Since the purpose of this work is to predict future Halleyids 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 but 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 Orionids main peak intensity prior to 1996, around the 2006 and 2007 outbursts, and since 2010. ZHR values at other times may be underestimated, with the exception of 1998 apparitions that is significantly overestimated. In the IMO Visual Meteor DataBase (VMDB), we find the observed maximum Orionid rate in 1997 reached 27 meteors per hour 3 , 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 Figure 14). We therefore consider this model successful in reproducing the global yearly activity variations of the η-Aquariids and Orionids.

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 Orionids profiles, especially for solar longitudes 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 year periodicity of the Orionids main peak location  can not be reliably investigated with our simulations. The pre-maximum branch of activity is in general in good agreement with the observations, with some notable exceptions (e.g., in 2002-2003, 2008 and 2014), suggesting it is relatively composed of 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, Orionids in 2006 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 well the intensity, duration and shape of the observed activity. 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 and Watanabe (2007) and Sekhar and Asher (2014)   The average shower radiant location in one degree solar longitude 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 . Simulated radiants computed for particles ejected before and after 1450 BCE are plotted with cyan and dark blue circles, respectively. 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 well the radiant structure around CAMO measurements, 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 meteoroids 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).

Forecast
The comparison of our new Halleyids 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 and 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 are 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 the η-Aquariids and Orionids characteristics 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.

Annual variation
The main peak ZHR of each η-Aquariid and Orionid apparition between 2020 and 2050 based on our calibrated model is shown in Figure 16. Consistent with Sekhar and Asher (2014), we expect no strong Orionid activity over the analyzed period. The predicted meteor rates vary between 20 and 40 meteors per hour, corresponding to the usual range of ZHR reported for the shower . In contrast, Figure 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.

Outbursts characteristics
Details of the predicted future η-Aquariids outbursts are presented in Figure 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°F 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 do 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. 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 solar longitude for a given year and the model cumulative mass distribution. The middle trail plots show each ejections relative contribution by circle size (small to large) and color following the scale bar with black showing trails with lowest contribution to yellow trails which provide 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 Figure 13 with estimates of the differential mass index s in the visual (red), video (green) and radar (blue) mass ranges provided for comparison.

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 ZHR 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 and Kiang 1981). To investigate the influence of 1P/Halley's ephemeris prior to 1404 BCE in our simulations, we have ejected 5000 particles of 1-10 mm in size from different clones of the comet  Figures  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 years in the case of the Orionids, and 10.72 years for the η-Aquariids. The relative contribution 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. lution to the model for each shower is presented in red in Figure  18.
The best PSO solution for the peak Orionid activity has a period of 11.83 years, close to the value 11.88 found in Egal (2020) and Jupiter's orbital period of 11.86 years. 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 Figure 18. There is a clear correlation between the activity level of the shower and the number of particles reaching the Earth in a given year which were ejected when 1P/Halley was in the 1:6 MMR. This suggests that the 1:6 mean motion resonance induces a ∼12-year cyclic variation of Orionid activity.
The annual η-Aquariids activity is better described with a sinusoid having a period of 10.72 years. The sinusoidal signal fit is in good agreement with the modeled peak ZHR since 2005, but does not reproduce well earlier measurements of the shower. When examining the possible contribution of 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 η-Aquariids annual activity is suggested by our analysis but not as convincing as the Orionid fits. This periodicity needs to be confirmed by future radar and optical observations of the shower.

Conclusions
We have developed a new model of the meteoroid streams released by comet 1P/Halley and used it to characterise 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 have 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 years of Halleyid observations, self-consistently applying a unique weighting scheme to the two showers.
From the model we find that the η-Aquariids 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 Orionids 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 have 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 −5 to 10 −7 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 years, while the Orionids are older. Future simulation efforts could better constrain the upper range to the Orionids 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 years.
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 into the 1:6 MMR resonance, similar to several previously observed η-Aquariid and Orionid outbursts.
Among the mean motion resonances with Jupiter and Saturn that can shape the long term evolution of the Halleyids meteoroids 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 and Watanabe 2007;Rendtel 2008;Sato and Watanabe 2014;Sekhar and 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 years. A potential periodicity of 10.7 years for the η-Aquariids needs to be confirmed by future observations of the shower. Visual, video and radar observations of each η-Aquariid and Orionid apparition are strongly encouraged to provide additional constraints on numerical models.