11institutetext: Department of Astronomy, Faculty of Mathematics, University
of Belgrade, Studentski trg 16, 11000 Belgrade, Serbia
11email<EMAIL_ADDRESS>22institutetext: ESA NEO Coordination
Centre, Largo Galileo Galilei, 1, 00044 Frascati, Italy 33institutetext:
Elecnor Deimos, Via Giuseppe Verdi, 6, 28060 San Pietro Mosezzo, Italy
44institutetext: Dipartimento di Matematica, Università di Pisa, Largo B.
Pontecorvo 5, 56127 Pisa, Italy

# Maps of secular resonances in the NEO region

Marco Fenucci 112233 Giovanni F. Gronchi 44 Bojan Novaković 11

(Received — / Accepted —)

###### Abstract

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 can not be used
for highly eccentric orbits and 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, even for high values of the
eccentricity.

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 can obtain suitable proper elements and proper frequencies.
Then, we computed the proper frequencies over a uniform grid in the proper
elements space. Secular resonances are thus 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 inside 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.

###### Key Words.:

Celestial mechanics - Minor planets, asteroids: general

##  1 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 $\varpi$
and of the longitude of the ascending node $\Omega$ of the orbits of the
planets are typically denoted with $g_{j}$ and $s_{j}$, respectively, where
the index $j=1,\dots,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. (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
the one of a planet. Linear secular resonances involve only one asteroid’s
frequency and one planet’s frequency, and resonances of type $g\simeq g_{j}$
are typically denoted with $\nu_{j}$, while those of type $s\simeq s_{j}$ are
denoted with $\nu_{1j}$.

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 $\nu_{5},\nu_{6}$, and $\nu_{16}$. Later on, a
simple analytical model was used by Yoshikawa (1987) to locate more precisely
the $\nu_{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 $\nu_{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 $\nu_{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.

##  2 Methods

###  2.1 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 $\mathtt{k}=\sqrt{\mathcal{G}m_{0}}$, where $\mathcal{G}$ is the
universal gravitational constant, and $m_{0}$ is the mass of the Sun. We also
set $\mu_{j}=m_{j}/m_{0},\ j=2,\dots,6$, where $m_{j}$ is the mass of the
$j$-th planet. The canonical Delaunay variables are defined as

$\begin{split}L=\mathtt{k}\sqrt{a},&\qquad\ell=M,\\\
G=\mathtt{k}\sqrt{a(1-e^{2})},&\qquad g=\omega,\\\
Z=\mathtt{k}\sqrt{a(1-e^{2})}\cos i,&\qquad z=\Omega,\end{split}$ (1)

where $a$ is the semi-major axis, $e$ is the eccentricity, $i$ is the
inclination, $M$ is the mean anomaly, $\omega$ is the argument of the
pericenter111Note that the symbol $g$ is also used to denote the secular
frequency of the longitude of perihelion $\varpi$, as mentioned in the
Introduction. However, it is clear from the context when we use $g$ to refer
to the argument of the pericenter $\omega$, or to refer to the frequency of
$\varpi$., and $\Omega$ is the longitude of the ascending node. The
Hamiltonian describing the motion of the asteroid is

$\mathcal{H}=\mathcal{H}_{0}+\varepsilon\mathcal{H}_{1},\qquad\varepsilon=\mu_{5},$
(2)

where $\mathcal{H}_{0}$ is the unperturbed Keplerian part, and
$\mathcal{H}_{1}$ accounts for the perturbations of the planets:

$\begin{split}\mathcal{H}_{0}&=-\frac{\mathtt{k}^{4}}{2L^{2}},\\\
\mathcal{H}_{1}&=-\mathtt{k}^{2}\sum_{j=2}^{6}\frac{\mu_{j}}{\mu_{5}}\bigg{(}\frac{1}{|\mathbf{r}-\mathbf{r}_{j}|}-\frac{\mathbf{r}\cdot\mathbf{r}_{j}}{|\mathbf{r}_{j}|^{3}}\bigg{)}.\end{split}$
(3)

In Eq. (3), $\mathbf{r}$ and $\mathbf{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 $\ell$ and $\ell_{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.

$\overline{\mathcal{H}}=-\varepsilon\mathtt{k}^{2}\sum_{j=2}^{6}\frac{\mu_{j}}{\mu_{5}}\int_{0}^{2\pi}\int_{0}^{2\pi}\frac{1}{|\mathbf{r}-\mathbf{r}_{j}|}\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). Note that the indirect perturbation is not
present in the averaged Hamiltonian $\overline{\mathcal{H}}$, since its double
integral vanishes. In this model, $L$ is a first integral, hence the semi-
major axis $a$ remains constant. Moreover, $\Omega$ 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{\mathcal{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 $\overline{\mathcal{H}}$ has a first order polar
singularity when the orbit of the asteroid crosses that of a planet (Gronchi &
Milani, 1998). By denoting with $\smash{y\in\\{G,Z,g,z\\}}$ one of the
Delaunay elements, the derivatives $\partial\overline{\mathcal{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 $\mathcal{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 $\overline{\mathcal{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).

###  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{\mathcal{H}}$ is
computed numerically and, if $\omega$ circulates, we can find the times
$t_{0}$ and $t_{\pi/2}$ such that $\omega(t_{0})=0,\ \omega(t_{\pi/2})=\pi/2$.
We also set $\Omega_{0}=\Omega(t_{0}),\ \Omega_{\pi/2}=\Omega(t_{\pi/2})$. The
Hamiltonian $\overline{\mathcal{H}}$ is even and $\pi$-periodic in $\omega$,
therefore the frequency $g-s$ of the argument of perihelion $\omega$ is given
by

$g-s=\frac{2\pi}{4(t_{\pi/2}-t_{0})},$ (5)

while the frequency $s$ of the longitude of the node $\Omega$ is given by

$s=\frac{\Omega_{\pi/2}-\Omega_{0}}{t_{\pi/2}-t_{0}}.$ (6)

The dynamical state in which $\omega$ librates is also called Kozai resonance.
In this work we assume that $\omega$ librates around $0^{\circ}$, therefore we
take $t_{0}$ such that $\omega(t_{0})=0$ and $\dot{\omega}(t_{0})>0$. Then, we
compute the time $\Delta t$ after which $\omega$ vanishes again and we denote
by $\Delta\Omega$ the change in $\Omega$ during this timespan. The frequency
$g-s$ of $\omega$ vanishes, but we can still define a libration frequency lf,
and the frequency $s$ for $\Omega$, as

$\text{lf}=\pm\frac{2\pi}{2\Delta t},\quad s=\frac{\Delta\Omega}{\Delta t}.$
(7)

The formula for the libration frequency lf is justified by the symmetry
properties of $\overline{\mathcal{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 $\omega$ 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 $\Omega$ is set to 0, but
any other value would not change the result since $\Omega$ does not appear in
the Hamiltonian $\overline{\mathcal{H}}$. The initial argument of pericenter
$\omega$ 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 $\nu_{j},\,\nu_{1j},\ j=2,\dots,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,\dots,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.

Table 1: Secular frequencies of the planets from Venus to Saturn, as determined by Laskar et al. (2011). Planet |  $g$ (${}^{\prime\prime}/$yr) |  $s$ (${}^{\prime\prime}/$yr)
---|---|---
Venus | 7.453 | -7.06
Earth | 17.368 | -18.848
Mars | 17.916 | -17.751
Jupiter | 4.257482 | 0.0
Saturn | 28.2449 | -26.347841

##  3 Results

###  3.1 Secular resonance maps at fixed inclination

Figure 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 $\nu_{j},j=2,3,4,5,6$, while the right column
shows the locations of $\nu_{1j},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^{\circ},30^{\circ},20^{\circ}$, and $10^{\circ}$, 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$.  Figure 2: Map of the $\nu_{j},j=2,3,4,5,6$
secular resonances in the $(a,e_{0})$-plane for $a$ between 0.5 and 3 au, at
$i_{0}=5^{\circ}$ (left panel), and at $i_{0}=2^{\circ}$ (right panel). Dashed
curves correspond to $q=a_{2},a_{3},a_{4}$ and $Q=a_{2},a_{3},a_{4}$. We also
reported the location of the Kozai librators in grey, identified by $g-s=0$.

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_{0}=2^{\circ},5^{\circ},10^{\circ},20^{\circ},30^{\circ},$ and
$40^{\circ}$.

The left column of Fig. 1 shows the location of the $\nu_{j},j=2,\dots,6$
secular resonances, while the right column shows that of $\nu_{1j},j=2,3,4,6$,
for values of proper inclination $i_{0}$ equal to
$40^{\circ},30^{\circ},20^{\circ}$, and $10^{\circ}$. Figure 2 shows the
location of the $\nu_{j},j=2,\dots,6$ resonances at lower inclination, namely
for $i_{0}=5^{\circ},2^{\circ}$. Figure 1 and 2 also show the curves
$\smash{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 $\nu_{j},\,j=2,3,4,5$ (i.e.
resonances with the inner planets, and with Jupiter) always appear inside the
NEO region, while the $\nu_{6}$ secular resonance with Saturn appears to a
large extent at proper inclinations $i_{0}$ smaller than $30^{\circ}$. The
resonances $\nu_{3}$ with the Earth and $\nu_{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^{\circ}$, all the
secular resonances are confined at values of eccentricity $e_{0}$ larger than
0.2, and at semi-major 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^{\circ}$ to
$i_{0}=20^{\circ}$, there is a significant change in the location of the
$\nu_{2}$ and $\nu_{5}$ resonances. For inclination larger than $30^{\circ}$
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^{\circ}$ 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^{\circ}$ and
$40^{\circ}$. Additionally, the $\nu_{5}$ resonance shows a second branch at
$a<1$ au for $i_{0}=20^{\circ}$. The same happens to $\nu_{3}$ and $\nu_{4}$
for $i_{0}=10^{\circ}$.

The secular resonance maps for $\smash{i_{0}=2^{\circ}}$ and 5∘ (see Fig. 2)
become more complex. At $i_{0}=5^{\circ}$, the $\nu_{3}$ and $\nu_{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 $\nu_{2}$ and $\nu_{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 $\nu_{6}$ resonance does not change
much with respect to the one obtained at $\smash{i_{0}=10^{\circ}}$. 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^{\circ}$ 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 $\nu_{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
$\smash{\nu_{1j},\,j=2,3,4,6}$ involving the longitude of the node. The
$\nu_{12}$ resonance basically never appears in the region 0.5 au $<a<$ 3 au
for inclination $i_{0}$ up to $40^{\circ}$. This was already noticed by Michel
& Froeschlé (1997), who found that $\nu_{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
$\nu_{13}$ and $\nu_{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^{\circ}$, $\nu_{13}$ and $\nu_{14}$ appear at $a<2$
au, while $\nu_{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}\approx 0.7$. Moreover, $\nu_{13}$ and $\nu_{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^{\circ}$, the qualitative picture is the same as that obtained at
$i_{0}=40^{\circ}$. At $i_{0}=20^{\circ}$, the $\nu_{16}$ resonance does not
change its qualitative shape. On the other hand, $\nu_{13}$ and $\nu_{14}$ are
composed by four different separated curves, confined at eccentricity smaller
than about 0.5. Similarly to the case of the $\nu_{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^{\circ}$, $\nu_{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 $\omega$
librates. In Fig. 1 and 2 the Kozai resonance region is identified in gray. At
$i_{0}=40^{\circ}$ Kozai librators are rare, and indeed the gray area is
small. At $i_{0}=30^{\circ}$ and $20^{\circ}$, 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}\leq 10^{\circ}$ 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.

###  3.2 Secular resonance maps at fixed eccentricity

Figure 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 $\nu_{j},j=2,3,4,5,6$, while the right column
shows the locations of $\nu_{1j},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$. Figure 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
$\nu_{j},j=2,3,4,5,6$, while the right column shows the locations of
$\nu_{1j},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$.

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^{\circ}<i_{0}<50^{\circ}$. The proper eccentricity
was fixed, and we computed the maps for values $e_{0}=0.1\cdot h,\
h=1,\dots,7$, and for low eccentricity $e_{0}=0.02,0.05$. The left column of
Fig. 3 shows the location of the $\nu_{j},j=2,\dots,6$ secular resonances,
while the right column shows the location of $\nu_{1j},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
$\nu_{j},\,j=2,\dots,6$ appear in the region 0.5 au $<a<$ 3 au. As expected,
$\nu_{3}$ and $\nu_{4}$ are always close to each other. At proper eccentricity
$e_{0}=0.1$, the $\nu_{3},\nu_{4}$, and $\nu_{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 $\nu_{6}$ resonance
with Saturn is confined at inclination smaller than about 20∘. The resonances
$\nu_{2}$ and $\nu_{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\approx 35^{\circ}$.

At proper eccentricity $\smash{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^{\circ}$, and $40^{\circ}$. 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\lesssim 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 $\nu_{1j},j=2,3,4,6$
secular resonances. As seen in Sec. 3.1, the resonance $\nu_{12}$ with Venus
is negligible. At eccentricity $e_{0}=0.1$, $\nu_{16}$ has three different
branches, while $\nu_{13}$ and $\nu_{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 $\nu_{13},\nu_{14}$,
and $\nu_{16}$ appear also in the Kozai resonance region. For larger value of
proper eccentricity, $\nu_{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 $\nu_{1j}$ 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.

###  3.3 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\times 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^{\circ}$, as a function of the proper eccentricity $e_{0}$. When
$e_{0}$ approaches the critical values $e_{0}^{(j)},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^{\circ}$. This feature is caused by the model itself, and it is due to the
orbit crossing singularity.

Figure 5: Frequency $g$ (black solid curve) computed on a grid in proper
eccentricity $e_{0}$, for a fixed value of semi-major axis $a=1.2$ au and
proper inclination $i_{0}=2^{\circ}$. The horizontal straight lines correspond
to the values $g_{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-e_{0})=a_{j},j=2,3$. Figure 6: Zoom of the map of the secular resonances
$\nu_{j},j=2,3,4,5,6$ at $i_{0}=2^{\circ}$, near $a=1.2$ au and $e_{0}=0.17$.

Figure 5 shows also an interesting feature. Between $\smash{e_{0}^{(3)}\approx
0.1667}$ and $\smash{e_{0}^{(2)}\approx 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_{0}^{(3)}$ from the right, the frequency $g$
increases sharply, and it crosses all the planetary frequencies
$\smash{g_{j},j=2,\dots,6}$ in a small interval of proper eccentricities
$e_{0}$. The same happens as $e_{0}$ tends to $e_{0}^{(2)}$ 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 $\nu_{j}$ would not appear in
the plot. This can be seen for instance in the maps at $i_{0}=2^{\circ}$,
presented in Fig. 2. Figure 5 shows that all the five resonances
$\nu_{j},j=2,\dots,6$ are crossed in a neighbourhood of $a=1.2$ au and
$e_{0}\approx 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})\in[1.15\text{ au},1.3\text{ au}]\times[0.1,0.3]$,
using a step of 0.004 au in semi-major axis and of $5\times 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 $\nu_{j}$ appear close to $q=a_{3}$, even though $\nu_{3},\nu_{4}$
and $\nu_{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.

###  3.4 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}\pm
1^{\prime\prime}/\text{yr}$ (or $s=s_{j}\pm 1^{\prime\prime}/\text{yr}$), for
a selected planet $j\in\\{2,3,4,5,6\\}$. The initial argument of pericenter
$\omega$ and the initial argument of the node $\Omega$ were both set to
$0^{\circ}$, while the initial mean anomaly $\ell$ was chosen randomly between
$0^{\circ}$ and $360^{\circ}$. 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
mercury222https://gemelli.spacescience.org/~hahnjm/software.html,
https://github.com/Fenu24/mercury 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
Horizons333https://ssd.jpl.nasa.gov/?horizons 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^{\circ}$. Figure 7 shows the evolution of a test asteroid selected with
the method described above, that is affected by the $\nu_{6}$ secular
resonance. The critical angle $\varpi-\varpi_{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 $\nu_{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 7: Evolution of a test particle inside the $\nu_{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 $\varpi-\varpi_{S}$.

Figure 8 shows a test asteroid affected by the $\nu_{5}$ secular resonance.
The critical angle $\varpi-\varpi_{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,
$\nu_{5}$ was able to push it to large eccentricity values, showing a
behaviour similar to that of $\nu_{6}$. By using numerical simulations,
(Gladman et al., 2000) showed that a combination of planetary close encounters
and the $\nu_{5}$ secular resonance is a route to Sun-grazing orbits at semi-
major axis $a<2$ au. Therefore, $\nu_{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 8: Same as Fig. 7, for a test particle in the $\nu_{5}$ secular
resonance. The last panel shows the evolution of the critical argument
$\varpi-\varpi_{J}$.

Figure 9 shows two test asteroids located near $\nu_{2}$ and $\nu_{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 $\nu_{3}$ shows an initial increasing of the inclination
up to $30^{\circ}$ caused by the effect of the resonance.

Figure 9: Same as Fig. 7, for two test particles in the $\nu_{2}$ (red curve)
and the $\nu_{3}$ (blue curve) secular resonances. The last panel shows the
evolution of the corresponding critical argument $\sigma$: $\varpi-\varpi_{V}$
(red curve) for $\nu_{2}$, and $\varpi-\varpi_{E}$ (blue curve) for $\nu_{3}$.

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^{\circ}$. 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.

##  4 Discussion

###  4.1 Possible dynamical paths inside the NEO region

Figure 10: Known NEOs close to some secular resonances, in the phase space
$(a,e_{0})$ (left panel) and $(a,i_{0})$ (right panel). Gray dots represent
the background population of NEOs, while red, yellow, green, blue, and violet
correspond to the $\nu_{2},\nu_{3},\nu_{5},\nu_{6}$ and $\nu_{16}$,
respectively. Black dots denote objects in the Kozai resonance.

Farinella et al. (1994) showed that the $\nu_{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 $\sim$0.5 Myr. Foschini et al. (2000); Gladman et al. (2000) showed some
numerical example in which the $\nu_{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 $\nu_{5}$ resonance
inside the NEO region was not determined. In Sec. 3.4, we also showed an
example where $\nu_{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
Centre444https://neo.ssa.esa.int/ (NEOCC) on 24${}^{\text{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 $\nu_{j},j=2,3,5,6$ and $\nu_{16}$ secular resonances, as well as Kozai
librators. NEOs near $\nu_{4}$ were not plotted, because they almost overlap
with $\nu_{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 $\nu_{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 $\nu_{j},j=2,\dots,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
$\nu_{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
$\nu_{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 $\nu_{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 $\nu_{16}$ secular resonance could be
responsible for the production of NEOs with inclination larger than
$30^{\circ}$. 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 $\nu_{2},\nu_{5}$, and $\nu_{16}$,
suggesting that they all could be capable of increasing the inclination.
Indeed, the maps of Fig. 3 show that $\nu_{2}$ and $\nu_{5}$ extend towards
values of the inclination larger than 50∘ for proper eccentricity of 0.6.
Additionally, $\nu_{12},\nu_{13}$, and $\nu_{16}$ all extend towards high
inclination values, regardless of the value of the proper eccentricity.
Therefore, even if not reported in Fig. 10, also $\nu_{12}$ and $\nu_{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.

###  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. 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
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 $\nu_{1j}$, since they always appear at
$a<2.5$ au, while they may somewhat affect the location of $\nu_{3},\nu_{4},$
and $\nu_{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.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 $\sim$0.02 au width for the 1:1 MMR with the
Earth, and of $\sim$0.008 au width for the 1:1 MMR with Venus.

Figure 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^{\circ}$.

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.

##  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 permits to propagate the secular evolution of objects beyond orbit-
crossing 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 $\nu_{j},j=2,\dots,6$ secular resonances all appear well inside the NEO
region, especially at inclinations smaller than about $30^{\circ}$. On the
other hand, the $\nu_{1j}$ 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 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
larger than 30∘.

###### Acknowledgements.

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

## References

  * Bottke et al. (2002) Bottke, W. F., Morbidelli, A., Jedicke, R., et al. 2002, Icarus, 156, 399
  * Chambers (1999) Chambers, J. E. 1999, MNRAS, 304, 793
  * Farinella et al. (1994) Farinella, P., Froeschle, C., Froeschle, C., et al. 1994, Nature, 371, 315
  * Fenucci et al. (2022) Fenucci, M., Gronchi, G. F., & Saillenfest, M. 2022, Celestial Mechanics and Dynamical Astronomy, 134, 23
  * Fenucci & Novaković (2022) Fenucci, M. & Novaković, B. 2022, Serbian Astronomical Journal, 204, 51
  * Foschini et al. (2000) Foschini, L., Farinella, P., Froeschlé, C., et al. 2000, A&A, 353, 797
  * Froeschle & Scholl (1986) Froeschle, C. & Scholl, H. 1986, A&A, 166, 326
  * Froeschle & Scholl (1989) Froeschle, C. & Scholl, H. 1989, Celestial Mechanics and Dynamical Astronomy, 46, 231
  * Gallardo (2006) Gallardo, T. 2006, Icarus, 184, 29
  * Gallardo et al. (2012) Gallardo, T., Hugo, G., & Pais, P. 2012, Icarus, 220, 392
  * Gladman et al. (2000) Gladman, B., Michel, P., & Froeschlé, C. 2000, Icarus, 146, 176
  * Granvik et al. (2018) Granvik, M., Morbidelli, A., Jedicke, R., et al. 2018, Icarus, 312, 181
  * Granvik et al. (2017) Granvik, M., Morbidelli, A., Vokrouhlický, D., et al. 2017, A&A, 598, A52
  * Gronchi & Milani (1998) Gronchi, G. F. & Milani, A. 1998, Celestial Mechanics and Dynamical Astronomy, 71, 109
  * Gronchi & Milani (2001) Gronchi, G. F. & Milani, A. 2001, Icarus, 152, 58
  * Gronchi & Tardioli (2013) Gronchi, G. F. & Tardioli, C. 2013, Discrete and Continuous Dynamical Systems - B, 18, 1323
  * Knežević (2020) Knežević, Z. 2020, MNRAS, 497, 4921
  * Knežević (2022) Knežević, Z. 2022, Serbian Astronomical Journal, 1
  * Knežević & Milani (2000) Knežević, Z. & Milani, A. 2000, Celestial Mechanics and Dynamical Astronomy, 78, 17
  * Knežević & Milani (2003) Knežević, Z. & Milani, A. 2003, A&A, 403, 1165
  * Knežević & Milani (2019) Knežević, Z. & Milani, A. 2019, Celestial Mechanics and Dynamical Astronomy, 131, 27
  * Knežević et al. (1991) Knežević, Z., Milani, A., Farinella, P., Froeschle, C., & Froeschle, C. 1991, Icarus, 93, 316
  * Kozai (1962) Kozai, Y. 1962, AJ, 67, 591
  * Laskar (1988) Laskar, J. 1988, A&A, 198, 341
  * Laskar (1990) Laskar, J. 1990, Icarus, 88, 266
  * Laskar (2005) Laskar, J. 2005, Frequency map analysis and quasiperiodic decompositions (Cambridge Scientific Pub), 99–134
  * Laskar et al. (2011) Laskar, J., Fienga, A., Gastineau, M., & Manche, H. 2011, A&A, 532, A89
  * Laskar et al. (2004) Laskar, J., Robutel, P., Joutel, F., et al. 2004, A&A, 428, 261
  * Michel & Froeschlé (1997) Michel, P. & Froeschlé, C. 1997, Icarus, 128, 230
  * Milani & Knežević (1990) Milani, A. & Knežević, Z. 1990, Celestial Mechanics and Dynamical Astronomy, 49, 347
  * Milani & Knežević (1994) Milani, A. & Knežević, Z. 1994, Icarus, 107, 219
  * Moons & Morbidelli (1995) Moons, M. & Morbidelli, A. 1995, Icarus, 114, 33
  * Morbidelli (2002) Morbidelli, A. 2002, Modern celestial mechanics : aspects of solar system dynamics (Taylor and Francis)
  * Morbidelli & Henrard (1991a) Morbidelli, A. & Henrard, J. 1991a, Celestial Mechanics and Dynamical Astronomy, 51, 131
  * Morbidelli & Henrard (1991b) Morbidelli, A. & Henrard, J. 1991b, Celestial Mechanics and Dynamical Astronomy, 51, 169
  * Morbidelli & Moons (1993) Morbidelli, A. & Moons, M. 1993, Icarus, 102, 316
  * Nobili et al. (1989) Nobili, A. M., Milani, A., & Carpino, M. 1989, A&A, 210, 313
  * Saillenfest et al. (2016) Saillenfest, M., Fouchard, M., Tommei, G., & Valsecchi, G. B. 2016, Celestial Mechanics and Dynamical Astronomy, 126, 369
  * Scholl & Froeschle (1986) Scholl, H. & Froeschle, C. 1986, A&A, 170, 138
  * Thomas & Morbidelli (1996) Thomas, F. & Morbidelli, A. 1996, Celestial Mechanics and Dynamical Astronomy, 64, 209
  * Tsirvoulis & Novaković (2016) Tsirvoulis, G. & Novaković, B. 2016, Icarus, 280, 300
  * Wetherill (1988) Wetherill, G. W. 1988, Icarus, 76, 1
  * Williams (1969) Williams, J. G. 1969, PhD thesis, University Of California, Los Angeles
  * Williams & Faulkner (1981) Williams, J. G. & Faulkner, J. 1981, Icarus, 46, 390
  * Xu et al. (2022) Xu, Y.-B., Zhou, L., Lhotka, C., Zhou, L.-Y., & Ip, W.-H. 2022, arXiv e-prints, arXiv:2206.06114
  * Yoshikawa (1987) Yoshikawa, M. 1987, Celestial Mechanics, 40, 233
  * Zhou et al. (2019) Zhou, L., Xu, Y.-B., Zhou, L.-Y., Dvorak, R., & Li, J. 2019, A&A, 622, A97


# AET-EFN: A Versatile Design for Static and Dynamic Event-Based Vision

Chang Liu  Xiaojuan Qi  Edmund Lam  Ngai Wong
The University of Hong Kong
<EMAIL_ADDRESS>{xjqi, elam<EMAIL_ADDRESS>
###### Abstract

The neuromorphic event cameras, which capture the optical changes of a scene,
have drawn increasing attention due to their high speed and low power
consumption. However, the event data are noisy, sparse, and nonuniform in the
spatial-temporal domain with extremely high temporal resolution, making it
challenging to design backend algorithms for event-based vision. Existing
methods encode events into point-cloud-based or voxel-based representations,
but suffer from noise and/or information loss. Additionally, there is little
research that systematically studies how to handle static and dynamic scenes
with one universal design for event-based vision. This work proposes the
Aligned Event Tensor (AET) as a novel event data representation, and a neat
framework called Event Frame Net (EFN), which enables our model for event-
based vision under static and dynamic scenes. The proposed AET and EFN are
evaluated on various datasets, and proved to surpass existing state-of-the-art
methods by large margins. Our method is also efficient and achieves the
fastest inference speed among others.

##  1 Introduction

Recently, the neuromorphic event cameras have drawn increasing attention from
researchers. Unlike traditional cameras, which take a whole frame
periodically, event cameras asynchronously capture the per-pixel optical
changes instead of the absolute brightness. Event cameras enjoy the advantages
of low latency (in microsecond level), high dynamic range, low power
consumption, and low bandwidth cost [10]. Because of the appealing
characteristics of event cameras, event-based vision is gaining more attention
in many research areas and applications, such as robotics, image deblurring
and auto-driving [3].

The outputs of event cameras are event flows defined by $(x,y,t,p)$, which are
pixel coordinates, timestamp, and polarity. The event flows denote
instantaneous increases or decreases of the light intensity, and are
incompatible with typical practices such as convolutional neural networks
(CNNs). Besides the unique data structure, event data are challenging to
process also because they do not contain information about the absolute
intensity but many noises. Therefore, to fully leverage the advantages, new
algorithms that can efficiently process event data are highly desired [3].

(a) Downsampling Event Data with FPS algorithm. The outlier noises are kept
during downsampling and the shape collapses.

(b) Events accumulated for 1/50, 1/100, 1/500 total time. Proper duration of
accumulation is needed for clear image.

Figure 1: Visualization of event data samples.

To leverage event cameras for vision tasks, some existing approaches regard
the event data as point clouds [27, 15], and directly apply point-cloud-based
models. However, the point-cloud-based methods such as [15, 11] originally
handle small clouds with around 1024 points, hence are not efficient for large
clouds such as event data with hundreds of thousands of events per second.
Moreover, the downsampling algorithms, such as the Furthest Point Sampling
used in the point-cloud-based methods [16, 11], are prone to noises commonly
found in event data, as shown in Fig. 1(a). Consequently, point-cloud-based
methods have limitations in accuracy and efficiency for event-based vision.

Besides point-cloud-based methods, the voxel-based approaches compress the
event data into voxels, bridging event cameras with conventional CNN
techniques [4, 9, 21, 2]. The voxel-based approaches such as [4] and [2]
achieve the state-of-the-art (SOTA) performance on several event-based
datasets, while enjoying higher efficiency than point-cloud-based methods, due
to the good regularity for parallel computing on GPUs.

Nonetheless, the voxel-based methods have to compress the data drastically to
meet the memory and computational requirements, and hence suffer from
information loss. According to the observation in our early-stage experiments,
the reduced temporal resolution also causes problems such as motion blur and
accuracy drop for dynamic scenes. Moreover, to fit the neural networks
originally designed for image data, they usually stack the different frames
along the channel axis [4, 2], which makes the temporal information tangle
with channel information. Therefore, voxel-based methods are limited by
information loss and flexibility under scenes at different speeds.

To this end, we propose an efficient representation for event data, which can
alleviate the information loss problem of voxel-based methods, while also
being robust under both static and dynamic scenes. We name it Aligned Event
Tensor (AET). In AET, the event data are quantized and voxelized in high
temporal resolution to minimize information loss. The voxels are then
compressed in local spatial-temporal regions by convolutional layers to
alleviate the misalignment problem across different frames. We also present a
two-branch neural network design, the Event Frame Net (EFN), which is
versatile for static and dynamic event-based vision. In EFN, a CNN encodes the
AET frames into vector representations, which are then processed independently
for frame-level predictions, or as a whole sequence for video-level
prediction.

Both AET and EFN can significantly increase the prediction accuracy for event-
based vision in the experiments. They are lightweight and enjoy the fastest
inference speed among the fellows. Our main contributions are:

  * •

We propose the efficient AET representation for event data, which improves
speed and accuracy for event data classification tasks under different speed
scenes.

  * •

The Event Frame Net (EFN) can process event data under static and dynamic
scenes, achieving the SOTA performances in all classification experiments.

  * •

Detailed experiments and speed analysis are conducted. Our method
significantly improves the accuracy while being the fastest among existing
methods. Visualization also empirically justifies our design.

##  2 Related Works

###  2.1 Reconstructive approaches

Several studies such as E2VID and EventSR [18, 24] reconstruct events back to
normal images for conventional vision. The reconstruction methods usually
utilize recursive modules such as ConvLSTM [26] or structures like UNet [19]
to reconstruct verisimilar images based on the events [18, 20]. However, the
massive pretraining workload is costly, and the cumbersome two-stage framework
is limited for real-time applications. It is also challenging to get
simultaneous ground-truth frames in the real world for training the
reconstructors. Simulators like the ESIM [17] are used to generate events from
normal videos [18], which still differ from real-world data. Due to these
challenges, end-to-end models are also widely studied for downstream tasks
(e.g., classification) [4, 2], where decisions are more interested than the
frames.

###  2.2 Point-cloud-based methods

Among the end-to-end methods, some researchers tend to consider event data as
3D point clouds in the spatial-temporal space, and directly apply point-cloud
networks. These point-cloud-based methods, such as [15, 16, 28, 11], are
originally designed to process spatial 3D point clouds with thousands of
points. The works in [16], [28] and [11] gradually aggregate local neighboring
points into centroids, hence extracting comprehensive information from the
point cloud. The schemes in [27] and [25] adopt the point cloud networks for
event data classification and demonstrate good results on the gesture
classification dataset DVS128 [1].

However, point-cloud-based methods are usually limited by the scale of input
clouds and are prone to noise. The event data contain hundreds of thousands of
events per second, which are much larger than typical point clouds studied by
the above methods (e.g., 1024 points). Therefore, the neighbor querying
algorithms in the point-cloud-based models, such as K-Nearest-Neighbours [27]
or Ball Query [11], cannot process the event clouds efficiently due to the
$O(N^{2})$ complexity. Moreover, point cloud downsampling algorithms such as
Furthest-Point-Sampling [16] are prone to outlier noise in event data, as
illustrated in Fig. 1(a). For these reasons, point-cloud-based methods are not
widely used for event-based vision.

###  2.3 Voxel-based methods

Voxel-based methods transform the event data into voxel form representations,
which are compatible with conventional CNN techniques. The representations are
either end-to-end trainable or based on heuristics. Two pioneering works are
in [21] and [9] that present the time-surface concept to convert the sparse
events into a 2D representation. Later the Event Spike Tensor (EST) [4]
assigns each event some weights, generated by a multilayer perceptron (MLP)
kernel, and projects the events to different anchor frames. The EST is a
milestone work for event-based vision and achieves SOTA performance at its
time.

Inspired by [4], Matrix-LSTM [2] proposes an encoding process with paralleled
long short-term memory (LSTM[6]) cells. The history of events that have
happened over each pixel is encoded into a vector by the LSTM, hence
generating a multi-channel image. It improves the prediction accuracy but has
a slower speed than EST.

However, for dynamic scenarios, such as human activity recognition, the above
approaches encounter motion blur problems due to the drastic compression in
the compromise of memory cost (e.g., 30k times compression to 9 frames). The
accuracy also drops because of information loss during the compression from
spatial-temporal 3D to 2D. In addition, we argue that since events are
triggered by the moving objects in the real world, the pixel-wise operation in
[4, 2] will mistakenly compress events from different objects (or parts) into
one element, causing the misalignment problem.

##  3 Methodology

In this section, we will introduce the proposed Aligned Event Tensor (AET) and
the Event Frame Net (EFN). AET is an efficient voxel-based representation for
event data, which is end-to-end trainable and alleviates problems of
information loss, motion blur and misalignment, while EFN is a versatile
network for handling both static and dynamic scenes for event-based vision.

Figure 2: Aligned Event Tensor workflow. For a better visual effect,
$\hat{M}=6$ and one Conv layer with $G=3$ are shown in this example. The
events are quantized along the time axis, and then accumulatively voxelized
(each voxel’s value is the sum of all historical events). After that, the
frames are grouped and compressed with convolutional layers.

###  3.1 Aligned Event Tensor

AET comprises three parts, namely, Timestamp Quantization, Accumulative
Voxelization and Aligned Compression, as shown in Fig. 2.

####  3.1.1 Timestamp Quantization

Event data are spatial-temporal points with an extremely high temporal
resolution, which theoretically implies a full resolution tensor of shape
$M\times H\times W$, where $M$ is the number of virtual frames, usually in the
order of millions. Due to the large size and high sparsity, it is inefficient
to process the whole tensor directly.

We also observed that accumulation of events during a proper period is needed
to form a recognizable shape. The duration can neither be too short nor too
long, otherwise the shape would be incomplete or blurred, as in Fig. 1(b).

Therefore, we propose to quantize the timestamp with a properly high temporal
resolution $\hat{M}$ (e.g., 100) to meet the computational resource
constraints. Although the optimal value for $\hat{M}$ could be sample-varying,
to facilitate batch training, a uniform time resolution is adopted for
simplicity. For each sample, the timestamps are normalized into the range of
$[0,1]$, and then evenly quantized into $\hat{M}$ bins:

$Q(\mathcal{T}_{i}^{(j)})=\max\left(\left\lceil\frac{\hat{M}\times\left(\mathcal{T}_{i}^{(j)}-\min(\mathcal{T}_{i})\right)}{\max(\mathcal{T}_{i})-\min(\mathcal{T}_{i})}\right\rceil,1\right)$
(1)

where $\mathcal{T}_{i}$ is the set of timestamps of the $i$th sample’s events
and the superscript $j$ corresponds to the $j$th event.

By Timestamp Quantization, events with adjacent timestamps are aggregated
together to form a recognizable shape. The high temporal resolution also
alleviates information loss and motion blur than low-resolution methods such
as EST [4] (9 frames). The procedure is in the left of Fig. 2.

Figure 3: Optical intensity curve, with the triggered events. The accumulated
values are shown in the bottom blue box.

####  3.1.2 Accumulative Voxelization

After quantization, the events now imply an $\hat{M}\times H\times W$ tensor,
mapped by their coordinates and quantized timestamps. Inspired by the
intrinsic mechanism of the event cameras shown in Fig. 3, we design the
Accumulative Voxelization method to determine the value of each voxel. Recall
that event cameras track optical intensity per pixel and generate a positive
or negative event whenever the change exceeds the threshold, the summation of
events at each position implies an absolute intensity change. At voxel
$(m,y,x)$, we sum the past events by their polarities:

$\mathcal{V}(m,y,x)=\sum_{j}{\mathbb{I}(x_{j}=x,y_{j}=y,\tau_{j}\leq m)p_{j}}$
(2)

where $\mathcal{V}$ is the intermediate voxel representation of the events,
$\tau_{j}=Q(\mathcal{T}_{i}^{(j)})$ is the quantized timestamp defined in
Equation (1), $\mathbb{I}$ is the indicator function and $p_{j}\in\\{-1,+1\\}$
is the polarity of the $j$th event. After the Accumulative Voxelization,
$\mathcal{V}\in\mathbb{R}^{(1\times)\hat{M}\times H\times W}$ now denotes
$\hat{M}$ frames with $C=1$ channels each.

The Accumulative Voxelization method efficiently tracks a long history of the
absolute intensity changes. Previous methods such as EST [4] and Matrix-LSTM
[2] use trainable kernels to encode histories. But they are more
computationally intensive and memory consuming, thus not feasible for high
temporal resolution as in AET, resulting in information loss. Moreover, in our
method, we record the rough intensity changes with regard to the first frame,
by summing positive and negative events together. Compared with EST alike
methods, the moving traces fades more directly. The frames in $\mathcal{V}$
will be compressed in next step.

####  3.1.3 Aligned Compression via Conv2D

During quantization and voxelization, the temporal resolution $\hat{M}$ is
kept high to alleviate information loss. However, it would be burdensome for
the neural network to process all $\hat{M}$ frames, and the contents in
adjacent frames could be similar and duplicate each other. Therefore, we
compress the voxels by merging nearby frames. For moving scenes, adjacent
frames are not completely aligned. Hence methods such as [4, 2] that conduct
pixel-wise compression suffer from blurs, lacking a mechanism for the events
to communicate in spatial.

To tackle the misalignment problem during compression, we propose using
convolutional layers to aggregate voxels within local spatial-temporal
regions. The trainable convolutional layers can project the voxel at
$(t,x+\Delta x,y+\Delta y)$ to $(t^{\prime},x,y)$ on the target frame, with a
spatial offset. Hence, the voxels corresponding to the same real-world point
can be potentially aligned to the same pixel on the anchor frame.

As shown in the right of Fig. 2, for each step, the frames are grouped with a
group size $G$, and concatenated on channel, resulting in $M^{*}=\hat{M}/G$
groups (or new frames) with $G\times C_{in}$ channels. Then they pass through
a shared Conv2d which compresses the channels to $1\times C_{out}$, meaning
that within each group, the original $G$ frames with $C_{in}$ channels each
are now compressed into $1$ frame with $C_{out}$ channels. This Conv2d
operation is equivalent to a strided Conv3d to merge frames, but is better
accelerated by PyTorch [14].

In Fig. 6, we show the comparison between our AET and EST [4]. Since this work
focuses on end-to-end training, we do not explicitly impose extra loss
function to train for sharper images. However, we observed that AET generates
sharper outputs with clearer textures. We hypothesis the spatial-temporal
compression in AET leaves the room for the model to auto-align the events, and
hence learn a sharper representation for higher classification accuracy. In
Sec. 4.5.2, the ablation studies also verify the effectiveness of AET and its
sub-modules respectively.

###  3.2 Event Frame Net

The Event Frame Net (EFN) processes the AET encoded event data for
classification tasks in three steps, CNN encoding, Frame- and Video-branch
classification, and Result synthesis. The EFN is designed to be versatile
under both static and dynamic scenes. In our ablation study, EFN is verified
to be able to boost accuracy effectively.

Figure 4: Event Frame Net workflow. The AET frames are encoded into vector
representations by the backbone CNN, and then processed in two branches.
Frame-branch makes a prediction for each frame with independent classifiers.
Video-branch treats the vectors as a sequence and process it with 1D
convolutional layers. During training, the predictions are averaged as the
final output. For inference, a synthesis module described in Fig. 5 and Sec.
3.2.4 is used to merge the results.

####  3.2.1 CNN Encoder

AET encodes event data into tensors of shape $C\times M^{*}\times H\times W$,
denoting $M^{*}$ frames with $C$ channels each. For each individual frame, we
use a shared CNN backbone model to extract geometrical information into a
vector. Following [4] and [2], we use the ImageNet-pretrained ResNet [5] with
the last decision-making layer removed.

####  3.2.2 Frame branch

As shown at the bottom of Fig. 4, the frame branch generates per-frame
predictions with an independent fully-connected layer for each individual
frame. In practice, we use a grouped $1\times 1$ convolutional layer with
$g=M^{*}$ groups to facilitate parallel computation. Also, a sliding window is
usually used for real-time event-based vision [27], so that the number of
frames, $M^{*}$, can be fixed and made compatible with our design under
different scenarios. In fact, one may also use a shared classifier for smaller
model size, but a 0.2% accuracy drop was observed.

For static scenes, the frames are similar and contain almost all the
information needed for decision making. Additionally, for dynamic tasks such
as human action recognition, one may also distinguish different actions by the
body figure from just one frame. Therefore, the frame-based classifier can
perform well under static and even some dynamic scenes.

####  3.2.3 Video branch

For certain tasks such as moving direction recognition, it is critical to make
decision based on a series of frames, rather than individuals. Therefore, we
propose the video-branch of EFN, illustrated in the upper part of Fig. 4,
where the frame representations are concatenated into a sequence and then
Conv1d modules generate the video-level prediction.

Techniques such as LSTM [6] and Transformer [23] are usually considered more
powerful for sequential data processing. However, from experiment results, we
found that Conv1d can also achieve similar accuracy with much fewer
parameters. We conjecture the sequence length $M^{*}$ is short enough so that
Conv1d can also handle well.

Figure 5: Synthesis Module. The results are merged by a weighted sum, where
weights are derived from each classifier’s performance on the validation set.

####  3.2.4 Result Synthesis

To merge the multiple predictions from video-branch and frame-branch, a
separate Synthesis module is used during inference of the EFN. The module
records per-class accuracy in the validation set, and then merges the
predictions with a weighted sum, as shown in Fig. 5.

Intuitively, since EFN is proposed to handle tasks under both static and
dynamic scenes, the reliability of classifiers in the frame-branch and video-
branch may also fluctuate correspondingly. Subsequently, we employ a
statistical method to weight these predictions. We evaluate the trained model
on the validation set while recording the per-class accuracy $Acc(p,q)$, which
is the accuracy when the classifier $p$ makes a prediction of class $q$. The
per-class accuracy is stored as a matrix shown on the left of Fig. 5.

The final output of the synthesis module is derived by a weighted sum of the
prediction logits. Firstly, each classifier generates a prediction vector,
indicating a class via ArgMax. Then the per-class accuracy is queried, with
keys being the classifier and its prediction. The per-class accuracies are
used as the weights for weighted-summing the classifiers’ predictions, which
can be formulated as:

$\textrm{out}=\sum_{p}{Acc(p,\textrm{argmax}(l_{p}))\times l_{p}}$ (3)

where $l_{p}$ is the prediction logits from classifier $p$.

##  4 Experiments

Dataset | Type | #classes | Total samples | Avg. events (kEv) | Avg. length (s) | Avg kEv/s
---|---|---|---|---|---|---
N-Caltech101 | Static | 101 | 8,709 | 112.76 | 0.30 | 375.24
N-Cars | Static | 2 | 24,029 | 3.09 | 0.10 | 30.93
N-Caltech256 | Static | 257 | 30,607 | 104.96 | 1.00 | 105.65
DVS128 | Dynamic | 11 | 1,464 | 367.98 | 6.52 | 55.62
UCF50 | Dynamic | 50 | 6,681 | 1030.83 | 6.80 | 160.25
Table 1: Statistics of the datasets used in the experiments.

Detailed experiments are conducted on our approach to evaluate the overall
performance, each module’s contribution and the efficiency. Our method
consistently surpasses all baselines under both static and dynamic scenes.

###  4.1 Datasets and Preprocessing

We evaluate our approach on datasets of static and dynamic scenes, whose
statistics are shown in Table 1. The N-Caltech101 [13] and N-Cars [21] are
benchmark datasets used in many previous studies, such as [4, 2, 21]. They
contain static daily objects and cars. Besides, recently Ref. [7] proposes a
new static dataset N-Caltech256 with many more samples and classes, and
therefore more challenging. Hence, we believe consistent improvements on this
dataset would be solid and convincing.

For dynamic scenes, the DVS128 dataset [1] records human gestures under
different illumination environments, which has been adopted in [27]. A new
dynamic and challenging dataset, the UCF-50 in event vision, is also newly
proposed by [7], which contains sports actions.

The datasets are split into training, validation, and test sets following
official practices. The DVS128 dataset is preprocessed by cutting the event
flow with a window size of 750ms and 100ms steps according to [27]. Among
them, the N-Caltech101, N-Cars, and DVS128 are recorded in the real-world,
while the other datasets are converted from conventional RGB datasets by
simulators.

###  4.2 Experimental Settings

We train our network on the five datasets with the Adam optimizer at a
learning rate of 1e-4. The learning rate is adjusted step by step with a
cosine scheduler, with a 10% linear warm-up. The model is trained for 200
epochs for N-Caltech256 and 120 epochs for other datasets.

For AET, we quantize the timestamps with $\hat{M}=100$ and use 2-step Conv2d
for the compression. Kernel size is 5, $G=10$ (2, 5 respectively) and
LeakyReLU is used for activation. $C_{in}=1,C_{out}=3$, and intermediate
channel = 4. For the EFN video-branch, 2 convolutional layers with kernel size
$k_{1}=5,k_{2}=3$ are used. We also use a max-pooling layer of size 2,
resulting in 1 element in the end. By default, Res18 is used for the backbone.

In the training process, predictions from different classifiers are averaged
and CrossEntropyLoss is used. The best model with the highest average accuracy
on the validation set is selected for testing. During inference, the Synthesis
module is used to combine the predictions.

###  4.3 Results on Static Datasets

Model | Backbone | N-Cars | N-Cal101 | N-Cal256
---|---|---|---|---
HOTS [9] | - | 62.4 | 21.0 | -
HATS [21] | Res34 | 90.9 | 69.1 | -
Two-Channel [12] | Res34 | 86.1 | 71.3 | -
Voxel Grid [29] | Res34 | 86.5 | 78.5 | -
EST [4] | Res34 | 92.5 | 83.7 | 63.11
Matrix-LSTM [2] | Res18 | 94.37 | 84.31 | -
MLSTM+E2VID[2] | Res34 | (95.65) | (85.72) | -
MLSTM+E2VID[2] | Res18 | (95.80) | (84.12) | -
AET-EFN (Ours) | Res18 | 95.91 | 89.25 | 67.86
AET-EFN (Ours) | Res34 | 95.28 | 89.95 | 68.51

Table 2: Accuracy (%) on static datasets. Second best in italic. Our method
performs strongly even with the lighter Res18. MatrixLSTM+E2VID [2] is
reconstruction-based hence not in the scope of this research, which focuses on
end-to-end training, but we still surpass their results.

As shown in Table 2, our proposed method surpasses current methods by a large
margin, including HOTS [9], HATS [21], Two-Channel Image [12], Voxel Grid
[29], Matrix-LSTM [2] and EST [4]. Our approach achieves the SOTA performance
in both static N-Caltech101 and N-Cars. Compared with the second-best
solution, AET-EFN increases the prediction accuracy by 1.54% and 5.64%,
respectively.

The N-Caltech256 is a new dataset that has not been used in previous works. We
re-implemented EST [4] on N-Caltech256 as the baseline, because EST is among
the best methods on the previous two datasets and is open-sourced. With the
new dataset, AET-EFN can still improve the accuracy by 5.4% over EST.
Considering the baselines can benefit from the large amount of data, the
improvement in N-Caltech256 strongly proves the superiority of our method.

Moreover, as in Table 4 and Sec 4.5.1, we found that more powerful backbones
can further boost the accuracy, even with fewer parameters. Specifically, AET-
EFN can achieve an impressive 91.95% accuracy on the N-Caltech101 dataset
using an EfficientNet backbone [22].

###  4.4 Results on Dynamic Datasets

Model | Type | DVS128 | UCF50
---|---|---|---
Amir et al.[1] | PointCloud | 94.4 | -
PointNet [15] | PointCloud | 88.8 | -
PointNet++ [16] | PointCloud | 95.2 | -
PAT [27] | PointCloud | 96.0 | -
EST [4] (Res18) | Voxel | 96.91 | 67.99
AET-EFN (Res18) | Voxel | 97.38 | 84.03
Table 3: Accuracy (%) on dynamic datasets.The methods listed in the table are
either point-cloud-based or voxel-based. Res18 backbone is used for both EST
and AET-EFN.

Table 3 shows the results on dynamic datasets, wherein the proposed method
outperforms all baselines and improves the accuracy to 97.38% for DVS128 and
84.03% for UCF-50. It surpasses point-cloud-based methods, such as Amir et al.
[1], PAT [27], PointNet [15] and PointNet++ [16]. While the methods in Table 2
only consider static scenes and have not been applied to dynamic datasets, we
re-implemented the EST method [4] as the representative of the voxel-based
class. Surprisingly, the EST method achieves an accuracy of 96.91%. We
hypothesize since DVS128 contains obviously different gesture classes with
little noise, it could be an easy task. Moreover, the recordings have been
clipped into short windows [27], hence the motion blur problem is also
mitigated. Nonetheless, our approach still achieves a new SOTA accuracy on
DVS128.

The UCF-50 is also a new dataset on which no existing work has been applied.
Different from DVS128, the dataset contains much longer recordings and complex
actions, hence less reasonable for clipping. The point-cloud-based methods are
not applicable due to the cloud size (1030k events per sample on average) and
noise. Therefore, we re-implemented the EST, which is also the best baseline
for DVS128, on UCF-50. For this challenging dataset, EST drops to 67.99%,
whereas our method achieves 84.03%, viz. a substantial improvement of 16.04%.

###  4.5 Ablation Study

We conduct ablation studies to analyze the contributions of different modules
in the following.

Backbone CNN | #params (Mb) | NCal101 Acc (%)
---|---|---
Res18 [5] | 11.18 | 89.25
Res34 [5] | 21.28 | 89.95
Dense121 [8] | 7.48 | 90.12
ENet_lite0 [22] | 4.03 | 90.75
ENetB2 [22] | 8.42 | 91.95
Table 4: Backbone ablation study on N-Cal101.

####  4.5.1 Backbone CNN

The EFN adopts a CNN backbone as the encoder to learn a vector representation
for each frame of AET. Previous tensor-based methods, such as [4, 2, 21],
adopt the popular ResNet [5] as their backbone. To study the effect of
different CNN models on the accuracy of EFN, we run ablation studies on the
N-Caltech101 dataset. As shown in Table 4, EFN benefits from more powerful CNN
backbones consistently. Compared with EFN + Res18 configuration, the adoption
of the EfficientNet_B2 [22] increases the accuracy by 2.7% to 91.95%. With
similar or smaller model size. Other choices, such as EfficientNet_lite and
DenseNet [8] can also improve the accuracy. These results confirm that the
improvement brought by AET-EFN is universal. Even better, the EFN framework is
orthogonal to its backbone and can harvest from more advanced CNN models.

Method | Variant | NCal101 Acc (%)
---|---|---
AET + EFN | - | 89.25
AET + Res18 | - | 86.90
EST[4] \+ EFN | - | 86.50
| Spike | 88.57
AET variants | Spike+Accum | 88.22
\+ EFN | Avg compress | 87.19
| Quantize only | 86.61
AET | no video-branch | 88.45
+EFN variants | no frame-branch | 87.24
Table 5: Ablation study on AET and EFN, the results prove that both AET and
EFN can bring significant improvements independently. For AET variations,
‘Spike’ means directly voxelizing the quantized events without Accumulative
Voxelization. ‘Avg compress’ averages adjacent frames for compression instead
of the original Aligned Compression. ‘Quantize only’ quantizes the events into
10 frames directly without the later procedures. Discussion in 4.5.2.

####  4.5.2 Ablation on AET

To quantify the contributions of the AET module and its three sub-modules, we
conduct ablation studies in Table 5. Compared with the original AET+EFN design
(89.25%), replacing the AET module with EST [4] degrades the accuracy by 2.75%
to 86.50%, showing the significance of the AET module.

Moreover, in the ‘AET variants’ part, we show the contribution of the AET sub-
modules. In this part, except the specific changes, other settings are the
same as AET. The last row ‘Quantize only’ directly quantizes events into 10
bins without the two sub-modules, Accumulative Voxelization and Aligned
Compression, resulting in a 2.64% drop. The first row ‘Spike’ omits the
Accumulative Voxelization step. Hence each voxel represents spikes of the
events rather than accumulative intensity changes. The ‘Spike+Accum’ denotes a
hybrid of the two voxelization methods. However, both of them decrease the
accuracy by 0.68% and 1.03%, respectively, showing the contribution of the
Accumulative Voxelization. Finally, the ‘Avg compress’ in the third row
replaces the Aligned Compression by averaging the frames in groups, where the
accuracy also decreases by 2.06%.

####  4.5.3 Ablation on EFN

Next, we also study the contribution of the EFN framework. As listed in Table
5, the accuracy drops by 2.35% to 86.90% without the EFN. As for the two
branches in EFN, the accuracies are 0.8% and 2.01% lower if the video-branch
and the frame-branch are removed, respectively. These results indicate that
EFN can effectively improve the accuracy and both branches contribute to the
improvement.

###  4.6 Influence of Hyper-Parameters

To study the impact of the hyper-parameters $\hat{M}$ and $G$, we conduct
experiments on N-Caltech101 as in Table 6. According to the results, the
choice of $\hat{M}$ for Timestamp Quantization and $G$ for Auto-aligned
Compression has only little impact on the final accuracy, within a proper
range. Except when the group size equals 50, the accuracy drops by 1.2%, which
is probably due to the information loss caused by the high compression rate.
We also test different kernel sizes (K) for the Auto-Aligned Compression in
AET. A larger kernel appears to benefit the prediction, which is reasonable
since it implies a larger receptive field for alignment.

$\hat{M}$ | Acc(%) |  | G | Acc(%) |  | K | Acc(%)
---|---|---|---|---|---|---|---
50 | 88.86 |  | 5 | 89.25 |  | 3 | 88.23
100 | 89.25 |  | 10 | 89.25 |  | 5 | 89.25
150 | 89.26 |  | 20 | 88.91 |  | 9 | 89.43
200 | 88.91 |  | 50 | 88.05 |  |  |
Table 6: Impact of hyper-parameters on N-Caltech101, detailed discussion in
Sec. 4.6

###  4.7 Efficiency of our approach

Method | Asynch. | Time (ms) | Speed (kEv/s)
---|---|---|---
Gabor-SNN [21] | Yes | 285.95 | 14.15
HOTS [9] | Yes | 157.57 | 25.68
HATS [21] | Yes | 7.28 | 555.74
Matrix-LSTM [2] | No | 8.25 | 468.36
EST [4] | No | 6.26 | 632.9
AET-EFN (Res34) | No | 4.64 | 933.51
AET-EFN (Res18) | No | 3.18 | 1194.20

Table 7: Inference Speed of AET-EFN. The second column indicates whether the
method processes new events asynchronously or synchronously in packages. For
fairness, we list the results for both Res18 and Res34. NCars is used
following the baseline methods.

Besides accuracy improvement, the proposed method also outperforms other
methods in terms of efficiency. Following [4, 2], we use the inference time
and the average event processing speed as metrics. We use one GTX 1080Ti GPU
and predict one sample at a time, same as in [2] for fairness. The methods are
categorized into asynchronous and synchronous classes. The asynchronous
methods, such as Gabor-SNN [21], HOTS [9] and HATS [21] process the input
events asynchronously upon their arrival. On the other hand, synchronous
methods such as Matrix-LSTM [2] and EST [4] process the events in packets.

From Table 7, the efficiency of AET-EFN is the highest among its fellows.
Following [4, 2], we conduct the efficiency test on the N-Cars dataset. To
fairly evaluate the speed of EFN, we use Res34 as the CNN backbone, same as
other baselines. Benefiting from the efficient AET representation, our
approach achieves the lowest inference latency and fastest event processing
speed, at 4.64ms and 933.51kEvents/s, respectively. Moreover, with the Res18
backbone, the AET-EFN is even faster at 3.18ms, equivalent to an impressive
314Hz frame rate.

##  5 Visualization

Figure 6: Visualizations of AET (left) and EST [4] (right) in pairs. Our
method enjoys sharper edges, clearer texture and finer details.

In Fig. 6 we visualize our AET in comparison with the EST representations.
Following [4], we take the average of the channels to reduce them to 3 for RGB
display. It is commonly observed that AET enjoys perceptibly better quality
images, which provides higher interpretability to our method (e.g., under
noisy background in the first row and clearer human faces in the third row).
Moreover, differences are even more obvious when there are texts, as shown in
the second row. Some interesting patterns, such as affluent contrast and sharp
edges have been observed, which further confirms the effectiveness of the
proposed AET. Together with the results in Sec. 4.5.2, we believe this
visualization can justify our AET as a better representation.

##  6 Conclusion

This work has proposed the Aligned Event Tensor (AET) and Event Frame Net
(EFN) for event-based vision under both static and dynamic scenes. Extensive
experiments have confirmed the integration of AET and EFN consistently
outperforms all other SOTA methods by large margins. The contributions of AET
and EFN have been extensively quantified through ablation studies. The AET
appears to be a generically better voxel-based representation for event data
under different scenes, and may be adopted for various event-vision tasks. The
versatile EFN also improves the prediction accuracy further. Compared with
existing approaches, AET-EFN is provably superior in inference time and event
processing speed, rendering it promising for real-time online tasks.

## References

  * [1] Arnon Amir, Brian Taba, David Berg, Timothy Melano, Jeffrey McKinstry, Carmelo Di Nolfo, Tapan Nayak, Alexander Andreopoulos, Guillaume Garreau, Marcela Mendoza, et al.  A low power, fully event-based gesture recognition system.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 7243–7252, 2017.
  * [2] Marco Cannici, Marco Ciccone, Andrea Romanoni, and Matteo Matteucci.  Matrix-lstm: a differentiable recurrent surface for asynchronous event-based data.  arXiv preprint arXiv:2001.03455, 2020.
  * [3] Guillermo Gallego, Tobi Delbruck, Garrick Orchard, Chiara Bartolozzi, Brian Taba, Andrea Censi, Stefan Leutenegger, Andrew Davison, Jörg Conradt, Kostas Daniilidis, et al.  Event-based vision: A survey.  arXiv preprint arXiv:1904.08405, 2019.
  * [4] Daniel Gehrig, Antonio Loquercio, Konstantinos G Derpanis, and Davide Scaramuzza.  End-to-end learning of representations for asynchronous event-based data.  In Proceedings of the IEEE International Conference on Computer Vision, pages 5633–5643, 2019.
  * [5] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 770–778, 2016.
  * [6] Sepp Hochreiter and Jürgen Schmidhuber.  Long short-term memory.  Neural computation, 9(8):1735–1780, 1997.
  * [7] Yuhuang Hu, Hongjie Liu, Michael Pfeiffer, and Tobi Delbruck.  Dvs benchmark datasets for object tracking, action recognition, and object recognition.  Frontiers in neuroscience, 10:405, 2016.
  * [8] Gao Huang, Zhuang Liu, Laurens Van Der Maaten, and Kilian Q Weinberger.  Densely connected convolutional networks.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 4700–4708, 2017.
  * [9] Xavier Lagorce, Garrick Orchard, Francesco Galluppi, Bertram E Shi, and Ryad B Benosman.  Hots: a hierarchy of event-based time-surfaces for pattern recognition.  IEEE transactions on pattern analysis and machine intelligence, 39(7):1346–1359, 2016.
  * [10] Patrick Lichtsteiner, Christoph Posch, and Tobi Delbruck.  A 128$\times$128 120 db 15$\mu$s latency asynchronous temporal contrast vision sensor.  IEEE journal of solid-state circuits, 43(2):566–576, 2008.
  * [11] Yongcheng Liu, Bin Fan, Shiming Xiang, and Chunhong Pan.  Relation-shape convolutional neural network for point cloud analysis.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 8895–8904, 2019.
  * [12] Ana I Maqueda, Antonio Loquercio, Guillermo Gallego, Narciso García, and Davide Scaramuzza.  Event-based vision meets deep learning on steering prediction for self-driving cars.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 5419–5427, 2018.
  * [13] Garrick Orchard, Ajinkya Jayawant, Gregory K Cohen, and Nitish Thakor.  Converting static image datasets to spiking neuromorphic datasets using saccades.  Frontiers in neuroscience, 9:437, 2015.
  * [14] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, et al.  Pytorch: An imperative style, high-performance deep learning library.  In Advances in neural information processing systems, pages 8026–8037, 2019.
  * [15] Charles R Qi, Hao Su, Kaichun Mo, and Leonidas J Guibas.  Pointnet: Deep learning on point sets for 3d classification and segmentation.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 652–660, 2017.
  * [16] Charles Ruizhongtai Qi, Li Yi, Hao Su, and Leonidas J Guibas.  Pointnet++: Deep hierarchical feature learning on point sets in a metric space.  In Advances in neural information processing systems, pages 5099–5108, 2017.
  * [17] Henri Rebecq, Daniel Gehrig, and Davide Scaramuzza.  Esim: an open event camera simulator.  In Conference on Robot Learning, pages 969–982, 2018.
  * [18] Henri Rebecq, René Ranftl, Vladlen Koltun, and Davide Scaramuzza.  High speed and high dynamic range video with an event camera.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 2019\.
  * [19] Olaf Ronneberger, Philipp Fischer, and Thomas Brox.  U-net: Convolutional networks for biomedical image segmentation.  In International Conference on Medical image computing and computer-assisted intervention, pages 234–241. Springer, 2015.
  * [20] Cedric Scheerlinck, Henri Rebecq, Daniel Gehrig, Nick Barnes, Robert Mahony, and Davide Scaramuzza.  Fast image reconstruction with an event camera.  In The IEEE Winter Conference on Applications of Computer Vision, pages 156–163, 2020.
  * [21] Amos Sironi, Manuele Brambilla, Nicolas Bourdis, Xavier Lagorce, and Ryad Benosman.  Hats: Histograms of averaged time surfaces for robust event-based object classification.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 1731–1740, 2018.
  * [22] Mingxing Tan and Quoc V Le.  Efficientnet: Rethinking model scaling for convolutional neural networks.  arXiv preprint arXiv:1905.11946, 2019.
  * [23] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  In Advances in neural information processing systems, pages 5998–6008, 2017.
  * [24] Lin Wang, Tae-Kyun Kim, and Kuk-Jin Yoon.  Eventsr: From asynchronous events to image reconstruction, restoration, and super-resolution via end-to-end adversarial learning.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 8315–8325, 2020.
  * [25] Qinyi Wang, Yexin Zhang, Junsong Yuan, and Yilong Lu.  Space-time event clouds for gesture recognition: from rgb cameras to event cameras.  In 2019 IEEE Winter Conference on Applications of Computer Vision (WACV), pages 1826–1835. IEEE, 2019.
  * [26] SHI Xingjian, Zhourong Chen, Hao Wang, Dit-Yan Yeung, Wai-Kin Wong, and Wang-chun Woo.  Convolutional lstm network: A machine learning approach for precipitation nowcasting.  In Advances in neural information processing systems, pages 802–810, 2015.
  * [27] Jiancheng Yang, Qiang Zhang, Bingbing Ni, Linguo Li, Jinxian Liu, Mengdie Zhou, and Qi Tian.  Modeling point clouds with self-attention and gumbel subset sampling.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 3323–3332, 2019.
  * [28] Hengshuang Zhao, Li Jiang, Chi-Wing Fu, and Jiaya Jia.  Pointweb: Enhancing local neighborhood features for point cloud processing.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 5565–5573, 2019.
  * [29] Alex Zihao Zhu, Liangzhe Yuan, Kenneth Chaney, and Kostas Daniilidis.  Unsupervised event-based learning of optical flow, depth, and egomotion.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 989–997, 2019.