© ESO 2019

1 Introduction

The surfaces of periodic comets in the orbits with a perihelion in the region of terrestrial planets are periodically heated by the Sun when they are close to perihelion, and are cooled when they move near aphelion. The heating causes a release of volatile gases, and the gas jets lift the dusty particles.

Since the ejection speed of a meteoroid is typically a few orders of magnitude lower than the heliocentric speed of the comet nucleus from which it originates, and its radius vector and that of the parent comet at the moment of ejection are practically identical, a meteoroid initially moves in an orbit that is very similar to the orbit of its parent body. We can say that all this matter moves in a common orbital corridor.

The non-gravitational effects and perturbations of planets, especially if the meteoroid orbits are situated near the resonances with the planets, can move the initial orbits of a subset of meteoroids to another orbital corridor or corridors. If there are several corridors and the Earth’s orbit passes through them, two or more meteor showers originating from the same parent body can be detected. Sometimes, there are many substreams, and the whole stream of the particular parent body can be regarded as the meteoroid complex associated with the parent.

Among the parent bodies of known meteoroid streams, comet 2P/Encke is the parent body of a very structured and extensive meteoroid complex. It was thought to be the parent of the Taurid complex (Whipple 1940). However, as we discuss below, there is no unique relationship between the comet and this complex. The researchers have thus tried to find another origin of the Taurid complex. Steel et al. (1991a) investigated the origin from a giant comet, a progenitor of the complex and 2P itself, which entered the inner solar system some time in the past 10 000−20 000 yr. It appeared, however, that this origin is possible only if high velocities (above those feasible in conventional ejection) of the particles escaping from the comet nucleus are assumed. Steel et al. (1991b) recognized that the Taurid complex is a continuous complex observed at four different crossings of the Earth’s orbit. According to these authors, comet 2P as the parent body of the whole complex is problematic because of a large difference of meteoroid velocity at perihelion, up to ~ 3 km s⁻¹. As an explanation, they suggested a splitting of large fragments originating from the original comet.

According to the Taurid swarm theory (Asher & Clube 1993; Asher 1994; Asher et al. 1994), Asher & Izumi (1998) found observational evidence that the long-term (~ 10 kyr) action of a mean-motion resonance with Jupiter can produce structure in a meteoroid stream in the form of a dense swarm, i.e., a non-uniform distribution of meteoroids along the orbit. Their comparison of predicted enhanced meteor and fireball activity from a Taurid swarm in the 7:2 resonance with the observational data collected in Japan over several decades showed a discernible correlation.

The additional parent bodies of Taurids were also searched for among the asteroids in the orbits approaching the Earth’s orbit. Asher et al. (1993) showed that a statistically significant number of Earth-crossing asteroids are part of the Taurid complex. They also identified another group of objects in the adjacent orbital phase space aligned with asteroid (2212) Hephaistos and suggested that the group could originate from the same progenitor as the Taurid complex or, alternatively, the Hephaistos and Taurid complex groups are similar because of the existence of a preferred “entry corridor” (Napier 1983; Wetherill 1991) for the capture into sub-Jovian orbits of Jupiter-family comets.

Babadzhanov (1999) proposed a set of asteroids that contributed to the complex. The set was then corrected (Babadzhanov 2001; Babadzhanov et al. 2008) because the original criteria for their relationship were found to be imperfect. However, we cannot be sure that the new, again subjective, criteria are fully satisfactory. This circumstance means that some of these proposed asteroids are still uncertain parent bodies of the Taurid complex. Babadzhanov et al. (2008) further concluded that some of the asteroids are related, with Northern and Southern Taurids and also with β-Taurids and ζ-Perseids. Porubčan et al. (2006) pointed out a possible association between seven Taurid filaments and nine near-Earth asteroids; they regarded the Southern Piscids and o-Orionids as Taurid filaments as well.

An important aspect in the dynamics of Taurid meteoroids was noticed by Asher (1994). The orbital period of these meteoroids can be close to 3.39 yr. If this is the case, then this period is very close to 2/7 of the orbital period of Jupiter (close to 11.86 yr) and the perturbations by this planet act in such a way that the long-term average value of the period is exactly 3.39 yr. The formation of swarms of meteoroids in the Taurid complex can be explained in principle with the help of the 3:1, 4:1, 10:3, and 11:3 mean-motion resonances with Jupiter (Asher 1994). The swarms occur due to non-uniform distribution of the meteoroids in the mean anomaly.

The difference between the orbits of comets 2P and the Taurids, noticed at the end of the 20th century, is a problem that could perhaps be explained via numerical modeling and following a further dynamical evolution of the 2P stream. In this work, we present the modeling of the 2P stream and follow its dynamical evolution. The specific aim of our modeling is to predict the maximum number of filaments of the 2P meteoroid complex that cross the Earth’s orbit. Also, we want to identify the meteor showers predicted in this way with their real counterparts as recorded in several databases considered. We also identify the predicted showers with those in the IAU-MDC list of all showers whose mean orbits are well-known. We note that the extent of the databases of meteor orbits has considerably increased in recent decades, and this circumstance provides us with a better chance than in the second half of the 20th century to find or disprove the agreement between theory and observation.

Our study should mainly answer two questions: Is comet 2P the parent of meteor showers and what are these showers? Are the Southern and Northern Taurids the main parts of the 2P meteoroid complex, or is there only a rough relationship between the Taurids and 2P? Since the stream of 2P alone is largely structured, we confine our study exclusively to the complex of 2P, i.e., we do not deal with the whole Taurid complex and its further possible parent bodies in this work.

Before starting the description of our work let us clearly define the sense of our usage of several basical terms. We use the words “complex” and “stream” as synonyms referring to the whole structure of meteoroids released from the parent body. The complex often consists of several “filaments”, i.e., the groups of meteoroids. The particles of each filament orbit the Sun in a single, common corridor. If some particles of the filament hit the Earth, they can cause a meteor “shower” corresponding to this filament.

2 Nominal orbit of comet 2P/Encke and its evolution

In our study, we consider the orbit of comet 2P/Encke with the orbital elements published in the JPL small-body browser (Giorgini et al. 1996)¹. The elements of this orbit and the non-gravitational parameters are listed in Table 1. Hereafter, the orbit is referred to as the nominal orbit.

To evaluate the determination uncertainty of the nominal orbit of 2P, we created 100 clones using the method by Chernitsov et al. (1998). In this method the elements of nominal orbit are written in the form of covariant 6 × 1 matrix y_o. The equallystructured covariant matrix with the elements of the orbit of jth clone, y_j, can be calculated as (1)

where A is a triangle matrix such that the product AA^T equals the covariance matrix related to the process of nominal-orbit determination¹. Since we have six elements, A is the 6 × 6 matrix; η is the 1 × 6 matrix with each element being a value randomly chosen from a Gaussian distribution peaked in zero and having the dispersion σ = 1; η^T is the transposed matrix η.

To numerically integrate the orbit of the parent comet and its clones together with the perturbing planets, we use integrator RA15 (Everhart 1985) within the software package MERCURY (Chambers 1999). The gravitational perturbations of the eight planets, Mercury to Neptune, are taken into account. We consider the heliocentric equatorial coordinate frame to perform all integrations.

The evolution of the four orbital elements of the nominal orbit and the evolution of the orbits of the clones are shown in Fig. 1. If the comet was influenced only by the gravity of the Sun and the planets, its nominal orbit has been reliable during the last millennium. However, there are no large deviations between the nominal orbit and orbits of clones, with a single exception, either during the whole followed period of 50 kyr.

An uncertainty of the nominal orbit due to the non-gravitational jet effect is estimated considering the known values of non-gravitational parameters. The acceleration of a comet due to the jet effect is (Marsden et al. 1973) (2)

where A_i are the non-gravitational parameters determined from the observations and function g(r) was found in the form (3)

where a = − 2.15, b = 5.093, h = −4.6142, r_o = 2.808 au, and α_ng = 0.1113.

The acceleration component a₁ is oriented in the direction of the heliocentric radius-vector of the comet; a₂ is perpendicular to a₁, lies in the comet’s orbital plane, and is oriented in the sense of the comet’s heliocentric motion; and component a₃ is perpendicular to both a₁ and a₂, whereby system a₁ − a₂ − a₃ is right-handed (the O(rtp) coordinate frame).

We use integrator RA15, which assumes rectangular components of acceleration in the coordinate frame O(xyz). The acceleration vector due to the jet effect is, however, determined in the coordinate frame O(rtp) (having the same origin as O(xyz)). Therefore, we need to transform the components of the acceleration due to the jet effect from O(rtp) to O(xyz). This can be done by using the common transformation mathematics, and the resultant components of the acceleration vector in O(xyz) are (4) (5) (6)

where r = (x,  y,  z) and v = (v_x,  v_y,  v_z) are radius and velocity vectors, respectively, of the particle in the coordinate frame O(xyz), and the components ṽ_y and ṽ _z of auxiliary velocity vector ṽ = (ṽ_x,  ṽ_y,  ṽ_z) are (7) (8)

Vectors r and v in O(xyz) are known in each integration step; therefore, the components of acceleration due to the jet effect can be calculated using Eqs. (4)−(8), and these components are added to the corresponding components of the gravitational acceleration due to all massive bodies (the Sun and perturbing planets) resulting from the RA15 algorithm.

The values of non-gravitational parameters, A₁ to A₃, can be found on the web pages of the JPL small-body browser and are also given in Table 1. We denote the determination errors of parameters A₁ and A₂ as ΔA₁ and ΔA₂. We construct another set of the clones of nominal orbit (orbit with A₁ = A₂ = A₃ = 0) with the non-zero values of A₁, A₂, and A₃. Since the determination error of A₃ is unknown, we do not vary this parameter. Parameters A₁ and A₂ acquire, in the cloned orbits, all combinations of values A₁ − ΔA₁, A₁, and A₁ + ΔA₁ for the first parameter, and A₂ − ΔA₂, A₂, and A₂ + ΔA₂ for the second parameter, and so there are nine combinations.

The evolution of perihelion distance, semi-major axis, eccentricity, and inclination of the comet’s nominal orbit and orbits of the second set of clones (which are influenced by the non-gravitational effects) is shown in Fig. 2. In this figure, we can see that the effects start to be significant in a farther past than about five millennia. And, they do not cause any radical change until 50 kyr in the past. The nominal orbit still appears to be one of the most appropriate orbits to model the stream of 2P.

The minimum distance between the orbits of the comet and Earth, often referred to as the minimum orbit intersection distance (MOID), is also investigated. The behavior of this distance is shown in Fig. 3 for the period of the last 50 kyr. We considerthe orbit influenced only by the gravity of the Sun and planets and the orbit that was also influenced by the jet effect. Both the post-perihelion and pre-perihelion arcs of the comet’s orbit, whether it was or was not influenced by the effect, have periodically approached the Earth’s orbit at an extremely short distance, often shorter than 0.01 au. At the present, the post-perihelion (pre-perihelion) arc is at a distance of 0.17 au (0.19 au) from the orbit of our planet. It means that the meteoroids that were released very recently from the comet surface will not collide with the Earth at the present time.

The evolution of the orbital nodes of 2P, when only the gravitational perturbations are considered, is shown in Fig. 4. If we also take into account the jet effect, described with the currently valid non-gravitational parameters, then no significant difference can be detected. Both ascending and descending nodes of the 2P orbit have crossed the Earth’s orbit many times during the past several 10 000 yr. The crossing points are distributed along the entire orbit of our planet. This is the reason why the complex of 2P is expected to have a largely extended and complicated structure that can collide with our planet during a relatively long time.

Fig. 1 Behavior of perihelion distance (panel a), semi-major axis (panel b), eccentricity (panel c), and inclination to the ecliptic (panel d) of the initially nominal orbit of comet 2P/Encke (red curve) and the orbits of the first set of its clones (green curves) mapping the uncertainty of the determination of the nominal orbit. The evolution is reconstructed backward for 50 000 yr. A non-gravitational force is ignored.

Fig. 2 Behavior of perihelion distance (panel a), semi-major axis (panel b), eccentricity (panel c), and inclination to the ecliptic (panel d) of the initially nominal orbit of comet 2P/Encke, which is (violet curve) or is not (red curve) influenced by the non-gravitational force (jet effect). The evolution of the corresponding elements of the orbits of the second set of its clones (green curves), which delimit the uncertainty of non-gravitational parameters A1 and A2, are also shown. The evolution is reconstructed backward for 50 000 yr.

Fig. 3 Evolution of the minimum distance between the orbital arcs of comet 2P/Encke and Earth from 50 000 yr before the present up to the present. The minimum distance of the post-perihelion (pre-perihelion) arc is shown by the violet (cyan) curve when the evolution of the 2P orbit was followed ignoring the non-gravitational force. When this force is taken into account, the minimum distance is shown by the green (brown) curve.

Fig. 4 Positions of the orbital nodes of comet 2P/Encke during the last 50 000 yr. The green circle indicates the orbit of the Earth. The red (blue) curve shows the positions of ascending (descending) node.

3 Stream modeling

In reality, the individual meteoroids of the 2P stream are obviously released from the surface of this parent comet at various times, wherever in its orbit (not only at perihelion), and have various sizes and shapes. In addition, there can be a difference in their chemical compositions and internal structures. All these detailsare largely unknown. It is therefore impossible to create a completely realistic model. Hence, we create a set of models mapping some parts of the real stream from the dynamical aspect.

In addition, there is also an uncertainty in the past evolution of the comet’s orbit. Since the orbital period of 2P is the shortest of all known comets, it passes its perihelion most frequently and this circumstance implies a large net effect from non-gravitational forces. Since the non-gravitational forces are stochastic, their influence on the dynamical evolution of the comet is not predictable and the parameters of its orbit are rather uncertain.

Becauseof a relatively large variety of the modifications of the comet’s orbit in the past, it is practically impossible to study all evolutionary possibilities. Hence, we start our study of the meteoroid stream of this comet considering its nominal (catalog) orbit. In themodeling, we do not consider the non-gravitational force (jet effect) acting on the comet since this force is variable, and thus unknown on the timescale of the constructed models.

A specific part of the theoretical stream of 2P is modeled in the way that was first used by Neslušan (1999). We slightly modified the procedure later on (Tomko & Neslušan 2012). Since the modeling has been described several times (e.g., Hajduková & Neslušan 2017), we only briefly recall its main steps here.

At first, the considered orbit of the parent body is integrated backwards in time to the moment of the body’s perihelion passage that happened closest to an arbitrarily chosen time t_ev. In more detail, t_ev is regarded as the evolutionary period (of stream particles) accounted from the past to the present; its value is positive. In the perihelion, a cloud of 10 000 test particles, representing the meteoroids, is modeled around the parent body. To model the cloud we simply assume an ejection of the particles in all directions uniformly and with the same ejection speed, v_ej, which equals 10⁻³v_q, where v_q is the heliocentric speed of the parent comet in the perihelion, i.e., at the moment of modeling. The values of the ejection speed were v_ej = 70.3, 73.4, 75.3, 82.6, and 92.8 m s⁻¹ in the models for t_ev = 1, 2, 4, 8, and 16 kyr, respectively.

A given specific part of the stream consists of particles ejected only in a single passage through the perihelion of the parent. Of course, the real meteoroids are not ejected only in the perihelion and not uniformly in all directions and all with the same speed. Our modeling does not aim to perfectly follow the real process of the ejection of meteoroids. We only aim to fill in the appropriate orbital phase space with the particles. At this filling, we assume that the dynamical characteristics of their initial motion, when they were ejected, are forgotten after a relaxation time. Our experiments with a non-uniform release of particles (asymmetric release with various speeds) indicate that after about 30 orbital revolutions of the parent body (~ 100 yr in the case of 2P), the particles move in practically the same stream regardless of the way they were released from the parent body.

After the ejection of the particles in a given stream part, we integrate their motion, as well as the motion of the parent comet and perturbing planets, in forward time. The integration is ended in the present time. At the end of the integration, the particles that move in the orbits that approach the Earth’s orbit within 0.05 au are selected (i.e., their MOID is smaller than or equal to 0.05 au). The characteristics of these particles correspond to those of the expected shower, and thus they enable us to predict the mean characteristics of the shower. As stated below, we create a set of models, each for a given combination of input parameters. A real shower can be predicted with the help of more than one of these models.

In our modeling, the numerical integration of all particles and gravitating massive bodies is performed, again by using integrator RA15 (Everhart 1985) within the software package MERCURY (Chambers 1999). The gravitational perturbations of the eight planets, Mercury to Neptune, are taken into account.

In the models, the dynamics of the particles is assumed to be influenced by the Poynting–Robertson (P–R) effect. The term P–R effect is used here to refer to the action of radial electromagnetic radiation pressure and to the effect of solar wind on the meteoroid particles. Specifically, we use the formula to calculate the acceleration of assumed spherically symmetric particle due to the P–R effect derived by Klačka (2014). When the dynamics of an object in the solar system is calculated where only some average properties of local solar wind are well known, Klačka recommends using formula (68) in his paper: (9)

In this relation, β is the ratio of the non-gravitational force due to the P–R effect (i.e., force due to the solar electromagnetic-radiation and solar-wind pressures) to the gravitational force by the Sun, G is the gravitational constant, M_⊙ is the mass of the Sun, r is the heliocentric distance of the meteoroid particle, andv = (v_r,  v_t,  v_p) is its heliocentric velocity in the O(rtp) coordinate system with the components radial (v_r), transverse (v_t), and perpendicular (v_p). Component v_t is assumed to lie in the particle’s orbital plane and v_p is perpendicular to both radial and transverse directions. Also, c is the speed of light in a vacuum, is the dimensionless efficiency factor of the radiation pressure averaged over the solar spectrum, and η₁ and η₂ are the parameters characterizing the solar wind. In Eq. (9), quantity u is the magnitude of the velocity of the individual solar wind particle with respect to the Sun. If the solar wind is not time-variable (as we assume in our calculations), term can be neglected. Vector e_R is the unit vector in the direction of the coordinate r-axis. Its components are e_R = (1,  0,  0).

We use the denotation dv∕dt = (a_r,  a_t,  a_p). With its help, vectorial Eq. (9) can be re-written to the formulas giving explicitly the components of the acceleration due to the P–R effect, (10) (11)

and the perpendicular component a_p ≐0 for the spherically symmetric particle. Since we consider the acceleration due to the P–R effect, the last term in Eq. (9), (−GM_⊙∕r²)e_R, giving the gravitational acceleration by the Sun, is omitted. We neglected term and calculated v ⋅e_R = v_r. According to Klačka (2014), the recommended values for η₁ and η₂ are η₁ = 1.1 and η₂ = 1.4 in the solar system. We also assumed theefficiency factor .

If η₁ = η₂ = 0 the action of solar wind is ignored, and Eqs. (10) and (11) are reduced to those giving the acceleration due to the solar radiation, which were well known a long time before the Klačka’s article was published. In this case, parameter β equals (12)

where L_⊙ is the solar luminosity and m is the mass of the meteoroid particle. Using the well-known values of quantities L_⊙, G, and M_⊙, the last relation can be given as (13)

Here, ρ is the mean mass density of meteoroid particle in kilograms per cubic meter and R is its average radius in microns.

Since we again use integrator RA15 giving the gravitational acceleration in the coordinate system O(xyz), we need to know the rectangular components of the acceleration due to the P–R effect in this coordinate frame. The acceleration vector due to the P–R effect is simply added to the vector of gravitational acceleration. The former vector, which is well known in the coordinate frame O(rtp), can be transformed to O(xyz) in the way already used and described in Sect. 2 in the case of the jet effect. Denoting dv/dt = (a_x,  a_y,  a_z) in the O(xyz), the individual components are (14) (15) (16)

where components ṽ_y and ṽ_z of auxiliary velocity vector ṽ = (ṽ_x,  ṽ_y,  ṽ_z) are again given by Eqs. (7) and (8).

To derive Eqs. (14)−(16), we utilized that v_p = 0. Vectors r and v in O(xyz) are again known in each integration step; therefore, the components of acceleration due to the P–R effect can be calculated using Eqs. (14)−(16) and these components are added to the corresponding components of the gravitational acceleration due to all massive bodies. In this integration, the jet effect is not considered.

The P–R effect parameter β depends on the properties of particles, such as the size, density, light scattering efficiency (then on albedo and light absorption ability), and mass. In the case of meteoroids, these properties are rather uncertain. So, we instead regard β as a free parameter. We consider a series of β values to search, as far as possible, for such a value, which results in a good match between the characteristics of predicted filaments and their observed counterparts, the meteor showers.

We create a series of models for various combinations of specific values of evolutionary time t_ev and parameter β. In reality, parameter β ranges over a wider interval of values since the stream consists of the particles of various sizes and densities. In addition, the particles arereleased at various times; therefore, their evolutionary times must be different and can also acquire the values from a wider interval. Hence, the particular model we created does not represent a whole stream. The model of a whole stream is a composition of partial models that give a good prediction of the corresponding real showers.

Specifically, we created the models for all combinations of values t_ev = 1, 2, 4, 8, and 16 kyr and β = 0.00001, 0.0001, 0.001, 0.003, 0.005, 0.007, 0.009, and 0.011.

4 Ways to identify predicted showers with their real counterparts

After the numerical integration of the particles in a given model was finished and the Earth-orbit approaching particles selected, we searched for a corresponding shower in the meteor databases. Specifically, we used one photographic, three video, and one radio meteor database. The details of all these databases are specified in Table 2.

Basically, we proceed in two ways to identify a predicted shower with its real counterpart: (1) we attempt to determine if there is a real counterpart to the predicted shower in each database considered, (2) we search for the real counterpart in the IAU-MDC list of all showers² (Jopek & Kaňuchová 2014), which contains the mean characteristics of the showers already found by other authors. In the following, we describe these ways in more detail.

To select these meteors from the given database, we used the “break-point method” (Neslušan et al. 1995, 2013). The method is based on an analysis of the dependence of the number of selected shower meteors on the limiting value of the Southworth-Hawkins (Southworth & Hawkins 1963) D_SH discriminant, D_lim³. Considering a value D_lim, we calculated, for every meteor in the data, the D_SH discriminant evaluating a similarity between the orbit of the meteor and mean orbit of the shower. If D_SH ≤ D_lim, the meteor is regarded as a member of the shower we investigated. After a scanning the whole database, we obtained N shower meteors selected at given value D_lim. Then, the orbits of these meteors were used to calculate a new mean orbit of the shower and the whole procedure, or iteration, is repeated until the difference between two consecutive mean orbits is negligible. The meteors separated by the last mean orbit are regarded as the members of the shower for the considered D_lim. Furthermore, this selection is repeated for a series of the values of D_lim to obtain the dependence N = N(D_lim).

The procedure described above within the break-point method can be performed by using an arbitrary initial orbit that enters the iteration. Hence, the procedure can sometimes provide us with a N = N(D_lim) dependence in the orbital phase space where no shower exists or is recorded in the data (even though we predicted it). If a real shower we search for is present in the data set, then the dependence N = N(D_lim) has a convex behavior with a constant or almost constant part, a plateau. The break-point is just the point at the beginning of the plateau (when we proceed from lower to higher values of D_lim). The value D_lim corresponding to the break point is the most suitable limiting value for the selection of the densest part of the real shower. We recall that a similar method, emphasizing a cutoff value chosen to reflect the strength of the shower compared to the local sporadic background, was also suggested by Moorhead (2016).

As mentioned above, in the break-point method the mean orbit of the real-shower meteors from the previous iteration step is used to go on with the iteration in the next step. At the proper beginning of the iteration, we have no previous orbit, however. The iteration to find the actual mean orbit of real shower starts instead with a theoretically predicted mean orbit. In our work, this initial orbit was determined as the mean orbit of a particular shower that was predicted in a given simulation or, for example, if we know the parent body and search for a shower consisting of meteoroids around its orbit, then just the orbit of the parent can be used as the initial orbit. In the case of the stream of comet 2P with a widely dispersed filaments, it was appropriate to change the break-point method in this step. Instead of the predicted mean orbit, we used the mean orbit of the most concentrated swarm of meteors in a given filament. Toselect these meteors, we divided the whole sky into small angular areas each delimited by the angle of 1° in right ascension and 1° in declination. Between these areas we then searched for those with the local maximum of surface number density, N, of the radiants of the predicted shower. This number density was calculated as the ratio of the number of radiants in the area and area’s space-angle size. In some models the maximum surface number density N of the radiants was relatively small. In such a case, we enlarged the size of both angles delimiting the particular area up to 3°.

The maximum found N was dividedinto eight levels. In Fig. 5, the areas with the number density n at various levels are plotted. For a given model, there can be several isolated zones of a higher concentration of radiants. Each concentration corresponds to the swarm of meteors. In a further processing, we considered only the swarms consisting of the areas with the maximum n obeying n ≥N∕4. For such a swarm, the mean orbit of the meteors with the radiants in these areas was calculated. In addition, just this mean orbit was the initial orbit entering the iteration in the break-point method.

Proceeding in this way, we found several concentrations of radiants and, consequently, of orbits in each model where the particles approached the Earth’s orbit at the end of simulation. A given concentration often occurred in several models where the corresponding mean orbit was not exactly the same; the mean orbit of the concentration was sometimes obtained with several modifications. When we subsequently separated a related real shower from a database, we considered all these modifications of the mean orbit, i.e., we repeated the separation of a given shower several times, for several initial orbits serving as the input into the iteration.

If the corresponding real shower was a well-defined compact shower, the break-point method would result in the same real shower after performing all the above-mentioned repetitions. Unfortunately, the separated real shower often appears to be structured and/or infrequent; therefore, we obtained several modifications, most probably of the same shower in these repeated separations, which led us to merge these modified showers into one single shower. The reason for the merging and how it was done is described in Sect. 5.2.

The mutual positions of predicted and observed radiants were compared, in some works, in the sun-centered ecliptical coordinates (see, e.g., Jenniskens 2017). The sun-centered ecliptical coordinates differ from the standard ecliptical coordinate system by the longitude, λ′, which can be calculated as λ′ = λ − λ_⊙, where λ is the common ecliptical longitude of the meteor radiant and λ_⊙ is the solar longitude of the meteor.

The symmetry of radiant areas of various filaments of the same stream used to be studied in the ecliptical coordinate system centered to the apex of the heliocentric motion of Earth. The longitude of meteor radiant in these coordinates, λ₂, equals λ₂ = λ − (λ_⊙ + 270°). Since there is a constant difference of 270° between λ′ and λ₂, the relative difference between the radiants is the same in both apex-centered and sun-centered coordinate systems. Therefore, we used the apex-centered ecliptical coordinate system to compare between the theoretical and real radiants and, at the same time, to study the symmetry of the radiant areas of various filaments.

The found real showers were also identified with the showers in the IAU-MDC list of all showers. Namely, some predicted showers might not be recorded in the databases considered, but these could have been found by other authors in other data on real meteors. In the process of this identification, we assumed that many real showers in the IAU-MDC list, which were not yet classified as the established showers, could be the actual showers, which could become established in future. Hence we did not omit the non-established showers in our analysis, though the established showers (included in the list) are regarded as a more reliable basis for the conclusions. Since only the mean characteristics of the showers are available in the list, we identified the theoretical filaments with the showers by using the D_SH discriminant in this case. We established an ad hoc criterion that the value of D_SH between the mean orbit of the filament and that of the recorded shower is equal to or smaller than 0.10.

Fig. 5 Panel a: number density of the theoretical radiants within the space angles of size 1° × 1°. With respect to the maximum number density, N, the number densities are divided into eight levels, corresponding to intervals delimited by 0, N∕8, 2N∕8, …, N; the seven most numerous intervals are color-coded according to the legend in the plot. The model for tev = 8 kyr and β = 0.007 is used as an example. Panel b: radiants of all Earth-approaching particles in the model.

5 Results

5.1 Predicted filaments of the 2P stream

The theoretical particles approached the Earth’s orbit within the chosen limit of 0.05 au in the models for all considered values of theevolutionary time, t_ev (see last paragraph of Sect. 3). However, the approach did not occur in every model, characterized with a given pair of values t_ev and β, because some values of β appeared inappropriate. No single approach was observed in the models for β = 0.00001, 0.0001, 0.001, and 0.003 and t_ev = 1 kyr. In addition, we omitted four models with the combinations of t_ev = 8 and 16 kyr and β = 0.009 and 0.011, because the radiants and the orbits of the approaching particles were widely dispersed. The radiants were situated in large parts of the northward and southward ecliptical hemispheres. The surface number density of the radiants was low, with no apparent maximum (in the models for t_ev = 16 kyr, the density was lower than 4 in whatever 1° × 1° space angle of sky). An attempt to distinguish between the filaments typically resulted in many predicted showers, but each with fewer than ten particles.

In the given model, for a given pair of t_ev and β values, the particles often approached the Earth’s orbit in more than a single filament; up to ten filaments in the model for t_ev = 16 kyr and β = 0.00001. In total, 87 filaments resulted from all the considered models together.

The geophysical characteristics of all these filaments are given in Table A.1 and their orbital characteristics in Table A.2. We give all these extensive data in the appendix since there are many real showers and their filaments in the adjacent orbital phase space and the result of our modeling (the data in the extensive tables) can help, in future, to recognize whether a particular shower is related to comet 2P. In addition to the geophysical and orbital characteristics of the predicted showers, the appendix also contains the orbital characteristics of the separated real showers that were identified as related to 2P (Table A.3). These data are supplementary to those in Table 3 which are discussed in the main text. The appendix also contains the found associations of the predicted filaments listed in Tables A.1 and A.2 and (i) the real showers separated from the databases used (Table A.4) as well as (ii) the real showers given in the IAU-MDC list of all showers (Table A.5).

In Tables A.1 and A.2, the individual filaments are distinguished with the working codes, each ending with the capital Roman numeral. The code together with t_ev and β enable the unique cross-identification of the particular filament in various tables. The filaments obtained from the various models are often almost identical (see Tables A.1 and A.2); therefore, many of the 87 filaments do not correspond to an independent shower, and thus fewer than 87 showers were predicted.

In the tables, the number of particles that approach the Earth’s orbit, N_s, is given for every filament. We recall that the total number of particles considered in each model is 10 000. One can see that the number varies from 13 to 1781; therefore, relatively weak to very strong showers are predicted.

The models for t_ev = 1 kyr (2, 4, 8, and 16 kyr) and all values of β resulted in the sum of 9 (14, 17, 21, and 26) Earth-orbit crossing filaments. All these filaments consisted of 4690, 13 364, 6090, 8912, and 6885 particles approaching the Earth for t_ev = 1, 2, 4, 8, and 16 kyr, respectively.It implies that the age of the largest number of the 2P stream particles at the Earth’s orbit is about two millennia.

The total number of filaments in all models (with all considered values of t_ev) for β = 0.00001 was 17. The number of Earth-orbit approaching particles in all these models was 5534. For β = 0.0001, 0.001, 0.003, 0.005, and 0.007, the numbers of filaments (Earth-orbit approaching particles) were 9, 8, 9, 13, and 15 (5534, 3998, 6747, 5152, 5308, and 5787), respectively, which means that there is practically no significant preference with respect to the size of particles among the particles constituting the showers associated with comet 2P.

Of the 87 predicted showers, 65 were identified with a corresponding real shower. Information about the identification of the particular predicted shower is given in the last column in Table A.1; several filaments were sometimes identified with the same real shower. Most of the 22 unidentified showers are (i) daytime showers, with the radiant area situated in close angular vicinity to the Sun on the sky at the maximum of their activity and/or (ii) weak, infrequent showers. Such showers are hardly detectable. In addition, most of the long-time evolving (models for t_ev = 16 kyr) and small-particle (with sizes corresponding to β ≥ 0.003) showers were also not identified. In all likelihood, the physical lifetime of the small particles in very short-period orbits, like that of comet 2P, typically does not reach 16 millennia.

More detailed information about the identification of the predicted showers with the real showers and descriptions of the found real showers are given in Sect. 5.2.

5.2 Identified real showers related to comet 2P

We said that using the break-point method resulted in the selection of a real counterpart to a majority of the predicted showers. The counterpart was often selected from several databases considered. Specifically, we selected 31, 31, 34, 28, and 51 real showers from the F, C, S, E, and R databases, respectively. However, the differences between the mean characteristics of some of these showers were within the interval of 1–σ of their determination errors. Therefore, these showers were obviously only the various modifications of the same shower. Therefore, we merged all these modifications from a given database to a single shower. After the merging, the resultant shower consisted of the set of meteors that were the unification of the meteors in all its modifications.

The iteration procedure yields the mean orbit, which is the average of the most concentrated set of the meteoroid orbits in given orbital phase space. If the initial orbit (predicted by a model) is an orbit at the border of a meteor concentration/shower or outside the border, the final mean orbit that results from the iteration can be significantly different from the initial orbit. Actually, the iteration sometimes provides a mean orbit of real shower which was, in the orbital phase space, far from the corresponding prediction. These showers were no longer considered.

As mentioned above, we merged the showers that were clearly similar (i.e., whose mean parameters were within the determination error). Since the meteoroid complex of 2P appeared to be active within a relatively long period during year, this merging could not yet provide the whole shower; a shower could still be divided into several parts with the different mean solar longitude, which were active in different time intervals. Because of this possibility, we determined the daily motion of the mean radiant of every separated shower from a given database. The equatorial coordinates of a radiant of a given meteor, α(λ) and δ(λ), corrected to the daily motion can be calculated as (17) (18)

where λ is the solar longitude of the given meteor and is the mean solar longitude of the shower or its analyzed part.

After the coefficients A_α, B_α, A_δ, and B_δ were determined, we illustrated the motion of the mean radiant of the showers separated from a given database. This was done separately for the northern and southern showers. In more detail, we plotted a path of the mean radiant on the sky during the period of shower activity (delimited by the meteors with the minimum and maximum solar longitudes). This illustrated path was an approximation of the actual motion of mean radiant on the sky since we approximated the behavior of each of the two radiant coordinates as the linear function of solar longitude. The illustration of the paths (see an example in Fig. 6) enabled us to see the motion of the mean radiant of all the showers we had at this step, and to evaluate a similarity.

In the example of the paths of mean radiant (Fig. 6), we noticed that the shower denoted C-N02 appeared to be a part of that denoted C-N01. So, these showers were finally unified and regarded as a single shower (C-N01 in tables). Showers C-N05 and C-N07 are regarded as other independent showers, and showers C-N03, C-N04, and C-N05 were much different from the originally predicted showers, and so they were no longer considered (C-N04 would be regarded as a part of C-N06, otherwise).

The candidate showers for the last merging, done on the basis of similarity of mean-radiant motion, were further analyzed if their differences in mean perihelion distance and mean eccentricity are smaller than the determination errors (1–σ). The showers with the differences that obeyed these limits were finally merged. Coefficients A_α, B_α, A_δ, and B_δ for the definitive showers are given in Table 3, and their mean orbital characteristics in Table A.3.

A general review of the real showers after the final merging of the parts can be seen in Table 4. The showers were separated from all the databases considered, i.e., F, C, S, E, and R (this letter is the first in the shower’s designation). The group of the separated showers constitutes the meteor complex of comet 2P/Encke. Altogether, it comprisesnine showers with the radiants on the northern sky and twelve showers with the radiants on the southern sky. Of these, one northern and four southern showers are new, not yet included in the IAU-MDC list of all showers.

The reliably determined showers in the IAU-MDC list are classified as “established”. Among the separated showers that are the parts of the 2P complex, only two northern and three southern showers belong to the established showers. Specifically, the northern showers are the Northern Taurids (#17) and Daytime ζ-Perseids (#172). The southern showers are the Southern Taurids (#2), Southern Daytime May Arietids (#156), and Daytime β-Taurids (#173). The numbers of meteors in the separated and predicted showers (see Table 4) imply that the Northern and the Southern Taurids are abundant in both predicted filaments and separated real showers. Showers #172 and #173 are predicted to be relatively abundant, but their real counterparts appear to be only diffuse, infrequent showers. This is likely a consequence of the fact that the radiants of these daytime showers are in the close vicinity of the Sun, and thus are hardly detectable. Shower #156 is infrequent. However, it was also predicted as an infrequent shower.

In addition to the established showers, the 2P complex also consists of showers #215, #629, #632, #634, #635, #726 (northern showers), and #485, #624, #626, #628, #636 (southern showers). In addition, we separated another northern (working denotation S-N11) and four southern showers (working denotations C-S09, S-S09, E-S05, and E-S08). We named these showers S-N11 ≡ η-Piscids, C-S09 ≡ γ-Geminids, S-S09 ≡ ω2-Aquariids, E-S05 ≡ 74 Orionids, and E-S08 ≡ 29 Piscids. In Table 4 we can see that showers #632, #635, #626, and #628 are predicted, as well as separated as the abundant showers. Hence, our study supports a suggestion that these showers should be regarded as actually real.

The meteor complex of 2P has a filamentary structure, which is reflected in the distribution of radiants of the related meteors. The positions of the radiant areas on the sky are shown in Fig. 7. In panel a, the predicted radiants of particles in the filaments identified with the corresponding showers in F, C, S, E, or R databases are shown in the equatorial coordinate system, which is often used for this purpose. In panel b, the radiants of real meteors in the corresponding showers are shown, again in the equatorial coordinate system. The real showers are distinguished with different symbols and colors. The same radiants in the apex-centered ecliptical coordinate system are shown in panels c and d: panel c corresponds to a, and panel d to b.

All showers of the 2P complex are the ecliptical showers since their radiant areas are situated near the ecliptic. The radiants are in all cardinal directions of ecliptic showers: helion and anti-helion, such that both have the northern and southern strands. Specifically, the radiant areas in the northern (southern) helion direction have Daytime ζ-Perseids (#172; Southern Daytime May Arietids (#156), and Daytime β-Taurids (#173)). The showers with the northern (southern) anti-helion radiants are Northern Taurids (#17), Northern δ-Piscids (#215), A2 Taurids (#629), November η-Taurids (#632), τ-Taurids (#634), A1 Taurids (#635), and December ɛ-Geminids (#726; Southern Taurids (#2), November ζ-Taurids (#485), ξ-Arietids (#624), λ-Cetids (#626), s Taurids (#628), and m-Taurids (#636)) according to our analysis.

Several showers were separated from more than a single database of real meteors. For example, the Northern Taurids (#17) were separated from the F, E, and R databases (Table 4). This was expected since a well-defined shower should be recorded in every comprehensive data set. Ideally, each shower should be recorded in every considered database. Inthe case of the 2P complex, the showers recorded in more than a single database agree with each other when they were separated from two (or more) of the F, C, S, and E databases. The centers of radiant areas were also situated near a common point on sky in this case. The differences in the individual parameters exceed 1–σ uncertainty when we compare the mean orbits of showers separated from R and another database.

A given shower can often be less numerous in one database than in other, in which case its radiant area is smaller. The example of the showers separated from four or three databases is shown in Fig. 8 in terms of the radiants. In panel a, the December ɛ-Geminids (#726), are separated from F, C, S, and E databases (colored symbols). Similarly, in panel b, the radiants of the s Taurids (#628), are shown. The latter was separated from the F, S, and R data sets (colored symbols, again).

In Fig. 8 the theoretical radiants of filaments, which correspond to showers #726 and #628, are also shown (various black symbols). In general, the prediction was not always perfect. In many cases, a difference like that in panel b occurred. The above-mentioned criteria for the identification of theoretical and real showers were, however, always obeyed.

The real showers R-N05, R-N06, and R-N10 are all related to the established shower #172 (see Table 4). The reason why three showers, which were separated as individual showers, are related to the same IAU-MDC shower is the following. The mean characteristics of shower #172 were determined independently by four authors (or author teams) on the basis of various observational data. Consequently, the characteristics of the shower as given by various authors are not exactly the same. This meant that the infrequent showers R-N05, R-N06, and R-N10 were not merged to a single shower in our separation, but could be related to the modifications of the same IAU-MDC shower.

As mentioned in the previous section, the mean orbits of the theoretical filaments were compared to the available mean orbits in the IAU-MDC list of all showers. The orbital similarity was evaluated with the help of the Southworth-Hawkins (1963) D_SH discriminant, where the identification was positive if D_SH ≤ 0.10. This D_SH discriminant-based identification is given in Table A.5. The list of showers in Table A.5 appears to be a little different than that in Table A.4. The set of IAU-MDC showers identified with the showers separated from the databases considered, in Table A.4, also contains showers #156, #172, #173, #215, and #485, which were not identified on the basis of mean orbits, and are therefore missing from Table A.5. On contrary, showers #256, #257, #631, and #633 were identified on the basis of the latter and are in Table A.5, but are missing from Table A.4. Of theses showers only the Southern χ-Orionids (#257) are the established shower. The difference between the lists of showers in Tables A.4 and A.5 reflects the still persisting problem that the characteristics of especially minor showers are poorly determined, dependent on the particular database and method of shower separation.

Fig. 6 Motion of the mean radiant of the parts of northern showers separated from the IAU-MDC CAMS database during the period of the activity. The equatorial coordinates of the mean radiant are approximated with the linear functions of time.

Fig. 7 Panels a and c: overall positions on the sky of radiants of theoretical filaments successfully identified with at least one real shower separated from any of considered databases. Panels b and d: radiants of meteors of individual real showers. These radiants are distinguished by various color symbols for different showers. Panels a and b: radiants are shown in the Hammer-Aitoff projection of the celestial sphere in the common equatorial coordinates. The sinusoid-like curve is the ecliptic. Panels c and d: same radiants are shown in the apex-centered ecliptical coordinates. The radiant areas are in the northern anti-helion (NAH), southern anti-helion (SAH), northern helion (NH), and southern helion (SH) cardinal directions.

Fig. 8 Positions of the radiants of predicted (various black symbols) and real meteors separated from the C (blue asterisks), S (green crosses), E (orange empty squares), R (violet triangles), and F (red plus-symbols) databases. Panels a and b: radiants of showers #726 and #628, respectively, and the predicted filaments corresponding to these showers are shown. The apex-centered ecliptical coordinate frame is used.

6 Conclusions

We created a grid of the models of meteoroid stream released from the short-period comet 2P/Encke. This modeling resulted in a proof that the comet is the parent body of 16 showers from the IAU-MDC list of all showers and of another 5 showers that are not in the list at the moment.

The radiants of the showers related to 2P are situated in all cardinal directions of ecliptic showers, helion and anti-helion, whereby both have the northern and southern strands. Specifically, the radiant area of shower #172 (#156 and #173) is situated in the northern (southern) helion direction. The radiants of showers #17, #215, #629, #632, #634, #635, #726, and a shower not yet listed in the IAU-MDC list of all showers that we named η-Piscids (#2, #485, #624, #626, #628, #636, and four previously unknown showers that we named as γ-Geminids, ω2-Aquariids, 74 Orionids, and 29 Piscids) are situated in the northern (southern) anti-helion direction. Only five ofthese showers are classified as established: Norther Taurids (#17), Daytime ζ-Perseids (#172; northern showers), Southern Taurids (#2), Southern Daytime May Arietids (#156), and Daytime β-Taurids (#173; southern showers).

The radiant areas in the northern helion (southern helion) cardinal direction are symmetrical to those in the southern anti-helion (northern anti-helion) direction with respect to the apex of the Earth’s motion around the Sun.

It is possible to conclude that comet 2P is most probably the significant source of the showers belonging to the Taurid complex. In the process of the separation of real showers, on the basis of predicted mean orbits, we also separated several showers that were too different from their predicted counterparts, and thus we did not include them in the 2P stream. These showers indicate the presence of other parent bodies in the adjacent orbital phase space.

Acknowledgements

This article was supported by the realization of the Project ITMS No. 26220120029, based on the supporting operational Research and development program financed by the European Regional Development Fund. The work was also supported in part by the Slovak Grant Agency for Science, VEGA, grant No. 2/0037/18, and by the Slovak Research and Development Agency under contract No. APVV-16-0148.

Appendix A Additional tables

References

All Tables

All Figures

Fig. 1 Behavior of perihelion distance (panel a), semi-major axis (panel b), eccentricity (panel c), and inclination to the ecliptic (panel d) of the initially nominal orbit of comet 2P/Encke (red curve) and the orbits of the first set of its clones (green curves) mapping the uncertainty of the determination of the nominal orbit. The evolution is reconstructed backward for 50 000 yr. A non-gravitational force is ignored.

In the text

Fig. 2 Behavior of perihelion distance (panel a), semi-major axis (panel b), eccentricity (panel c), and inclination to the ecliptic (panel d) of the initially nominal orbit of comet 2P/Encke, which is (violet curve) or is not (red curve) influenced by the non-gravitational force (jet effect). The evolution of the corresponding elements of the orbits of the second set of its clones (green curves), which delimit the uncertainty of non-gravitational parameters A1 and A2, are also shown. The evolution is reconstructed backward for 50 000 yr.

In the text

Fig. 3 Evolution of the minimum distance between the orbital arcs of comet 2P/Encke and Earth from 50 000 yr before the present up to the present. The minimum distance of the post-perihelion (pre-perihelion) arc is shown by the violet (cyan) curve when the evolution of the 2P orbit was followed ignoring the non-gravitational force. When this force is taken into account, the minimum distance is shown by the green (brown) curve.

In the text

	Fig. 4 Positions of the orbital nodes of comet 2P/Encke during the last 50 000 yr. The green circle indicates the orbit of the Earth. The red (blue) curve shows the positions of ascending (descending) node.
In the text

Fig. 5 Panel a: number density of the theoretical radiants within the space angles of size 1° × 1°. With respect to the maximum number density, N, the number densities are divided into eight levels, corresponding to intervals delimited by 0, N∕8, 2N∕8, …, N; the seven most numerous intervals are color-coded according to the legend in the plot. The model for tev = 8 kyr and β = 0.007 is used as an example. Panel b: radiants of all Earth-approaching particles in the model.

In the text

	Fig. 6 Motion of the mean radiant of the parts of northern showers separated from the IAU-MDC CAMS database during the period of the activity. The equatorial coordinates of the mean radiant are approximated with the linear functions of time.
In the text

Fig. 7 Panels a and c: overall positions on the sky of radiants of theoretical filaments successfully identified with at least one real shower separated from any of considered databases. Panels b and d: radiants of meteors of individual real showers. These radiants are distinguished by various color symbols for different showers. Panels a and b: radiants are shown in the Hammer-Aitoff projection of the celestial sphere in the common equatorial coordinates. The sinusoid-like curve is the ecliptic. Panels c and d: same radiants are shown in the apex-centered ecliptical coordinates. The radiant areas are in the northern anti-helion (NAH), southern anti-helion (SAH), northern helion (NH), and southern helion (SH) cardinal directions.

In the text

Fig. 8 Positions of the radiants of predicted (various black symbols) and real meteors separated from the C (blue asterisks), S (green crosses), E (orange empty squares), R (violet triangles), and F (red plus-symbols) databases. Panels a and b: radiants of showers #726 and #628, respectively, and the predicted filaments corresponding to these showers are shown. The apex-centered ecliptical coordinate frame is used.

In the text