© The Authors 2023

1 Introduction

In the dynamics of the N-body problem, with N > 2, the geometry of the orbits of the bodies changes over time due to their mutual gravitational interactions. In particular, the angles that specify the orientation of the osculating ellipse precess at a secular rate. In the context of our Solar System, the precession frequencies of the longitude of the perihelion, ϖ, and of the longitude of the ascending node, Ω, of the orbits of the planets are typically denoted with ɡ_j and S_j, respectively, where the index j = 1,…,8 provides labels for the planets from Mercury to Neptune. Accurate values of planetary secular frequencies were previously computed, for instance, by Nobili et al. (1989) and Laskar et al. (2004, 2011).

The presence of the planets also cause both short and long-term perturbations on the motion of small Solar System bodies, such as asteroids. The amplitude of the secular perturbations increases near a secular resonance, that occurs when the precession frequency, g, of the longitude of the perihelion, or the frequency S of the longitude of the node, of the asteroid is nearly equal to that of a planet. Linear secular resonances involve only one asteroid's frequency and one planet's frequency, and resonances of type ɡ ≃ ɡ_j are typically denoted with V_j, while those of type s ≃ S_j are denoted with $v_{1_j}$.

The problem of locating the secular resonances in the main belt has gone through a series of improvements and it now stands as a well-known topic in asteroid dynamics. The first attempt to produce maps of secular resonances in the main belt was done by Williams (1969), who used a semi-analytical method to derive the asteroids' proper elements and proper frequencies. The method was later refined by Williams & Faulkner (1981) to better locate the three strongest secular resonances, namely ν₅, ν₆, and ν₁₆. Later on, a simple analytical model was used by Yoshikawa (1987) to locate the ν₆ secular resonance more precisely. More accurate analytical theories, where the perturbing function is expanded in Fourier series of cosines of combinations of the angular variables and the amplitudes are expressed as polynomials in eccentricity and sine of inclination, were developed by Milani & Knežević (1990, Milani & Knežević 1994). The authors were able to produce maps of both linear and non-linear secular resonances in the main belt, for small proper eccentricity and proper inclination. The same method was used by Knežević et al. (1991) to map the location of all linear secular resonances with the outer planets, for heliocentric distances between 2 and 50 au. A semi-analytical theory that avoids the expansion in eccentricity and inclination has been developed by Morbidelli & Henrard (1991a,b), and it permitted the authors to locate the linear secular resonances at higher values of eccentricity and inclination. More recently, the introduction of a synthetic theory for the computation of proper elements and frequencies (see Knežević & Milani 2000, 2003) led to attempts to locate secular resonances based on purely numerical techniques. A first empirical method was used by Tsirvoulis & Novaković (2016), who identified a large population of synthetic asteroids in a secular resonance with (1) Ceres and (4) Vesta by matching their secular frequencies up to a certain threshold. A new, completely synthetic method was recently developed by Knežević & Milani (2019). This method was used by Knežević (2020) to significantly improve the location of the ν₅ resonance, and by Knežević (2022) to review the location of the linear and several non-linear secular resonances in the main belt.

On the contrary, little attention has been paid to the location of secular resonances in the region of near-Earth objects (NEOs), which are asteroids with a perihelion distance q smaller than 1.3 au. So far, the first and only attempt was undertaken by Michel & Froeschlé (1997), who used a semi-analytic model similar to that of Kozai (1962) to locate the linear secular resonances at semi-major axis a between 0.5 and 2 au, only for a value of proper eccentricity equal to 0.1. However, NEOs can reach much higher eccentricities, and they may cross the orbit of one or more planets. Orbit crossings were indeed the reason for which Michel & Froeschlé (1997) were not able to compute maps at higher proper eccentricities. In turn, locating the secular resonances inside the NEO region may help improve our understanding of the global dynamics of NEOs. From purely numerical simulations, it is now known that secular resonances play a fundamental role in the delivery of objects from the main belt to the NEO region (see, e.g., Scholl & Froeschle 1986; Froeschle & Scholl 1986, 1989; Granvik et al. 2017). In particular, the ν₆ secular resonance is the most important mechanism for the production of NEOs and it provides the largest fraction of Earth's impactor (Bottke et al. 2002; Granvik et al. 2018). Although planetary close encounters are one of the main perturbations in the motion of NEOs, secular resonances may still play a significant role in their dynamics (Foschini et al. 2000; Gladman et al. 2000).

In this work, we located the linear secular resonances with the planets from Venus to Saturn inside the NEO region, for a semi-major axis between 0.5 and 3 au. For this purpose, we used the semi-analytic method by Gronchi & Milani (2001) to propagate the secular dynamics of NEOs, from which the proper frequencies ɡ and s are extracted. Proper frequencies are then computed on a grid of points in the phase spaces of (a, e₀) or (a, i₀), where e₀ are i₀ are the proper eccentricity and proper inclination. Maps of the secular resonances with planets from Venus to Saturn were then determined by computing the level curves corresponding to planets' proper frequencies. The method by Gronchi & Milani (2001) permits us to compute the secular evolution even beyond an orbit crossing, thus enabling us to produce secular resonance maps even at high eccentricity.

The paper is structured as follows. In Sect. 2, we introduce the semi-analytical Hamiltonian model for the computation of the secular evolution of NEOs and we describe the method to locate the secular resonances. In Sect. 3, we show the location of the secular resonances obtained with our method and show some comparisons with full numerical simulations. In Sect. 4, we report on the implications of our results with respect to the global dynamics of NEOs and we discussed the effects of mean-motion resonances. Finally, in Sect. 5, we provide our conclusions.

2 Methods

2.1 Hamiltonian secular model

We used the averaged model by Gronchi & Milani (2001) to compute the secular evolution of an asteroid. The planets from Venus to Saturn are included in the model and they are placed on circular co-planar orbits. We denoted the Gauss constant as $\text{k}=\sqrt{𝒢m_0}$, where 𝒢 is the universal gravitational constant, and m₀ is the mass of the Sun. We also set μ_j = m_j/m₀, j = 2,…, 6, where m_j is the mass of the jth planet. The canonical Delaunay variables are defined as $$\begin{array}{rr}\hfill L=\text{k}\sqrt{a},& \hfill \ell =M,\\ \hfill G=\text{k}\sqrt{a\left(1-e^2\right)},& \hfill g=\omega ,\\ \hfill Z=\text{k}\sqrt{a\left(1-e^2\right)}\cos i,& \hfill z=\text{Ω,}\end{array}$$(1)

where a is the semi-major axis, e is the eccentricity, i is the inclination, M is the mean anomaly, ω is the argument of the pericenter¹, and Ω is the longitude of the ascending node. The Hamiltonian describing the motion of the asteroid is: $$\begin{array}{cc}\mathscr{H}={\mathscr{H}}_0+\epsilon {\mathscr{H}}_1,& \epsilon ={\mu }_5\end{array},$$(2)

where ℋ₀ is the unperturbed Keplerian part, and ℋ₁ accounts for the perturbations of the planets: $$\begin{array}{l}{\mathscr{H}}_0=-\frac{{\text{k}}^4}{2L^2},\hfill \\ {\mathscr{H}}_1=-{\text{k}}^2{\displaystyle \sum _{j=2}^6\frac{{\mu }_j}{{\mu }_5}\left(\frac{1}{\left|r-r_j\right|}-\frac{r\cdot r_j}{{\left|r_j\right|}^3}\right).}\hfill \end{array}$$(3)

In Eq. (3), r and r_j denote the heliocentric position of the asteroid and the jth planet, respectively.

Assuming that there are no mean-motion resonances between the asteroid and any of the planets, the mean anomalies ℓ and ℓ_j of the asteroid and the j-th planet evolve much faster than the other orbital elements. The secular evolution at first order in the planetary masses is then obtained by the double average of the Hamiltonian of Eq. (2) over the fast angles, namely, $$\overline{\mathscr{H}}=-\epsilon {\text{k}}^2{\displaystyle \sum _{j=2}^6\frac{{\mu }_j}{{\mu }_5}{\displaystyle {\int }_0^{2\pi }{\displaystyle {\int }_0^{2\pi }\frac{1}{\left|r-r_j\right|}\text{d}\ell \text{d}{\ell }_j.}}}$$(4)

With the assumptions of zero eccentricity and inclination for the orbits of the planets, this model can be seen as the dominant term of an expansion in powers of the planetary eccentricities and inclinations of the double-averaged Hamiltonian in the case of eccentric and inclined orbits of the planets (Thomas & Morbidelli 1996). We note that the indirect perturbation is not present in the averaged Hamiltonian, $\overline{\mathscr{H}}$, since its double integral vanishes. In this model, L is a first integral, hence, the semimajor axis, a, remains constant. Moreover, Ω does not appear in the Hamiltonian of Eq. (4) because the planets are placed on circular and co-planar orbits, and the conjugate momentum, Z, is also a first integral. The system defined by $\overline{\mathscr{H}}$ has therefore only one degree of freedom and it is integrable. Delaunay elements are classically used to study this problem (see e.g., Kozai 1962; Thomas & Morbidelli 1996; Gallardo et al. 2012; Saillenfest et al. 2016), and this is because the Hamiltonian in these coordinates can be naturally reduced to a Hamiltonian of one degree of freedom, depending only on (G, ɡ).

The averaged Hamiltonian $\overline{\mathscr{H}}$ has a first order polar singularity when the orbit of the asteroid crosses that of a planet (Gronchi & Milani 1998). By denoting with y ∈ {G, Z, ɡ, z} one of the Delaunay elements, the derivatives $\partial \overline{\mathscr{H}}/\partial y$ can be computed numerically simply by exchanging the derivative and the integral signs when there are no orbit crossings. On the other hand, at orbit crossings the double integral of the derivatives of ℋ₁ is divergent, and therefore the theorem of differentiation under the integral sign does not apply. Gronchi & Milani (1998), Gronchi & Tardioli (2013) showed that, in a neighborhood of an orbit-crossing configuration, the vector field defined by the averaged Hamiltonian has a two-fold definition that can be computed analytically. By using this property, Gronchi & Milani (2001) set up a numerically stable algorithm that allows us to compute a solution for Hamilton's equations for $\overline{\mathscr{H}}$ beyond an orbit-crossing configuration. In this work, we use this algorithm to propagate the secular evolution of an asteroid beyond an orbit-crossing singularity. Details about the method and the practical numerical implementation can be found in Gronchi & Milani (1998, 2001), Gronchi & Tardioli (2013), Fenucci et al. (2022).

2.2 Locating the secular resonances

Proper frequencies can be computed as in Gronchi & Milani (2001). A solution of the Hamiltonian secular system defined by $\overline{\mathscr{H}}$ is computed numerically and, if ω circulates, we can find the times, t₀ and t_π/2, such that ω(t₀) = 0, ω(t_π/2) = π/2. We also set Ω₀ = Ω(t₀), Ω_π/2 = Ω(t_π/2). The Hamiltonian, H, is even and π-periodic in ω, therefore the frequency, ɡ − s, of the argument of perihelion, ω, is given by: $$g-s=\frac{2\pi }{4\left(t_{\pi /2}-t_0\right)},$$(5)

while the frequency, s, of the longitude of the node, Ω, is given by: $$s=\frac{{\text{Ω}}_{\pi /2}-{\text{Ω}}_0}{t_{\pi /2}-t_0}.$$(6)

The dynamical state in which ω librates is also called Kozai resonance. In this work, we assume that ω librates around 0°, therefore, we take t₀ such that ω(t₀) = 0 and $\dot{\omega }\left(t_0\right)>0$. Then, we compute the time, Δt, after which ω vanishes again and we denote by ΔΩ the change in Ω during this timespan. The frequency, ɡ − s, of ω vanishes, but we can still define a libration frequency, lf, and the frequency, s, for Ω, as follows: $$\begin{array}{cc}\text{lf}=\pm \frac{2\pi }{2\text{Δ}t},& s=\frac{\text{ΔΩ}}{\text{Δ}t}\end{array}.$$(7)

The formula for the libration frequency, lf, is justified by the symmetry properties of $\overline{\mathscr{H}}$. Besides the proper frequencies, the proper eccentricity and the proper inclination can be defined as e₀ = e(t₀) and i₀ = i(t₀), respectively. This is a slightly different definition from that adopted by Gronchi & Milani (2001), who used the minimum eccentricity, e_min, and the corresponding maximum inclination, i_max, attained alongside a solution. However, the two definitions coincide in the case that the solution does not cross any planet's orbit (Michel & Froeschlé 1997). We note that ω may also librate around other values and, therefore, it may never vanish during a libration period (see, e.g., Michel & Froeschlé 1997; Gronchi & Milani 1998, 2001). While these regions certainly exist, they do not appear in our results, because they would need another special definition of the proper elements based on the center of their libration island.

The method to identify the secular resonances is similar to that of Michel & Froeschlé (1997). The initial longitude of the node Ω is set to 0, but any other value would not change the result since Ω does not appear in the Hamiltonian $\overline{\mathscr{H}}$. The initial argument of pericenter ω is also set to 0, so that the initial eccentricity and inclination also correspond to the proper values. Then, we go on to fix the initial value of the proper eccentricity e₀ (or the initial proper inclination i₀), and choose a grid in the (a, i₀)-plane (or in the (a, e₀)-plane). The proper frequencies g and s are computed for each point of the grid. The location of the secular resonances v_j, $v_{1_j}$, j = 2,…,6, in the (a, i₀)-plane or in the (a, e₀)-plane, is then given by the contour lines corresponding to the levels ɡ = ɡ_j and s = s_j, j = 2,…, 6, respectively. In this work, we used the values of the secular frequencies of the planets determined by Laskar et al. (2011), reported in our Table 1.

3 Results

3.1 Secular resonance maps at fixed inclination

We discretized the (a, e₀)-plane with a step of 0.01 au for the semi-major axis and 0.008 for the proper eccentricity. We took into account the region defined by 0.5 au < a < 3 au and 0.01 < e₀ < 0.8. The proper inclination was fixed and we computed the maps for i₀ = 2°, 5°, 10°, 20°, 30°, and 40°.

The left column of Fig. 1 shows the location of the v_j, j = 2,…,6 secular resonances, while the right column shows that of $v_{1_j}$, j = 2, 3, 4, 6, for values of proper inclination i0 equal to 40°, 30°, 20°, and 10°. Figure 2 shows the location of the v_j, j = 2,…,6 resonances at lower inclination, namely for i₀ = 5°, 2°. Figures 1 and 2 also show the curves q = a₂, a₃, a₄ and Q = a₂, a₃, a₄, where q = a(1 − e) and Q = a(1 + e) are the perihelion and the aphelion distance of the asteroid, respectively, and a_j denote the semi-major axis of the planets.

From the left column of Fig. 1 we can see that v_j, j = 2, 3, 4, 5 (i.e., resonances with the inner planets, and with Jupiter) always appear inside the NEO region, while the v₆ secular resonance with Saturn appears to a large extent at proper inclinations i0 smaller than 30^°. The resonances, V3 with the Earth and v₄ with Mars, are always close to each other because the precession rates of the longitude of the perihelion of these two planets are similar (see Table 1). At i₀ = 40^°, all the secular resonances are confined at values of eccentricity e₀ larger than 0.2, and at semi-major axis a larger than about 1.3 au. As the proper inclination i0 decreases, the secular resonances move towards smaller values of semi-major axis and they span an increasing interval of eccentricity. It is worth noting that when passing from i₀ = 30^° to i₀ = 20°, there is a significant change in the location of the v₂ and v₅ resonances. For inclinations larger than 30^°, these two resonances are present in the region q < a₄ and they are confined at rather low eccentricity, while for q > a₄, they extend towards larger eccentricity values, up to the limit of 0.8 we considered. For inclinations smaller than 30°, they disappear from the region q < a₄ and they remain confined at semi-major axis smaller than 1.5 au, where they assume values of eccentricity from near 0 up to 0.8. The remaining resonances do not show this change in inclination between 10^° and 40^°. Additionally, the V5 resonance shows a second branch at a < 1 au for i₀ = 20°; the same happens to v₃ and v₄ for i₀ = 10^°.

The secular resonance maps for i₀ = 2^° and 5^° (see Fig. 2) become more complex. At i₀ = 5^°, the V3 and V4 disappear from the region of q < a₄, while a new almost-circular branch appears in the region of a₃ < q < a₄. In the region of q > a₃, these two resonances are close to the curve q = a₃ for semi-major axis smaller than about 1.7 au, at which point they start extending towards larger values of eccentricity. A similar feature is seen for v₂ and v₅, which tend to follow the curves q = a₃ and q = a₂ for certain intervals of semi-major axis values. On the contrary, the v₆ resonance does not change much with respect to the one obtained at i₀ = 10^°. We also note that more branches of all the resonances appear at semi-major axis smaller than about 1.2 au and most of them follow the curves of constant, Q. The global picture for i₀ = 2^° is similar to the previous case, with the difference that the alignment along the curves q = a₂ and q = a₃ is more definite and it also appears for the ν₆ secular resonance. These results show that when the orbital elements get closer to the curves of perihelion distance (or, to a minor extent, aphelion distance) equal to that of a planet, the value of the proper frequency ɡ is significantly affected. Moreover, this effect increases as the proper inclination decreases.

The right column of Fig. 1 shows the location of the resonances $v_{1_j}$, j = 2, 3, 4, 6, involving the longitude of the node. The v₁₂ resonance basically never appears in the region 0.5 au < a < 3 au for inclination i₀ up to 40°. This was already noticed by Michel & Froeschlé (1997), who found that ν₁₂ appears only at semi-major axis smaller than 0.5 au, at least for proper eccentricity of 0.1. Our results confirm this feature even at larger eccentricity values, hence, this resonance is not expected to play any major role in the dynamics of NEOs. The resonances ν₁₃ and ν₁₄ involving the Earth and Mars are again close to each other, since the node of their orbits also precesses at a similar rate (see Table 1).

At an inclination, i₀ = 40^°, we can see that v₁₃ and v₁₄ appear at a < 2 au, while ν₁₆ is present in the region a < 2.5 au. These three resonances all span a large interval of eccentricities, up to e₀ ≈ 0.7. Moreover, ν₁₃ and ν₁₆ extend towards the limit of e₀ = 0.8 that we used in our computations. As the proper inclination decreases, the secular resonance curves shrink, and they move towards smaller values of semi-major axis and proper eccentricity. At i₀ = 30^°, the qualitative picture is the same as that obtained at i₀ = 40°. At i₀ = 20°, the ν₁₆ resonance does not change its qualitative shape. On the other hand, v₁₃ and ν₁₄ are composed by four different separated curves, confined at eccentricity smaller than about 0.5. Similarly to the case of the v_j resonances, the curves of constant q, Q equal to the semimajor axis of a planet play a role in determining the proper frequencies s, and the location of the secular resonances. At i₀ = 10°, v₁₆ has a different qualitative structure. It shows a closed curve in the region a₂ < q < a₃, and another branch appears at q < a₃. This is the only resonance appearing to a large extent, while all the other ones either disappeared or are confined to a very small region of the phase space. At lower values of proper inclination, i₀, we found that all the secular resonances involving the nodal longitude disappear.

The Kozai resonance does not involve the motion of the planets and it occurs when ɡ − s = 0 or, equivalently, when the argument of the perihelion ω librates. In Figs. 1 and 2, the Kozai resonance region is identified in gray. At i₀ = 40° Kozai librators are rare and, indeed, the gray area is small. At i₀ = 30^° and 20^°, the Kozai resonance appears to a large extent in three different regions: near a = a₂ and a = a₃ at small eccentricity and at the intersection between q = a₂ and Q = a₃. At inclinations i₀ ≤ 10°, it is only the Kozai librators near the semi-major axes of Venus and the Earth that are left, while a small gray area appears near a = a₄. We also note that the extent of the gray area decreases as the proper inclination value decreases.

Fig. 1 Map of secular resonances in the (a, e0)-plane for a between 0.5 and 3 au, for different values of proper inclination, i0. The left column shows the locations of vj, j = 2,3,4,5,6, while the right column shows the locations of $v_{1_j}$, j = 2,3,4,6. The color code of each resonance is indicated by the legend in the top row panels. From the top to the bottom, rows correspond to values of i0 equal to 40°, 30°, 20°, and 10°, respectively. The dashed curves correspond to q = a2, a3, a4 and Q = a2, a3, a4 where q = a(1 − e), Q = a(1 + e) are the perihelion and the aphelion distance of the asteroid, respectively. We also reported the location of the Kozai librators in grey, identified by g − s = 0.

Fig. 2 Map of the vj, j = 2,3,4,5,6 secular resonances in the (a, e0)-plane for a between 0.5 and 3 au, at i0 = 5° (left panel), and at i0 = 2° (right panel). Dashed curves correspond to q = a2, a3, a4 and Q = a2, a3, a4 We also reported the location of the Kozai librators in grey, identified by ɡ − s = 0.

3.2 Secular resonance maps at fixed eccentricity

We discretized the (a, i₀)-plane with a step of 0.01 au for the semi-major axis, and 0.5° for the inclination. We took into account the region defined by 0.5 au < a < 3 au and 0° < i₀ < 50°. The proper eccentricity was fixed, and we computed the maps for values e₀ = 0.1 · h, h = 1,…,7, and for low eccentricity e₀ = 0.02, 0.05. The left column of Fig. 3 shows the location of the v_j,_j = 2,…,6 secular resonances, while the right column shows the location of $v_{1_j}$, j = 2, 3, 4, 6, for values of proper eccentricity, e₀, equal to 0.1, 0.2, 0.4, and 0.6. We note that the first row of Fig. 3 is the same case considered by Michel & Froeschlé (1997). The results we obtained are in really good agreement with those presented by Michel & Froeschlé (1997) (see Fig. 2 and 3 therein), with the difference that we are able to extend the maps to the planet-crossing regions. Figure 4 shows the same maps obtained for e₀ = 0.05, 0.02.

From the left column of Fig. 3, we can see that all the resonances v_j, j = 2,…, 6 appear in the region of 0.5 au < a < 3 au. As expected, V3 and V4 are always close to each other. At a proper eccentricity of e₀ = 0.1, the v₃, v₄, and v₆ resonances all have four different branches, each one appearing in the regions a > a₄, a₃ < a < a₄, a₂ < a < a₃, and in a < a₂. The v₆ resonance with Saturn is confined at inclination smaller than about 20°. The resonances v₂ and v₅ have three different branches, each one appearing for a > a₃, a₂ < a < a₃, and in a < a₂. We can also notice the presence of Kozai librators around a = a₂ and a = a₃, which fill this area of the phase space up to i ≈ 35°.

At a proper eccentricity e₀ = 0.2, the branches in the region a₂ < a < a₃ have almost disappeared, while those that were appearing at a < a₂ are shifted to smaller values of semi-major axis, if compared with the previous case. For a > a₃, the map is qualitatively similar to the case e₀ = 0.1. The Kozai librators still appear near a₂ and a₃, but they fill a smaller area of the phase space, which is also moved towards larger inclinations. Moreover, another region filled by Kozai librators appears between Venus and the Earth, that corresponds to that seen in Fig. 1 for a proper inclination of 30^°, and 40^°. We found that the location of the secular resonances for e₀ = 0.3 is qualitatively similar to e₀ = 0.2 overall and, therefore, it is not shown.

The maps for e₀ = 0.4 and e₀ = 0.6 are also qualitatively similar to each other. The region at a ≲ 1.3 au is almost cleared out from all the secular resonances, which appear to a large extent at larger semi-major axis values, where they only have one branch. Moreover, they cover a larger interval of inclination values than in the previous cases. It is also worth noting that Kozai librators have disappeared, which is consistent with the results presented in Figs. 1 and 2.

The left column of Fig. 4 shows the maps obtained for e₀ = 0.05 and 0.02. We can see they are qualitatively similar to what we obtained for e₀ = 0.1. We can see that the Kozai resonance region is smaller than for the case e₀ = 0.1, and the secular resonances extend closer to the semi-major axes of Venus, of the Earth, and of Mars.

The right column of Fig. 3 shows the location of the $v_{1_j}$, j= 2, 3, 4, 6 secular resonances. As seen in Sect. 3.1, the resonance ν₁₂ with Venus is negligible. At an eccentricity of e₀ = 0.1, we can see that ν16 has three different branches, while v₁₃ and ν₁₄ have only one. These three resonances are confined at semimajor axis smaller than 2.3 au and they span a large interval of inclination values. We note that ν₁₃, ν₁₄, and ν₁₆ appear also in the Kozai resonance region. For a higher value of proper eccentricity, we can see that ν₁₆ also have a single branch. Moreover, as e0 increases, their location moves towards smaller values of the semi-major axis and they span a smaller interval of inclinations. Finally, the right column of Fig. 4 shows the location of the $v_{1_j}$ resonances for e₀ = 0.02 and 0.05. The maps are qualitatively similar to the case e₀ = 0.1, with the difference that the three resonances involving the Earth, Mars, and Saturn all extend towards larger semi-major axis values.

Fig. 3 Map of secular resonances in the (a, i0)-plane for a between 0.5 and 3 au, for different values of proper eccentricity, e0. The left column shows the locations of vj, j = 2, 3, 4, 5, 6, while the right column shows the locations of $v_{1_j}$, j = 2, 3, 4, 6. The color code of each resonance is indicated by the legend in the top row panels. From top to the bottom: rows correspond to values of e0 equal to 0.1, 0.2, 0.4, and 0.6, respectively. We also report the location of the Kozai librators in grey, identified by ɡ − s = 0.

3.3 Frequencies near perihelion equal to planets' semi-major axes

Figures 1 and 2 show that the secular resonances tend to follow the perihelion curves a(1 − e) = a_j, j = 2, 3, 4, and the trend is more prominent when the proper inclination, i0, decreases toward 0. To better understand the reason behind this effect, we computed the proper frequencies on a grid in eccentricity by fixing both the values of the semi-major axis and of the proper inclination. To gain a better picture of the behavior of the frequency, we also used a step in eccentricity of 5 × 10⁻⁴, which is smaller than that used to produce Figs. 1 and 2. Figure 5 shows the frequency ɡ computed for a = 1.2 au and i₀ = 2°, as a function of the proper eccentricity, e₀. When e₀ approaches the critical values $e_0^{\left(j\right)}$, j = 2, 3 for which $a\left(1-e_0^{\left(j\right)}\right)=a_{j,}j=2,3$, the frequency ɡ sharply increases. Moreover, the magnitude of the increase is larger and larger as the inclination tends to 0^°. This feature is caused by the model itself, and it is due to the orbit-crossing singularity.

Figure 5 shows also an interesting feature: between $e_0^{\left(3\right)}\approx 0.1667$ and $e_0^{\left(2\right)}\approx 0.3972$, the frequency, ɡ, has a minimum that is smaller than ɡ₅, which is the smallest planetary frequency. As e₀ tends to $e_0^{\left(3\right)}$ from the right, the frequency, ɡ, increases sharply, and it crosses all the planetary frequencies ɡ_j, j = 2,…, 6 within a small interval of proper eccentricities, e₀. The same happens as e₀ tends to $e_0^{\left(2\right)}$ from the right, with the difference that only ɡ₂ and ɡ₅ are attained near the critical value. This effect causes the secular resonances to follow the Earth perihelion curve and to be close to each other. In addition, the fact that ɡ grows so fast in such a small interval of eccentricity will make the resonance width very small, so that objects following the perihelion curves by effect of the secular resonances will be extremely rare to find, because they would be pushed outside of the resonance by other small perturbations. In fact, we did not find any object following a planet's perihelion curve in our purely numerical integrations (see Sect. 3.4).

The sharp increase of the frequency ɡ near planets' perihelion curves also causes some small numerical errors in the global plot of the maps presented in Figs. 1 and 2. Because we use a discrete grid in the (a, e₀)-plane, it may happen that a value larger than one of the planets' frequency, ɡ_j, is not achieved in the neighborhood of one of such perihelion curves. As a consequence, the corresponding secular frequency v_j would not appear in the plot. This can be seen, for instance, in the maps at i₀ = 2°, presented in Fig. 2. Figure 5 shows that all the five resonances v_j, j = 2,…, 6 are crossed in a neighborhood of a = 1.2 au and e₀ ≈ 0.17, while it is not clear that this would happen in the left panel shown in Fig. 2. This is caused by the fact that the grid is not dense enough in this tiny neighborhood. So, to prove this assumption, we recomputed the secular frequencies in the portion (a, e₀) ∈ [1.15 au, 1.3 au] × [0.1, 0.3], using a step of 0.004 au in semi-major axis and of 5 × 10⁻⁴ in proper eccentricity. Figure 6 shows the location of the secular resonances in this portion of the phase space. Here, we can see that all the five secular resonances, v_j, appear close to q = a₃, even though v₃,,v₄, and v₆ are barely distinguishable since they are very close to each other. Overall, we found these numerical effects to happen at low proper inclination, i₀, and only near q = a_j, j = 2, 3, 4, where the width of the resonance is very small in any case. Therefore, using a denser grid would not change the global picture of the secular resonances shown in Sect. 3.1, which was obtained using a less dense grid.

Fig. 4 Map of secular resonances in the (a, i0)-plane for a between 0.5 and 3 au, for different values of proper eccentricity, e0. The left column shows the locations of vj, j = 2, 3, 4, 5, 6, while the right column shows the locations of $v_{1_j}$, j = 2, 3, 4, 6. The color code of each resonance is indicated by the legend in the top-left panel. From top to the bottom: Rows correspond to values of e0 equal to 0.05, 0.02, respectively. We also report the location of the Kozai librators in grey, identified by ɡ − s = 0.

Fig. 5 Frequency, ɡ, (black solid curve) computed on a grid in proper eccentricity, e0, for a fixed value of semi-major axis of a = 1.2 au and proper inclination of i0 = 2°. The horizontal straight lines correspond to the values ɡj, j = 2, 3, 4, 5 of the proper frequencies of the planets. The vertical dashed lines correspond to the values of eccentricity, such that a(1 − e0) = aj, j = 2, 3.

Fig. 6 Zoom on the map of the secular resonances vj, j = 2, 3, 4, 5, 6 at i0 = 2°, near a = 1.2au and e0 = 0.17.

3.4 Comparison with numerical integrations

We performed purely numerical integrations to confirm the secular resonance locations predicted by the theory in Sects. 3.1 and 3.2. To account for the fact that the maps were computed in the proper elements space, we took initial orbital elements along the curves ɡ = ɡ_j (or S = S_j) and ɡ = ɡ_j ± 1" yr⁻¹ (or S = S_j ± 1" yr⁻¹), for a selected planet j ∈ {2, 3, 4, 5, 6}. The initial argument of pericenter, ω, and the initial argument of the node, Ω, were both set to 0^°, while the initial mean anomaly, ℓ, was chosen randomly between 0^° and 360^°. Numerical integrations were performed by using a hybrid symplectic scheme (Chambers 1999) that is able to handle close encounters with planets, that is included in the MERCURY² package (see also Fenucci & Novakovic 2022). The gravitational attraction of the Sun and of all the planets from Mercury to Neptune were included in the model and the initial conditions for planets at epoch 58800 MJD were taken from the JPL Horizons³ ephemeris system. We note that here we take into account a complete model of the Solar System to show that the simplified model presented in Sect. 2 is actually able to predict the location of secular resonances in a qualitative picture. We propagated the evolution of the asteroids and the planets for a total time of 2 Myr by using a timestep of one day.

We present here the results obtained by choosing an initial inclination of 20^°. Figure 7 shows the evolution of a test asteroid selected with the method described above, which is affected by the ν₆ secular resonance. The critical angle, ϖ − ϖ_S librates around 0^°, while the eccentricity decreases down to values of about 0.4 during the first 0.2 Myr of evolution, and then suddenly increases up to values close to 1, causing this object to impact with the Sun in less than 0.5 Myr of dynamical evolution. This is a well-known effect of the ν6 resonance, found for the first time by Farinella et al. (1994), which is able to deliver Sun impactors even starting from the main belt at a very low eccentricity. This simulation shows that, as expected, this behavior still occurs well inside the NEO region, even at higher values of the initial eccentricity.

Figure 8 shows a test asteroid affected by the ν5 secular resonance. The critical angle ϖ − ϖ_J librates around 0° during the 2 Myr of evolution and the eccentricity passes from an initial value of about 0.2 to values close to 0.6. The semi-major axis assumes values between 1.2 and 1.5 au, and it undergoes a random walk evolution, caused by close encounters with the inner planets. Despite the changes in the semi-major axis, the critical angle still librates and it is not removed from the secular resonance. Although this object did not end the evolution on a collision with the Sun, ν₅ was able to push it to large eccentricity values, showing a behavior similar to that of ν₆. By using numerical simulations, Gladman et al. (2000) showed that a combination of planetary close encounters and the ν₅ secular resonance is a route to Sun-grazing orbits at semi-major axis a < 2 au. Therefore, ν₅ could provide another significant mechanism for the delivery of asteroids impacting the Sun and the inner planets that acts well within the NEO region.

Figure 9 shows two test asteroids located near ν₂ and ν₃. These objects both show the corresponding critical argument librating around 0^° but, in contrast to the previous examples, these resonances do not increase the eccentricity in the 2 My integration time-span. On the other hand, the test asteroid affected by ν₃ shows an initial increasing of the inclination up to 30° caused by the effect of the resonance.

Here, we present only few examples corresponding to one of the maps of Fig. 1, however, all the numerical results obtained are in agreement with the locations computed for i₀ = 20^°. Therefore, we think that the maps obtained in Sects. 3.1 and 3.2 provide a reasonable estimation for the location of the secular resonances in the NEO region.

Fig. 7 Evolution of a test particle inside the ν6 secular resonance. The first three panels show the evolution of semi-major axis, eccentricity, and inclination, while the last panel shows the evolution of the critical angle ϖ − ϖS.

4 Discussion

4.1 Possible dynamical paths inside the NEO region

Farinella et al. (1994) showed that the ν6 secular resonance is very efficient in increasing the eccentricity of asteroids initially located in the main belt, causing e to pass from values near to 0 to values near to 1 in just ~0.5 Myr. Foschini et al. (2000), Gladman et al. (2000) showed some numerical example in which the ν₅ secular resonance increases the eccentricity of NEOs located at a < 2 au, causing asteroids to end up in a collision with the Sun. However, the location of the ν₅ resonance inside the NEO region was not determined. In Sect. 3.4, we also showed an example where ν₅ is able to pump up the eccentricity to values near to 1. Thus, some secular resonances are responsible for the delivery of asteroids impacting the planets or the Sun, even well inside the NEO region.

To search for possible dynamical paths, we took the population of known NEOs and identified those that are near the secular resonances. The nominal orbital elements were downloaded from the Near-Earth Objects Coordination Centre⁴ (NEOCC) on August 24, 2022. Then, we computed their proper frequencies ɡ, s, and their proper elements e₀, i₀, as we defined in Sect. 2. We plotted all the NEOs in the (a, e₀) and in the (a, i₀) planes and highlighted the objects near the v_j, j = 2, 3, 5, 6 and v₁₆ secular resonances, as well as Kozai librators. NEOs near v₄ were not plotted because they almost overlap with v₃. Since we did not attempt to compute the width of the resonances, we set an arbitrary threshold of 1" yr⁻¹ for the identification of resonant objects that we believe to be a reasonable choice (see also Gronchi & Milani 2001). Figure 10 shows the result of the plot. We note that secular resonances appear as a swarm of points in this figure, differently to what described in Sect. 3. This is because here we show the projections of the whole three-dimensional space of the proper elements on a two-dimensional plane, and also because we assumed a non-zero resonance width. From the plot in the (a, e₀)-plane, we can notice that objects near to v₆ cross the curve q = a₃, confirming again that this resonance is able to move objects onto Earth-crossing orbits (Bottke et al. 2002; Granvik et al. 2018). Moreover, NEOs near ν_j,·j = 2,…,6 extend to high eccentricity values, even surpassing 0.8. The fact that known NEOs are currently near these resonances and are placed at high eccentricity is an indication that all the v_j resonances may be able to pump up the eccentricity and deliver objects impacting the Sun. It is worth also noting that there are NEOs near v₁₆ between 1.5 au and 2 au that extend from the border of the NEO region at q = 1.3 au to q = a₃. Therefore, this resonance might also be able to move asteroids onto Earth-crossing orbits. We note that v₁₆ is near to the Hungaria region (see e.g. Froeschle & Scholl 1986), that is now recognized to be an important source region of NEOs (Granvik et al. 2018) producing a non-negligible fraction of Earth impactors.

Wetherill (1988) suggested that the v₁₆ secular resonance could be responsible for the production of NEOs with inclinations larger than 30°. Later, Gladman et al. (2000) showed with numerical simulations that this mechanism is actually possible. In the right panel of Fig. 10, we can see objects at high inclinations near v₂, v₅, and v₁₆, suggesting that they all could be capable of increasing the inclination. Indeed, the maps of Fig. 3 show that v₂ and v₅ extend towards values of the inclination larger than 50° for proper eccentricity of 0.6. Additionally, v₁₂, v₁₃, and v₁₆ all extend towards high inclination values, regardless of the value of the proper eccentricity. Therefore, even if not reported in Fig. 10, also v₁₂ and v₁₃ could be responsible for increasing the inclination to values higher than 30°.

In this work, we only gave an indication of possible dynamical paths inside the NEO region that are based on the secular motion and on the current position of known NEOs. Additional extensive numerical simulations with a full N-body model need to be performed to better investigate and establish these hypotheses.

Fig. 8 Evolution of a test particle inside the v5 secular resonance. The last panel shows the evolution of the critical argument, ϖ − ϖ.

Fig. 9 Evolution of a test particle in the v2 (red curve), and another one in the v3 (blue curve) secular resonances. The last panel shows the evolution of the corresponding critical argument σ: ϖ − ϖV (red curve) for v2, and ϖ − ϖE (blue curve) for v3.

Fig. 10 Known NEOs close to some secular resonances, in the phase space (a,e0) (left panel) and (a, i0) (right panel). Gray dots represent the background population of NEOs, while red, yellow, green, blue, and violet correspond to the v2, v3, v5, v6 and v16, respectively. Black dots denote objects in the Kozai resonance.

4.2 Effects of mean-motion resonances

The interaction between secular resonances and mean-motion resonances (MMRs) with Jupiter has been studied by Morbidelli & Moons (1993), Moons & Morbidelli (1995) in the main belt. The results presented in this paper were obtained with a semi-analytical model that does not take into account the effect of MMRs between the asteroid and a planet. One method for the computation of proper elements and frequencies of resonant NEOs, which extends that of Gronchi & Milani (2001), was recently developed by Fenucci et al. (2022). These authors showed that the secular evolution of asteroids that are placed well inside the strongest MMRs with Jupiter is significantly different from the non-resonant secular evolution, and that the frequencies, ɡ and s, are generally affected by these resonances. The same happens for some MMRs with the Earth, and Venus. On the contrary, MMRs with Mars have generally little effect on the long-term dynamics.

Some of the strongest MMRs with Jupiter, such as the 5:2, 7:3, and 2:1, are located beyond 2.5 au. The 3:1 MMR, that is the strongest one with Jupiter, is located at about 2.5 au. Therefore, they should not modify the results obtained for the nodal resonances $v_{1_j}$, since they always appear at a < 2.5 au, while they may somewhat affect the location of v₃, v₄, and v₆ at inclinations, i₀, larger than about 20°. On the other hand, the 7:2, 4:1, and 5:1 MMRs with Jupiter occur at smaller semi-major axis values, that is, at about 2.25 au, 2.06 au, and 1.78 au, respectively. Additionally, low-order resonances with Venus and the Earth mostly appear at a < 2.5 au (see e.g., Gallardo 2006), hence, they could all somewhat affect the secular resonances location.

Recently, Zhou et al. (2019), Xu et al. (2022) determined maps of secular resonances near the 1:1 MMR with the Earth and with Venus, computing the proper frequencies by using pure numerical integrations and frequency analysis (Laskar 1988, 1990, 2005). The maps obtained in these works do show some differences with respect to those presented in this paper. However, the effects of MMRs is only local and the secular resonances are affected only in a strip in semi-major axis of ~0.02 au width for the 1:1 MMR with the Earth and of ~0.008 au width for the 1:1 MMR with Venus.

To give a bigger picture of the strongest MMRs appearing at a < 2.5 au, we show their widths in Fig. 11. The width was computed with a simplified analytical method taken from Morbidelli (2002), for a fixed value of inclination of 15°. MMRs generally occupy a small volume of this part of the phase space and those with Venus and the Earth are especially narrow. This observation, coupled with the fact that MMRs significantly change the secular dynamics only when the asteroid is well inside the resonance (Fenucci et al. 2022), should limit the effect of MMRs on the results presented in this paper, that are aimed to give the general picture of the whole NEO region.

Fig. 11 Strongest mean-motion resonances at semi-major axis a < 2.5 au. The filled light-blue area represents the resonance width, computed with an analytic method (see Morbidelli 2002) for a fixed value of inclination of 15°.

5 Conclusions

In this paper, we determined the location of the secular resonances with the planets from Venus to Neptune inside the NEO region. Proper elements and proper frequencies of NEOs were computed by using a semi-analytical secular model that allow for the propagation of the secular evolution of objects beyond orbit-crossing singularities. Maps of the secular resonances were then obtained by keeping the proper inclination, i0, (or the proper eccentricity, e₀) at a fixed value, and we showed how the locations change by varying the value of the fixed proper element.

The v_j, j = 2,…,6 secular resonances all appear well inside the NEO region, especially at inclinations smaller than about 30^°. On the other hand, the $v_{1_j}$ secular resonances with the node of the planets mainly appear solely for j = 3, 4, 6, and they tend to be negligible at small inclinations. To confirm the computed locations, we performed full numerical N-body simulations that included all the planets of the Solar System. The results obtained showed that the maps computed with the semi-analytical approach provide a good estimate of the secular resonance locations. The current distribution of NEOs close to secular resonances suggests that some of these resonances could be able to produce Sun-grazing asteroids and Earth impactors -or they may increase the inclination to values bigger than 30^°.

Acknowledgements

The authors have been supported by the MSCA-ITN Stardust-R, Grant Agreement no. 813644 under the European Union H2020 research and innovation program. GFG also acknowledges the project MIUR-PRIN 20178CJA2B "New frontiers of Celestial Mechanics: theory and applications" and the GNFM-INdAM (Gruppo Nazionale per la Fisica Matematica).

References

All Tables

All Figures

Fig. 1 Map of secular resonances in the (a, e0)-plane for a between 0.5 and 3 au, for different values of proper inclination, i0. The left column shows the locations of vj, j = 2,3,4,5,6, while the right column shows the locations of $v_{1_j}$, j = 2,3,4,6. The color code of each resonance is indicated by the legend in the top row panels. From the top to the bottom, rows correspond to values of i0 equal to 40°, 30°, 20°, and 10°, respectively. The dashed curves correspond to q = a2, a3, a4 and Q = a2, a3, a4 where q = a(1 − e), Q = a(1 + e) are the perihelion and the aphelion distance of the asteroid, respectively. We also reported the location of the Kozai librators in grey, identified by g − s = 0.

In the text

	Fig. 2 Map of the vj, j = 2,3,4,5,6 secular resonances in the (a, e0)-plane for a between 0.5 and 3 au, at i0 = 5° (left panel), and at i0 = 2° (right panel). Dashed curves correspond to q = a2, a3, a4 and Q = a2, a3, a4 We also reported the location of the Kozai librators in grey, identified by ɡ − s = 0.
In the text

Fig. 3 Map of secular resonances in the (a, i0)-plane for a between 0.5 and 3 au, for different values of proper eccentricity, e0. The left column shows the locations of vj, j = 2, 3, 4, 5, 6, while the right column shows the locations of $v_{1_j}$, j = 2, 3, 4, 6. The color code of each resonance is indicated by the legend in the top row panels. From top to the bottom: rows correspond to values of e0 equal to 0.1, 0.2, 0.4, and 0.6, respectively. We also report the location of the Kozai librators in grey, identified by ɡ − s = 0.

In the text

Fig. 4 Map of secular resonances in the (a, i0)-plane for a between 0.5 and 3 au, for different values of proper eccentricity, e0. The left column shows the locations of vj, j = 2, 3, 4, 5, 6, while the right column shows the locations of $v_{1_j}$, j = 2, 3, 4, 6. The color code of each resonance is indicated by the legend in the top-left panel. From top to the bottom: Rows correspond to values of e0 equal to 0.05, 0.02, respectively. We also report the location of the Kozai librators in grey, identified by ɡ − s = 0.

In the text

	Fig. 5 Frequency, ɡ, (black solid curve) computed on a grid in proper eccentricity, e0, for a fixed value of semi-major axis of a = 1.2 au and proper inclination of i0 = 2°. The horizontal straight lines correspond to the values ɡj, j = 2, 3, 4, 5 of the proper frequencies of the planets. The vertical dashed lines correspond to the values of eccentricity, such that a(1 − e0) = aj, j = 2, 3.
In the text

	Fig. 6 Zoom on the map of the secular resonances vj, j = 2, 3, 4, 5, 6 at i0 = 2°, near a = 1.2au and e0 = 0.17.
In the text

	Fig. 7 Evolution of a test particle inside the ν6 secular resonance. The first three panels show the evolution of semi-major axis, eccentricity, and inclination, while the last panel shows the evolution of the critical angle ϖ − ϖS.
In the text

	Fig. 8 Evolution of a test particle inside the v5 secular resonance. The last panel shows the evolution of the critical argument, ϖ − ϖ.
In the text

	Fig. 9 Evolution of a test particle in the v2 (red curve), and another one in the v3 (blue curve) secular resonances. The last panel shows the evolution of the corresponding critical argument σ: ϖ − ϖV (red curve) for v2, and ϖ − ϖE (blue curve) for v3.
In the text

	Fig. 10 Known NEOs close to some secular resonances, in the phase space (a,e0) (left panel) and (a, i0) (right panel). Gray dots represent the background population of NEOs, while red, yellow, green, blue, and violet correspond to the v2, v3, v5, v6 and v16, respectively. Black dots denote objects in the Kozai resonance.
In the text

	Fig. 11 Strongest mean-motion resonances at semi-major axis a < 2.5 au. The filled light-blue area represents the resonance width, computed with an analytic method (see Morbidelli 2002) for a fixed value of inclination of 15°.
In the text