Maps of secular resonances in the NEO region

Context. From numerical simulations, it is known that some secular resonances may affect the motion of near-Earth objects (NEOs). However, the specific location of the secular resonance inside the NEO region is not fully known because the methods previously used to predict their location cannot be applied to highly eccentric orbits or the time when the NEOs cross the orbits of the planets. Aims. In this paper, we aim to map the secular resonances with the planets from Venus to Saturn in the NEO region, while including high eccentricity values as well. Methods. We used an averaged semi-analytical model that can deal with orbit-crossing singularities for the computation of the secular dynamics of NEOs, from which we were able to obtain suitable proper elements and proper frequencies. Then, we computed the proper frequencies over a uniform grid in the proper elements space. Secular resonances can thus be located by the level curves corresponding to the proper frequencies of the planets. Results. We determined the location of the secular resonances with the planets from Venus to Saturn, showing that they appear well within the NEO region. By using full numerical N -body simulations, we also showed that the location predicted by our method is fairly accurate. Finally, we provided some indications about possible dynamical paths inside the NEO region due to the presence of secular resonances.


Introduction
In the dynamics of the N-body problem, with N > 2, the geometry of the orbits of the bodies changes with 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 g j and s j , respectively, where the index j = 1, . . ., 8 labels the planets from Mercury to Neptune.Accurate values of planetary secular frequencies were computed, for instance, by Nobili et al. (1989) and Laskar et al. (2004Laskar et al. ( , 2011)).
The presence of the planets also cause both short and longterm 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 the one of a planet.Linear secular resonances involve only one asteroid's frequency and one planet's frequency, and resonances of type g ≃ g j are typically denoted with ν j , while those of type s ≃ s j are denoted with ν 1 j .
The problem of locating the secular resonances in the main belt underwent a series of improvements, and it is nowadays 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 asteroids' proper elements and proper frequencies.The method was later refined by Williams & Faulkner (1981) to better locate the three strongest secular resonances, namely ν 5 , ν 6 , and ν 16 .Later on, a simple analytical model was used by Yoshikawa (1987) to locate more precisely the ν 6 secular resonance.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, 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 attempt to locate secular resonances with 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 ν 5 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), that are asteroids with a perihelion distance q smaller than 1.3 au.So far, the first and only attempt was done 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 to better understand 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 ν 6 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 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 g and s are extracted.Proper frequencies are then computed on a grid of points in the phase spaces of (a, e 0 ) or (a, i 0 ), where e 0 are i 0 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 Section 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 Section 3 we show the location of the secular resonances obtained with our method, and show some comparison with full numerical simulations.In Section 4, we reported the implication of our results on the global dynamics of NEOs, and we discussed the effects of mean-motion resonances.Finally, in Section 5 we provide our conclusions.

The 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 denote the Gauss constant with k = √ Gm 0 , where G is the universal gravitational constant, and m 0 is the mass of the Sun.We also set µ j = m j /m 0 , j = 2, . . ., 6, where m j is the mass of the j-th planet.The canonical Delaunay variables are defined as 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 pericenter1 , and Ω is the longitude of the ascending node.The Hamiltonian describing the motion of the asteroid is where H 0 is the unperturbed Keplerian part, and H 1 accounts for the perturbations of the planets: In Eq. ( 3), r and r j denote the heliocentric position of the asteroid and the j-th 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, i.e.
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).Note that the indirect perturbation is not present in the averaged Hamiltonian H, since its double integral vanishes.In this model, L is a first integral, hence the semi-major 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 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 1-degree of freedom Hamiltonian depending only on (G, g).
The averaged Hamiltonian 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, g, z} one of the Delaunay elements, the derivatives ∂H/∂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 H 1 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 neighbourhood 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 permits to compute a solution of Hamilton's equations for H beyond an orbit crossing configuration.In this work we used 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).

Locating the secular resonances
Proper frequencies can be computed as in Gronchi & Milani (2001).A solution of the Hamiltonian secular system defined by H is computed numerically and, if ω circulates, we can find the times t 0 and t π/2 such that ω(t 0 ) = 0, ω(t π/2 ) = π/2.We also set Ω 0 = Ω(t 0 ), Ω π/2 = Ω(t π/2 ).The Hamiltonian H is even and π-periodic in ω, therefore the frequency g − s of the argument of perihelion ω is given by while the frequency s of the longitude of the node Ω is given by The dynamical state in which ω librates is also called Kozai resonance.In this work we assume that ω librates around 0 • , therefore we take t 0 such that ω(t 0 ) = 0 and ω(t 0 ) > 0.Then, we compute the time ∆t after which ω vanishes again and we denote by ∆Ω the change in Ω during this timespan.The frequency g − s of ω vanishes, but we can still define a libration frequency lf, and the frequency s for Ω, as The formula for the libration frequency lf is justified by the symmetry properties of H.Besides the proper frequencies, the proper eccentricity and the proper inclination can be defined as e 0 = e(t 0 ) and i 0 = i(t 0 ), 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 along a solution.However, the two definitions coincide in the case that the solution does not cross any planet's orbit (Michel & Froeschlé 1997).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 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 fix the initial value of the proper eccentricity e 0 (or the initial proper inclination i 0 ), and choose a grid in the (a, i 0 )-plane (or in the (a, e 0 )-plane).The proper frequencies g and s are computed for each point of the grid.The location of the secular resonances ν j , ν 1 j , j = 2, . . ., 6, in the (a, i 0 )-plane or in the (a, e 0 )-plane, is then given by the contour lines corresponding to the levels g = g 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), that are reported in Table 1.

Secular resonance maps at fixed inclination
We discretized the (a, e 0 )-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 < 0.8.The proper inclination was fixed, and we computed the maps for i The left column of Fig. 1 shows the location of the ν j , j = 2, . . ., 6 secular resonances, while the right column shows that of ν 1 j , j = 2, 3, 4, 6, for values of proper inclination i 0 equal to 40 • , 30 • , 20 • , and 10 • .Figure 2 shows the location of the ν j , j = 2, . . ., 6 resonances at lower inclination, namely for i 0 = 5 • , 2 • .Figure 1 and 2 also show the curves q = a 2 , a 3 , a 4 and Q = a 2 , a 3 , a 4 , 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 ν j , j = 2, 3, 4, 5 (i.e.resonances with the inner planets, and with Jupiter) always appear inside the NEO region, while the ν 6 secular resonance with Saturn appears to a large extent at proper inclinations i 0 smaller than 30 • .The resonances ν 3 with the Earth and ν 4 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 0 = 40 • , all the secular resonances are confined at values of eccentricity e 0 larger than 0.2, and at semimajor axis a larger than about 1.3 au.As the proper inclination i 0 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 0 = 30 • to i 0 = 20 • , there is a significant change in the location of the ν 2 and ν 5 resonances.For inclination larger than 30 • these two resonances are present in the region q < a 4 , and they are confined at rather low eccentricity, while for q > a 4 they extend towards larger eccentricity values, up to the limit of 0.8 we considered.For inclination smaller than 30 • they disappear from the region q < a 4 , 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 for inclination between 10 • and 40 • .Additionally, the ν 5 reso- The left column shows the locations of ν j , j = 2, 3, 4, 5, 6, while the right column shows the locations of ν 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 i 0 equal to 40 • , 30 • , 20 • , and 10 • , respectively.The dashed curves correspond to q = a 2 , a 3 , a 4 and Q = a 2 , a 3 , a 4 , 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.
nance shows a second branch at a < 1 au for i 0 = 20 • .The same happens to ν 3 and ν 4 for i 0 = 10 • .
The secular resonance maps for i 0 = 2 • and 5 • (see Fig. 2) become more complex.At i 0 = 5 • , the ν 3 and ν 4 disappear from the region q < a 4 , while a new almost-circular branch appears in the region a 3 < q < a 4 .In the region q > a 3 , these two resonances are close to the curve q = a 3 for semi-major axis smaller than about 1.7 au, at which they start extending towards larger values of eccentricity.A similar feature is seen for ν 2 and ν 5 , that tend to follow the curves q = a 3 and q = a 2 for certain intervals of semi-major axis values.On the contrary, the ν 6 resonance does not change much with respect to the one obtained at i 0 = 10 • .Note also 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 0 = 2 • is similar to the previous case, with the difference that the alignment along the curves q = a 2 and q = a 3 is more definite, and it appears also for the ν 6 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 g is significantly affected.Moreover, their role increases as the proper inclination decreases.
The right column of Fig. 1 shows the location of the resonances ν 1 j , j = 2, 3, 4, 6 involving the longitude of the node.The ν 12 resonance basically never appears in the region 0.5 au < a < 3 au for inclination i 0 up to 40 • .This was already noticed by Michel & Froeschlé (1997), who found that ν 12 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 should not play any major role in the dynamics of NEOs.The resonances ν 13 and ν 14 involving the Earth and Mars are again close to each other, since also the node of their orbits precesses at a similar rate (see Table 1).
At inclination i 0 = 40 • , ν 13 and ν 14 appear at a < 2 au, while ν 16 is present in the region a < 2.5 au.These three resonances all span a large interval of eccentricities, arriving up to e 0 ≈ 0.7.Moreover, ν 13 and ν 16 extend towards the limit of e 0 = 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 0 = 30 • , the qualitative picture is the same as that obtained at i 0 = 40 • .At i 0 = 20 • , the ν 16 resonance does not change its qualitative shape.On the other hand, ν 13 and ν 14 are composed by four different separated curves, confined at eccentricity smaller than about 0.5.Similarly to the case of the ν j resonances, the curves of constant q, Q equal to the semi-major axis of a planet play a role in determining the proper frequencies s, and the location of the secular resonances.At i 0 = 10 • , ν 16 has a different qualitative structure.It shows a closed curve in the region a 2 < q < a 3 , and another branch appears at q < a 3 .This is the only resonance appearing to a large extent, while all the other either disappeared or are confined to a very small region of the phase space.At lower values of proper inclination i 0 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 g − s = 0 or, equivalently, when the argument of the perihelion ω librates.In Fig. 1 and 2 the Kozai resonance region is identified in gray.At i 0 = 40 • Kozai librators are rare, and indeed the gray area is small.At i 0 = 30 • and 20 • , the Kozai resonance appears to a large extent in three different regions: near a = a 2 and a = a 3 at small eccentricity, and at the intersection between q = a 2 and Q = a 3 .At inclination i 0 ≤ 10 • only the Kozai librators near the semi-major axes of Venus and the Earth are left, while a small gray area appear near a = a 4 .
Note also that the extent of the gray area decreases as the proper inclination value decreases.

Secular resonance maps at fixed eccentricity
We discretized the (a, i 0 )-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 0 < 50 • .The proper eccentricity was fixed, and we computed the maps for values e 0 = 0.1 • h, h = 1, . . ., 7, and for low eccentricity e 0 = 0.02, 0.05.The left column of Fig. 3 shows the location of the ν j , j = 2, . . ., 6 secular resonances, while the right column shows the location of ν 1 j , j = 2, 3, 4, 6, for values of proper eccentricity e 0 equal to 0.1, 0.2, 0.4, and 0.6.Note that the first row of Fig. 3 is the same case considered by Michel & Froeschlé (1997).The results we obtained are in a 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 = 0.05, 0.02.
From the left column of Fig. 3, we can see that all the resonances ν j , j = 2, . . ., 6 appear in the region 0.5 au < a < 3 au.As expected, ν 3 and ν 4 are always close to each other.At proper eccentricity e 0 = 0.1, the ν 3 , ν 4 , and ν 6 resonances all have four different branches, each one appearing in the regions a > a 4 , a 3 < a < a 4 , a 2 < a < a 3 , and in a < a 2 .The ν 6 resonance with Saturn is confined at inclination smaller than about 20 • .The resonances ν 2 and ν 5 have three different branches, each one appearing for a > a 3 , a 2 < a < a 3 , and in a < a 2 .We can also notice the presence of Kozai librators around a = a 2 and a = a 3 , that fill this area of the phase space up to i ≈ 35 • .
At proper eccentricity e 0 = 0.2 the branches in the region a 2 < a < a 3 almost disappeared, while those that were appearing at a < a 2 are shifted to smaller values of semi-major axis, if compared with the previous case.For a > a 3 , the map is qualitatively similar to the case e 0 = 0.1.The Kozai librators still appear near a 2 and a 3 , but they fill a smaller area of the phase space, that 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 proper inclination of 30 • , and 40 • .We found that the location of the secular resonances for e 0 = 0.3 is overall qualitatively similar to e 0 = 0.2, and therefore it is not shown.
The maps for e 0 = 0.4 and e 0 = 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, that appear to a large extent at larger semi-major axis value, 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 disappeared, consistently with the results of Fig. 1 and  2.
The left column of Fig. 4 shows the maps obtained for e 0 = 0.05 and 0.02, and they are qualitatively similar to that obtained for e 0 = 0.1.We can see that the Kozai resonance region is smaller than for the case e 0 = 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 ν 1 j , j = 2, 3, 4, 6 secular resonances.As seen in Sec.3.1, the resonance ν 12 with Venus is negligible.At eccentricity e 0 = 0.1, ν 16 has three different branches, while ν 13 and ν 14 have only one.These three resonances are confined at semi-major axis smaller than 2.3 au, and they span a large interval of inclination values.Note that ν 13 , ν 14 , and ν 16 appear also in the Kozai resonance region.For shows the locations of ν j , j = 2, 3, 4, 5, 6, while the right column shows the locations of ν 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 e 0 equal to 0.1, 0.2, 0.4, and 0.6, respectively.We also reported the location of the Kozai librators in grey, identified by g − s = 0.
larger value of proper eccentricity, ν 16 also have a single branch.Moreover, as e 0 increases, their location move towards smaller values of semi-major axis, and they span a smaller interval of inclinations.Finally, the right column of Fig. 4 shows the location of the ν 1 j resonances for e 0 = 0.02 and 0.05.The maps are qualitatively similar to the case e 0 = 0.1, with the difference that the three resonances involving the Earth, Mars, and Saturn all extend towards larger semi-major axis values.

Frequencies near perihelion equal to planets' semi-major axes
Figures 1, 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 i 0 decreases towards 0.
To better understand the reason of this, 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 have a better picture of the behaviour of the frequency, we also used a step in eccentricity of 5 × 10 −4 , that is smaller than that used to produce Fig. 1 and 2. Figure 5 shows the frequency g computed for a = 1.2 au and i 0 = 2 • , as a function of the proper eccentricity e 0 .When e 0 approaches the critical values e ( j) 0 , j = 2, 3 for which a(1 − e ( j) 0 ) = a j , j = 2, 3, the frequency g 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 (3) 0 ≈ 0.1667 and e (2) 0 ≈ 0.3972 the frequency g has a minimum that is smaller than g 5 , that is the smallest planetary frequency.As e 0 tends to e (3) 0 from the right, the frequency g increases sharply, and it crosses all the planetary frequencies g j , j = 2, . . ., 6 in a small interval of proper eccentricities e 0 .The same happens as e 0 tends to e (2)  0 from the right, with the difference that only g 2 and g 5 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 g 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 Sec. 3.4).
The sharp increase of the frequency g near planets' perihelion curves causes also some small numerical errors in the global plot of the maps presented in Fig. 1 and 2. Because we use a discrete grid in the (a, e 0 )-plane, it may happen that a value larger than one of the planets' frequency g j is not achieved in the neighbourhood of one of such perihelion curves.As a consequence, the corresponding secular frequency ν j would not appear in the plot.This can be seen for instance in the maps at i 0 = 2 • , presented in Fig. 2. Figure 5 shows that all the five resonances ν j , j = 2, . . ., 6 are crossed in a neighbourhood of a = 1.2 au and e 0 ≈ 0.17, while this is not clear to 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 neighbourhood.To prove this, we recomputed the secular frequencies in the portion (a, e 0 ) ∈ [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 −4 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 ν j appear close to q = a 3 , even though ν 3 , ν 4 and ν 6 are barely distinguishable since they are very close to each other.Overall, we found these numerical effects to happen at low proper inclination i 0 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 Sec.3.1, that was obtained using a less dense grid.

Comparison with numerical integrations
We performed purely numerical integrations to confirm the secular resonance locations predicted by the theory in Sec.3.1 and 3.2.To account for the fact that the maps are computed in the proper elements space, we took initial orbital elements along the curves g = g j (or s = s j ) and g = g 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 mercury2 package (see also Fenucci & Novaković 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 Horizons3 ephemeris system.Note that here we take into account a complete model of the Solar System to show that the simplified model presented in Sec. 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 1 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, that is affected by the ν 6 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 the Sun in less than 0.5 Myr evolution.This is a well known effect of the ν 6 resonance, found for the first time by Farinella et al. (1994), that is able to deliver Sun impactors even starting from the main belt at very low eccentricity.This simulation shows that, as expected, this behaviour 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, ν 5 was able to push it to large eccentricity values, showing a behaviour similar to that of ν 6 .By using numerical simulations, (Gladman et al. 2000) showed that a combination of planetary close encounters and the ν 5 secular resonance is a route to Sun-grazing orbits at semi-major axis a < 2 au.Therefore, ν 5 could provide another significant mechanism for the delivery of asteroids impacting the Sun and the inner planets, that acts well inside the NEO region.
Figure 9 shows two test asteroids located near ν 2 and ν 3 .These objects both show the corresponding critical argument librating around 0 • but, contrary 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 ν 3 shows an initial increasing of the inclination up to 30 • caused by the effect of the resonance.Here we showed 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 0 = 20 • .Therefore, we think that the maps obtained in Sec 3.1 and 3.2 provide a reasonable estimation for the location of the secular resonances in the NEO region.

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 ν 5 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 ν 5 resonance inside the NEO region was not determined.In Sec.3.4, we also showed an example where ν 5 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 Centre4 (NEOCC) on 24 th August 2022.Then, we computed their proper frequencies g, s, and their proper elements e 0 , i 0 as we defined in Sec. 2. We plotted all the NEOs in the (a, e 0 ) and in the (a, i 0 ) planes, and highlighted the objects near the ν j , j = 2, 3, 5, 6 and ν 16 secular resonances, as well as Kozai librators.NEOs near ν 4 were not plotted, because they almost overlap with ν 3 .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.Note that secular resonances appear as a swarm of points in this figure, differently to what shown in Sec. 3.This is because here we show the projections of the whole 3-dimensional space of the proper elements on a 2-dimensional plane, and also because we assumed a non-zero resonance width.From the plot in the (a, e 0 )-plane, we can notice that objects near to ν 6 cross the curve q = a 3 , 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 up to large eccentricity values, even over 0.8.The fact that known NEOs are currently near these resonances, and are placed at high eccentricity, is an indication that all the ν 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 ν 16 between 1.5 au and 2 au, that extend from the border of the NEO region at q = 1.3 au to q = a 3 .Therefore, also this resonance could be able to move asteroids onto Earth-crossing orbits.Note that ν 16 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 ν 16 secular resonance could be responsible for the production of NEOs with inclination 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 ν 2 , ν 5 , and ν 16 , suggesting that they all could be capable of increasing the inclination.Indeed, the maps of Fig. 3 show that ν 2 and ν 5 extend towards values of the inclination larger than 50 • for proper eccentricity of 0.6.Additionally, ν 12 , ν 13 , and ν 16 all extend towards high inclination values, regardless of the value of the proper eccentricity.Therefore, even if not reported in Fig. 10, also ν 12 and ν 13 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.

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.A method for the computation of proper elements and frequencies of resonant NEOs, that extends the one by Gronchi & Milani (2001), has been recently developed by Fenucci et al. (2022).The 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 the frequencies g 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 longterm 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 ν 1 j , since they always appear at a < 2.5 au, while they may somewhat affect the location of ν 3 , ν 4 , and ν 6 at inclinations i 0 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, i.e. 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.  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(Laskar , 1990(Laskar , 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 width if 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.

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 permits to propagate the secular evolution of objects beyond orbitcrossing singularities.Maps of the secular resonances were then obtained by keeping the proper inclination i 0 (or the proper eccentricity e 0 ) at a fixed value, and we showed how the locations change by varying the value of the fixed proper element.
The ν 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 ν 1 j secular resonances with the node of the planets appear to a large extent only for j = 3, 4, 6, and they tend to be negligible at small inclination.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 semianalytical approach provide a good estimate of the secular reso-nance 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 larger than 30 • .

Fig. 1 .
Fig.1.Map of secular resonances in the (a, e 0 )-plane for a between 0.5 and 3 au, for different values of proper inclination i 0 .The left column shows the locations of ν j , j = 2, 3, 4, 5, 6, while the right column shows the locations of ν 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 i 0 equal to 40 • , 30 • , 20 • , and 10 • , respectively.The dashed curves correspond to q = a 2 , a 3 , a 4 and Q = a 2 , a 3 , a 4 , 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. 3 .
Fig.3.Map of secular resonances in the (a, i 0 )-plane for a between 0.5 and 3 au, for different values of proper eccentricity e 0 .The left column shows the locations of ν j , j = 2, 3, 4, 5, 6, while the right column shows the locations of ν 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 e 0 equal to 0.1, 0.2, 0.4, and 0.6, respectively.We also reported the location of the Kozai librators in grey, identified by g − s = 0.

Fig. 4 .Fig. 5 .
Fig.4.Map of secular resonances in the (a, i 0 )-plane for a between 0.5 and 3 au, for different values of proper eccentricity e 0 .The left column shows the locations of ν j , j = 2, 3, 4, 5, 6, while the right column shows the locations of ν 1 j , j = 2, 3, 4, 6.The color code of each resonance is indicated by the legend in the top-left panel.From the top to the bottom, rows correspond to values of e 0 equal to 0.05, 0.02, respectively.We also reported the location of the Kozai librators in grey, identified by g − s = 0.

Fig. 7 .
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 .

Fig. 8 .
Fig. 8. Same as Fig. 7, for a test particle in the ν 5 secular resonance.The last panel shows the evolution of the critical argument ̟ − ̟ J .

Fig. 11 .
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 • .