Mon. Not. R. Astron. Soc. 418, 1102–1114 (2011) doi:10.1111/j.1365-2966.2011.19576.x

On the Emmenthal distribution of highly inclined asteroids

V. Carruba� and J. F. Machuca
UNESP, Univ. Estadual Paulista, Grupo de dinâmica Orbital e Planetologia, Guaratinguetá, SP, 12516-410, Brazil

Accepted 2011 August 1. Received 2011 July 11; in original form 2011 March 25

ABSTRACT
Highly inclined asteroids are objects with sin (i) > 0.3. Among highly inclined asteroids, we
can distinguish between objects with inclinations smaller than that of the centre of the ν6 =
g − g6 secular resonance and objects at higher inclinations. Using the current mechanisms of
dynamical mobility, it is not easy to increase the values of an asteroid with an initial small
inclination to values higher than that of the centre of the ν6 resonance. The presence of highly
inclined objects might therefore be related to the early phases of the Solar system.

It has been observed that several dynamically stable regions are characterized by a very low
number density of objects, unlike low-inclined bodies that tend to occupy all the dynamically
viable regions. The distribution of asteroids at a high inclination in the domain of proper
elements in dynamically stable regions resembles an Emmenthal cheese, with regions of low
number density close to highly populated areas. While this phenomenon has been observed
qualitatively in the past, no quantitative study has yet been carried out on the extent and
long-term stability of these regions.

In this paper, we identify two dynamically stable regions characterized by very low values
of number density and permanence times of 100 Myr or more when the Yarkovsky force is
considered. We show that the low number density of objects in these areas cannot be produced
as a statistical fluctuation of any simple one-dimensional statistical distribution, such as the
Poissonian, uniform and Gaussian distributions, or of a tri-dimensional distribution, such as
the tri-variate normal distribution. The presence of unoccupied dynamically stable regions
could indicate that the primordial asteroidal population might not have reached all available
zones at high-i. This sets constraints on the scenarios for the early phases of the history of our
Solar system.

Key words: celestial mechanics – minor planets, asteroids: general.

1 IN T RO D U C T I O N

Highly inclined asteroids are objects with sin (i) > 0.3. For these
bodies, the analytical theory used to obtain the proper elements is
not very accurate. Highly inclined asteroids are also characterized
by a non-uniform distribution. It has been observed that several
dynamically stable regions are characterized by a very low number
density of objects, unlike low-inclined bodies that tend to occupy
all the dynamically viable regions. While this phenomenon has
been observed qualitatively in the past (Carruba 2009b, 2010b;
Michtchenko et al. 2010), no quantitative study has yet been carried
out on the extent and long-term stability of these regions. This is
important because it might provide information on the origin of the
highly inclined population. If dynamically stable regions are not
occupied, this could indicate that the primordial population might
not have reached all available zones at high-i. This sets constraints

�E-mail: vcarruba@feg.unesp.br

on the scenarios for the early phases of the history of our Solar
system. The void zones in the asteroid distribution (in our analogy,
the holes in the Emmenthal cheese), as well as the regions of higher
asteroid number density, might therefore still carry the memory of
events associated with the primordial phase of planetary migration.

In this paper, we investigate the long-term stability of underpop-
ulated areas of highly inclined asteroids by using synthetic proper-
element maps and simulations of highly inclined real and fictitious
objects, which include the Yarkovsky force. We then determine the
probability that an originally uniform distribution of asteroids, a
population following Poissonian or uni- and tri-variate Gaussian
statistical distributions, describes the currently observed asteroid
distribution.

This paper is arranged as follows. Having introduced the problem
of stable underpopulated and overpopulated areas among highly in-
clined objects in the introduction, we revise the current knowledge
on asteroid synthetic proper elements for highly inclined asteroids
in Section 2. In Section 3, we discuss the orbital distribution of large
asteroids and in Section 4 we consider statistically the probability

C© 2011 The Authors
Monthly Notices of the Royal Astronomical Society C© 2011 RAS

 at FundaÃ
§Ã

£o C
oordenaÃ

§Ã
£o de A

perfeiÃ
§oam

ento de Pessoal de N
Ã

­vel Superior on July 15, 2013
http://m

nras.oxfordjournals.org/
D

ow
nloaded from

 

http://mnras.oxfordjournals.org/


Distribution of highly inclined asteroids 1103

Figure 1. Panel A shows the (a, e) projection of highly inclined AstDyS asteroids. The blue asterisks display asteroids with a standard deviation on σ a between
0.0003 and 0.01, while red circles show asteroids with σ a larger than 0.01. The vertical red lines display the locations of mean-motion resonances, while the
blue lines show the positions of linear secular resonances. Panel B shows the [a, sin (i)] projection of the same asteroids.

that various asteroid distributions might replicate the observed dis-
tribution. In Section 5, we obtain the number density distributions
(density maps) of such objects for large (H < 12) and small objects,
and we identify potential dynamically stable regions with a low
number density of objects. In Section 6, we obtain synthetic proper-
element maps for the underpopulated areas identified in Section 3.
In Section 7, we study the long-term stability of these regions when
the Yarkovsky force is considered. Finally, in Section 8, we present
our conclusions.

2 SYNTHETIC PROPER ELEMENTS FOR
H I G H LY IN C L I N E D O B J E C T S

The first step to identifying underpopulated dynamically stable re-
gions is to have a reliable and updated set of asteroid proper ele-
ments. We start by revising the data available in the literature.

On 2010 December 8, the AstDyS site1 listed a set of 240 831
asteroids for which the synthetic proper elements (and their errors)
– obtained using the approach of Knežević & Milani (2003) – are
available. We selected 10 073 objects of high inclination character-
ized by a semimajor axis in the range of 2.0–3.4 au and by sin (i) >

0.3. We start by looking at projections in the (a, e) and [a, sin (i)]
planes.

Fig. 1 displays (a, e) (panel A) and [a, sin (i)] (panel B) projec-
tions of the AstDyS proper elements for the highly inclined objects.
Following Gil-Hutton (2006), we can divide the highly inclined pop-
ulation of asteroids into three regions: the Phocaea region (zone A
in Gil-Hutton 2006) between the 7J:-2A and 3J:-1A mean-motion
resonances, the Hansa and Pallas region (zone B) between the 3J:-
1A and the 5J:-2A mean-motion resonances and the Euphrosyne
region (zone C) between the 5J:-2A and the 2J:-1A mean-motion
resonances. The thick vertical lines are associated with the main
mean-motion resonances with Jupiter, while the thin lines show the
location of higher-order mean-motion resonances. The blue lines
show the location of the centre of the three main linear secular reso-
nances (ν6 = g − g6, ν5 = g − g5 and ν16 = s − s6), computed using

1 Knežević & Milani (2003) http://hamilton.dm.unipi.it/astdys

the second-order and fourth-degree secular perturbation theory of
Milani & Knežević (1992, 1994). Secular resonances occur when
there is a commensurability between the precession frequency of
an asteroid’s pericentre g or node s and that of a planet gi, si, where
the suffix i = 1, . . . , 8 indicates the planet number as a function of
the distance from the Sun (1 for Mercury, 2 for Venus, etc.). We
warn the reader that the secular perturbation theory of Milani &
Knežević (1992, 1994), which is based on an expansion in series
of eccentricities and inclination of the perturbing function, loses
accuracy at high inclination and close to mean-motion resonances.
Thus, in Fig. 1 and in the subsequent figures, we report the location
of the secular resonances given by this theory just as a qualitative
indication, but we do not expect it to match the distribution of actual
asteroids. The blue asterisks display asteroids with σ a (the standard
deviation of the values of the semimajor axis as computed with
the approach of Knežević & Milani 2003) between 0.0003 au (the
limit given by Knežević & Milani 2003 for ‘stable’ synthetic proper
elements) and 0.01 au (the limit for pathological cases). The red
circles show asteroids with σ a larger than 0.01 au.2 We remind the
reader that objects at very high eccentricities (larger than about 0.4)
are subject to close encounters with Mars and other planets, and are
unstable on time-scales of a few Myr.

As can be seen in the figure, we can observe a region of low
asteroid density near the ν6 resonance border, which is caused by
an increase in the eccentricity that asteroids on circulating and
librating orbits experience because of the resonance topology (for
a more in-depth discussion of the mechanism of asteroid depletion,
see also Carruba 2010a; Carruba, Machuca & Gasparino 2011b). An
exception to this mechanism is observed for asteroids in antialigned
librating states of the ν6 resonance, such as those of the Tina family
(Carruba & Morbidelli 2011a), whose resonant states shield them
from experiencing planetary close encounters. Besides the regions
influenced by the ν6 resonance, we can see a region of low asteroid

2 Knežević & Milani define stable, unstable and pathological proper ele-
ments as elements having ‘small’, ‘medium’ and ‘large’ relative errors. The
thresholds for the three cases with respect to each element are discussed in
the text.

C© 2011 The Authors, MNRAS 418, 1102–1114
Monthly Notices of the Royal Astronomical Society C© 2011 RAS

 at FundaÃ
§Ã

£o C
oordenaÃ

§Ã
£o de A

perfeiÃ
§oam

ento de Pessoal de N
Ã

­vel Superior on July 15, 2013
http://m

nras.oxfordjournals.org/
D

ow
nloaded from

 

http://mnras.oxfordjournals.org/


1104 V. Carruba and J. F. Machuca

Figure 2. Panel A shows the [e, sin (i)] projection of highly inclined AstDyS asteroids. The blue asterisks display asteroids with a standard deviation on σ e

between 0.003 and 0.1, while red circle show asteroids with σ e larger than 0.1. The blue lines show the positions of linear secular resonances. Panel B shows
the [e, sin (i)] projection of the same asteroids, but this time the blue asterisks display asteroids with a standard deviation on σ i between 0.001 and 0.03, while
the red circles show asteroids with σ i larger than 0.03.

Figure 3. Panel A shows the (a, g) projection of highly inclined AstDyS asteroids. The blue asterisks display asteroids with a standard deviation on σ g between
1 and 10, while the red circles show asteroids with σ g larger than 10. Panel B shows the (a, s) projection of the same asteroids, but this time the symbols refer
to σ s (the range of values for the errors is the same).

density between the ν5 and ν16 resonances in the [a, sin (i)] plane
in the region of the Pallas family (see Fig. 1, panel B, and Carruba
2010b). More strikingly, there is a low-density region of asteroids
between the 5J:-2A and 7J:-3A mean-motion resonances, which
seems otherwise to be relatively stable. We investigate these subjects
further in Sections 3 and 6.

In Fig. 2, we display [e, sin (i)] projections of asteroids with
values of the errors in e (panel A) and i (panel B) for ‘stable’ (black
dots), ‘unstable’ (blue circles) and ‘pathological’ (red circles) proper
elements, according to Knežević & Milani (2003). As in Fig. 1, the
blue lines display the locations of the main secular resonances in
the region. As can be seen, the errors in eccentricity and inclination
are larger for asteroids at high inclination, with very few objects
having large errors at low eccentricities.

Finally, Fig. 3 displays (a, g) (panel A) and (a, s) (panel B)
projections of highly inclined AstDyS asteroids. We can see the

departure from a linear dependence for g values as a function of a
near the 2J:-1A and 3J:-1A mean-motion resonances, as discussed
by Carruba & Michtchenko (2009). The negative values of g near the
Hansa family are related to objects with small proper eccentricities
(e < 0.0179), and to the problem of determining the correct proper
frequency g for asteroids whose elements (k, h) pass through zero
(Carruba 2010b). The ‘inclined’ dependence of s values as a function
of a is discussed by Carruba & Michtchenko (2007).

3 A STEROI D O RBI TAL DI STRI BU TI ON

When studying the current orbital distribution of objects, it is nec-
essary to distinguish between ‘large’ and ‘small’ asteroids. Large
objects have experienced very little orbital mobility, which is caused
by the Yarkovsky effect, and it can be assumed that their current po-
sition is not far from the position they reached during the last stages

C© 2011 The Authors, MNRAS 418, 1102–1114
Monthly Notices of the Royal Astronomical Society C© 2011 RAS

 at FundaÃ
§Ã

£o C
oordenaÃ

§Ã
£o de A

perfeiÃ
§oam

ento de Pessoal de N
Ã

­vel Superior on July 15, 2013
http://m

nras.oxfordjournals.org/
D

ow
nloaded from

 

http://mnras.oxfordjournals.org/


Distribution of highly inclined asteroids 1105

Figure 4. An [a, sin i(i)] plot of low-eccentricity (e < 0.175) H < 12 asteroids (panel A) and high-eccentricity (e > 0.175) H < 12 asteroids (panel B). The
red asterisks identify regions affected by the ν6 resonance, according to the model of Yoshikawa (1987). The dashed vertical lines display regions partially
depopulated of asteroids near the mean-motion resonances.

of planetary migration. However, small objects tend to be displaced
by the Yarkovsky force and, as a result, by mean-motion and secular
resonances. Their current distribution can provide hints about the
dynamical stability of the regions: underpopulated areas might be
associated with dynamical unstable regions (but not necessarily, as
we see later in this paper).

We define a ‘large’ object as an asteroid with an absolute magni-
tude H > 12. We choose this criterion because asteroids with this ab-
solute magnitude are not significantly affected by non-gravitational
effects. If we estimate the asteroid diameter using the relationship
(Carruba et al. 2003)

D = D0√
pV

× 10−0.2H , (1)

where D0 = 1329 km, pV is the geometric albedo and H is the
asteroid’s absolute magnitude. Even for very large values of pV =
0.5, an asteroid with H = 12 will have a diameter of 7.5 km. At this
size, the YORP effect is negligible. Using an obliquity of the spin
axis of ±90◦ and parameters typical of C-type asteroids (Carruba
et al. 2003), the mobility caused by the Yarkovsky force is at most
0.1635 au over 4.5 Gyr. Therefore, we should not expect the current
orbital locations of such bodies to have changed much since the
formation of the Solar system.3

To investigate the locations of stable areas characterized by a low
population of objects, we have obtained the distribution of the num-
ber density of asteroids in the [a, sin (i)] representative plane. We
concentrate on the [a, sin (i)] plane because we are interested in ar-
eas that were originally depleted in asteroids, and eccentricities are
more easily changed by mechanisms of dynamical mobility than by
proper inclinations (semimajor axes are affected by the Yarkovsky
force). Regions in the [a, sin (i)] plane with inclinations different

3 Of course, the orbital mobility of an object depends on several other param-
eters than the absolute magnitude, such as the spin axis obliquity, the thermal
inertia, the thermal conductivity, the surface density, the bond albedo, etc.
Thus, this criterion is approximated at best. However, because many of these
parameters are not known for large numbers of asteroids, it is helpful to have
a simple criterion that can be used for statistical considerations.

from those of the bodies originally injected in highly inclined or-
bits can therefore still be characterized by a relatively low number
density of objects.

To identify regions that were originally underpopulated, we need
to eliminate areas whose low number density is caused by the desta-
bilizing effects of the local dynamics. Based on several previous
studies of the dynamics of highly inclined objects (see, for instance,
Carruba 2009a,b, 2010a,b; Carruba et al. 2011b), we can distinguish
between destabilizing resonances (resonances that cause a change
in the asteroid eccentricity that is large enough to bring its pericentre
to regions of terrestrial planet encounters, causing the loss of at least
50 per cent of the bodies initially inside the resonance because of a
collision with the Sun or other planets during time-scales of 10 Myr
or less, in conservative simulations of the orbital evolution of an
asteroid subjected to the influence of the Sun and all the planets)
and diffusive resonances (resonances whose passage through causes
a significant change in the asteroid proper elements, but not large
enough to cause the loss of the body). We can safely assume that
bodies in regions affected by or near destabilizing resonances will
be lost on time-scales of up to 10 Myr, so that local low densities
do not reflect a non-uniform primordial distribution. The destabiliz-
ing effect of the resonances discussed are also investigated in some
detail later in this paper.

Among destabilizing mean-motion resonances, we can identify
from previous works the 7J:-2A, 3J:-1A, 5J:-2A and 2J:-1A mean-
motion resonances with Jupiter. Fig. 4 shows [a, sin (i)] projections
of asteroids with absolute magnitude H < 12 in the highly inclined
area. Following the example of Michtchenko et al. (2010), we di-
vided the population of objects into low-eccentricity (e < 0.175;
Fig. 4, panel A) and high-eccentricity (e > 0.175; Fig. 4, panel B)
asteroids. Both panels show vertical band (red dashed lines) that
are related to the boundaries of the unstable regions associated with
the destabilizing mean-motion resonances. To estimate the sizes
of these boundaries, we used the following procedure. First, we
computed the size of the resonance in the (a, e) plane using the
approach described in Morbidelli (2002). Then, we added 0.015 au,
the size of the depopulated region near the borders of the 3J:-1A
mean-motion resonance, as observed by Guillens, Vieira Martins
& Gomes (2002), to the value of the size in the semimajor axis

C© 2011 The Authors, MNRAS 418, 1102–1114
Monthly Notices of the Royal Astronomical Society C© 2011 RAS

 at FundaÃ
§Ã

£o C
oordenaÃ

§Ã
£o de A

perfeiÃ
§oam

ento de Pessoal de N
Ã

­vel Superior on July 15, 2013
http://m

nras.oxfordjournals.org/
D

ow
nloaded from

 

http://mnras.oxfordjournals.org/


1106 V. Carruba and J. F. Machuca

Table 1. Cell identification, the limits in a, e and sin (i), and the number of real H < 12 objects.

Cell identification Intervals in a Intervals in e Intervals in sin (i) No. of real H < 12 objects

1 a < a31 e < 0.175 0.3 < sin i < 0.45 1
2 a < a31 e < 0.175 sin i > 0.45 0
3 a < a31 0.175 < e < 0.6 0.3 < sin i < 0.45 15
4 a < a31 0.175 < e < 0.6 sin i > 0.45 2
5 a31 < a < a83 e < 0.175 0.3 < sin i < 0.45 8
6 a31 < a < a83 e < 0.175 sin i > 0.45 0
7 a31 < a < a83 0.175 < e < 0.6 0.3 < sin i < 0.45 7
8 a31 < a < a83 0.175 < e < 0.6 sin i > 0.45 4
9 a83 < a < a52 e < 0.175 0.3 < sin i < 0.45 8

10 a83 < a < a52 e < 0.175 sin i > 0.45 1
11 a83 < a < a52 0.175 < e < 0.6 0.3 < sin i < 0.45 3
12 a83 < a < a52 0.175 < e < 0.6 sin i > 0.45 5
13 a52 < a < a94 e < 0.175 0.3 < sin i < 0.45 13
14 a52 < a < a94 e < 0.175 sin i > 0.45 1
15 a52 < a < a94 0.175 < e < 0.6 0.3 < sin i < 0.45 15
16 a52 < a < a94 0.175 < e < 0.6 sin i > 0.45 1
17 a94 < a < a21 e < 0.175 0.3 < sin i < 0.45 115
18 a94 < a < a21 e < 0.175 sin i > 0.45 5
19 a94 < a < a21 0.175 < e < 0.6 0.3 < sin i < 0.45 47
20 a94 < a < a21 0.175 < e < 0.6 sin i > 0.45 2

calculated at e = 0.25, the mean value of the eccentricity of aster-
oids in the region. Asteroids in these areas can be destabilized by
mean-motion resonances either directly or indirectly via migration
caused by the Yarkovsky force.

Concerning secular resonances, we considered the effect of the ν6

resonance, and we computed the unstable regions using the criterion
described in Carruba (2010a) and Carruba et al. (2011b). We first
computed the k6 parameter, defined as

k6 = b − ν6, (2)

where ν6 = 28.243 arcsec yr−1 is the precession frequency of Sat-
urn’s pericentre and b is a frequency-like parameter computed using
the Yoshikawa (1987) model of the ν6 resonance. Asteroids with
k6 values in the interval to −2.55 < k6 < 2.55 arcsec yr−1 in
the internal main belt (values of a smaller than the centre of the
3J:-1A mean-motion resonance; Carruba 2010a) and in the inter-
val −3.5400 < k6 < 1.6815 arcsec yr−1 in the central main belt
(between the 3J:-1A and 5J:-2A mean-motion resonances; Carruba
et al. 2011b) are considered likely to have their eccentricity pushed
to terrestrial planet-crossing levels and to be destabilized on time-
scales of 10 Myr.4 In the external main belt, it is not possible to
use the Yoshikawa approach because of the perturbing effect of
the 2J:-1A mean-motion resonance, not accounted for in the model.
However, the effect of the ν6 resonance is more limited in this region
of the asteroid belt than in the internal part (see Section 6).

In order to have a sample of objects as close as possible to the
primordial population, we have eliminated from the list of 313
highly inclined H < 12 asteroids those objects that are members of
the following families: Phocaea, Barcelona, Hansa, Gersuind, Pal-
las, Gallia, Tina, Euphrosyne, Alauda, Luthera and (16708) (1995
SP1); see Carruba (2009b, 2010b) for the determination of the fam-
ily membership of highly inclined objects. This left us with a list

4 The exception is the Tina family members in the central main belt whose
orbits have k6 values in the range discussed in Carruba et al. (2011b),
but whose eccentricities do not reach planet-crossing levels because of the
protecting effect of the antialigned ν6 resonant configuration, as discussed
in Carruba & Morbidelli (2011a).

of 253 asteroids whose current orbits should not have been greatly
modified since the formation of the asteroid belt.

Then, we divided the highly inclined asteroid regions into 20 areas
according to these criteria. We have used the main mean-motion
resonances with Jupiter as ‘natural dynamical boundaries’ and we
have created five regions: a < a31, a31 < a < a83, a83 < a < a52, a52 <

a < a94 and a94 < a < a21.5 This division creates slightly unequal
cells in the semimajor axis, but we believe that this inconvenience
is more than compensated by the fact that such a division respects
the local dynamics and does not tend to incorporate in the same cell
objects belonging to different dynamical groups. We then divided
the asteroids into two large eccentricity cells, one for objects with
e < 0.175 and the other for objects with larger eccentricity. We
assume that, as observed, no asteroid has an eccentricity larger than
0.6, so that the limits of the cells at high eccentricity are 0.175 < e <

0.6. Finally, we also created two cells in inclination separated by the
value sin (i) = 0.45, which roughly corresponds to the inclination
of the centre of the ν5 resonance in the central main belt. Overall,
we created 20 cells in the [a, e, sin (i)] domain. Fig. 4 shows [a,
sin (i)] projections of low-eccentricity (e < 0.175; panel A) and
high eccentricity (e > 0.175; panel B) objects. Real asteroids in the
region, selected according to the criterion described earlier in this
section, are identified as blue asterisks, and the numbers in each cell
are related to the number of objects per cell. Table 1 displays the
cell numbers (from 1 to 20), the intervals in a, e and sin i of each
cell and the number of real objects in each cell.

From this preliminary analysis, we can see that the number of
H < 12 objects above the ν6 resonance is indeed limited, especially
at low eccentricities. We have identified 23 objects at low eccentric-
ities (e < 0.175) and 59 objects at high eccentricities [e > 0.175;
(2) Pallas is one of these objects]. These are characterized as having
values of proper g, the frequency of precession of the pericentre,
lower than 28.243 arcsec yr−1, a value that corresponds to the ν6 res-
onance. The limited number of large objects above the ν6 resonance

5 For the cell boundaries, we have used the centre of the resonance, and the
suffixes (e.g. 21) refer to the mean-motion resonances.

C© 2011 The Authors, MNRAS 418, 1102–1114
Monthly Notices of the Royal Astronomical Society C© 2011 RAS

 at FundaÃ
§Ã

£o C
oordenaÃ

§Ã
£o de A

perfeiÃ
§oam

ento de Pessoal de N
Ã

­vel Superior on July 15, 2013
http://m

nras.oxfordjournals.org/
D

ow
nloaded from

 

http://mnras.oxfordjournals.org/


Distribution of highly inclined asteroids 1107

when compared with the analogous population at lower inclinations
could suggest that the original population of highly inclined objects
was limited to begin with.

4 STATISTICAL CONSIDERATIONS

In the previous sections, we have seen that there are large variations
in the orbital distribution of asteroids in highly inclined regions. In
this section, we try to understand whether this phenomenon is com-
patible with normal statistical fluctuations of a uniformly spaced or
Gaussian-spaced population or whether we need some further hy-
pothesis on the mechanisms that formed the original highly inclined
populations. In our statistical analysis, we use two approaches. First,
we consider the one-dimensional probability distributions for the
H < 12 asteroid population using the cells defined in Section 3,
checking for possible correlations among cells. Secondly, we try
to see if the [a, e, sin (i)] tri-dimensional distribution follows a tri-
variate Gaussian distribution by performing a normality test. We
start by looking at the one-dimensional probability distributions.

4.1 One-dimensional probability distributions

When looking at one-dimensional probability distributions, there
are two possible approaches. We can look at the H < 12 asteroid
population and verify whether it follows simple one-dimensional
probability distributions, such as the Poissonian, uniform or Gaus-
sian distributions, in each cell. We can also check if there are cor-
relations among cells at same inclinations or in the same range of
semimajor axes. We start our analysis by looking at the Poissonian
statistical distribution.

Table 2 displays the cell identification, as given in Table 1, the
number of real objects in each cell, the weighted mean number of
objects per cell, the probability that this number is generated by
Poissonian statistics, the number of objects that would be produced
in each cell if we assume a uniform distribution of objects in [a,

e, sin (i)], which best fit the data, and the number of objects that
would be produced if we assume a Gaussian distribution in [a, e,
sin (i)], which best fit the data. We have found that the H < 12
asteroid population has mean values of a, e and sin (i) of 3.0270 au,
0.1516 and 0.3656, respectively. The variances in a, e and sin (i) are
0.0653 au, 0.0057 and 0.0030, respectively. We have estimated the
probability that a number of objects are produced by the Poissonian
distribution, assuming that the expected number of k occurrences in
a given interval is given by

f (k, λ) = λke−λ

k!
. (3)

Here, λ is a positive real number, equal to the expected number of
occurrences in the given interval. For our purposes, we used for
λ the mean values of objects per cell of real asteroids, weighted
to account for the different volume occupied by each cell. This is
given by

λ = Nast
Vcell

VTot
, (4)

where Nast is the total number of objects with H < 12 in stable
regions (253), Vcell is the volume in (a, e, i) space occupied by each
cell – see Section 3 for the limits of the cell in the (a, e, i) domain –
and VTot is the total volume in the region. Using standard statistical
terminology, we define the null hypothesis as the possibility that
the data are drawn from a given distribution. We can reject the null
hypothesis if it is associated with a probability lower than a threshold
that is usually of the order of 1 per cent or 0.01. Our results show
that 13 cells have negligible probability that their numbers can be
reproduced by a Poissonian distribution. The cells compatible with
a Poissonian distribution (cells 5, 7, 8, 9, 11, 12 and 13) are those
associated with the Hansa, Gallia and Brucato, Gersuind, Pallas and
Danae regions. Other cells have an orbital distribution of highly
inclined objects incompatible with a Poissonian distribution, either
because they are underpopulated (cells 1, 2, 4, 6, 10, 14, 16, 18 and
20) or overpopulated (cells 3, 15, 17 and 19). Because 14 out of

Table 2. Number of real objects, the weighted mean number of objects per cell, the probability of real objects
being a product of a Poissonian statistical distribution, the number of simulated objects with uniform and Gaussian
distributions as a function of the cell in the [a, e, sin (i)] domain.

Cell identification No. of Weighted mean number of Poissonian Uniform Gaussian
real objects objects per cell probability distribution distribution

1 1 11.598 1.1 × 10−4 15 2
2 0 18.761 7.1 × 10−9 15 8
3 15 7.543 0.0058 41 51
4 2 18.320 1.8 × 10−6 23 39
5 8 6.810 0.1264 4 1
6 0 16.548 6.5 × 10−8 4 2
7 7 4.431 0.0792 18 31
8 4 10.760 0.0119 4 27
9 8 3.274 0.0123 7 2
10 1 7.950 0.0028 6 2
11 3 2.128 0.1913 8 11
12 5 5.170 0.1750 7 3
13 13 8.520 0.0400 6 5
14 1 20.692 2.1 × 10−8 3 3
15 15 5.540 0.0004 19 10
16 1 13.457 1.9 × 10−5 14 13
17 115 8.394 <10 × 10−11 12 4
18 5 20.385 4.1 × 10−5 6 3
19 47 5.459 <10 × 10−11 20 14
20 2 13.257 1.5 × 10−4 21 22

C© 2011 The Authors, MNRAS 418, 1102–1114
Monthly Notices of the Royal Astronomical Society C© 2011 RAS

 at FundaÃ
§Ã

£o C
oordenaÃ

§Ã
£o de A

perfeiÃ
§oam

ento de Pessoal de N
Ã

­vel Superior on July 15, 2013
http://m

nras.oxfordjournals.org/
D

ow
nloaded from

 

http://mnras.oxfordjournals.org/


1108 V. Carruba and J. F. Machuca

20 of the considered cells do not satisfy the null hypothesis for the
Poissonian distribution, it seems unlikely that the observed changes
in orbital distribution can be explained as simple fluctuations of a
Poissonian distribution for the whole region.

We have checked whether the asteroid number distribution in
areas in the same ranges of semimajor axes, eccentricities and in-
clinations can be approximated by a Poissonian distribution. We
started by looking at asteroids in the same range of semimajor axes,
with a < a31 (Phocaea region), a31 < a < a83 (Hansa region), a83 <

a < a52 (Pallas region), a52 < a < a94 and a94 < a < a21 (Eu-
phrosyne region). Using equation (3), we computed the probability
that the asteroid proper inclination followed the Poissonian distri-
bution. Because the Poissonian distribution is only computed for
integer values of k, we multiply values of sin (i) by a factor of 100
and take only the first two digits. We best-fitted the value of λ in
equation (3) for the observed asteroids using the MATLAB routine
poissfit. We generated a random population of asteroids following
the Poissonian distribution using the MATLAB routine poissrnd. To
check if the simulated distribution was compatible with the ob-
served distribution, we performed a Pearson χ 2 test by computing
a χ 2-like variable, defined as

χ 2 =
nint∑
i=1

(qi − pi)2

qi

. (5)

Here, nint is the number of intervals of our variables, corresponding
to the number of asteroids in each region, qi is the number of real
objects in the ith interval and pi is the number of simulated objects
in the same interval. Values of χ 2 of the order of the number of
intervals are consistent with a good fit. For each region of inter-
est, we obtained 1000 random populations, and we computed the
value of χ 2 for each simulated population. Then, using the standard
χ 2 probability distribution (Press et al. 2001), we calculated the
probability that the real population might be following a Poissonian
distribution. The minimum value of χ 2 among the 1000 simulated
populations gave the best likelihood estimation that such hypotheses
were true.

The lowest χ 2 value (52.98) was observed for a set of 18 simulated
bodies in the Phocaea region. This corresponds to a probability that
the real distribution is Poissonian of just 8.2 × 10−9. All the other
regions have higher values of χ 2 and negligible probabilities that the
observed orbital distributions are Poissonian. We also checked the
distributions in the semimajor axes for low-inclined [sin (i) < 0.45]
and high-inclined [sin (i) > 0.45] asteroids, and for low-eccentricity
(e < 0.175) and high-eccentricity (0.175 < e < 0.600) objects. Once
again, we did not find values of χ 2 compatible with a Poissonian
distribution. With the possible exception of the Phocaea region,
large asteroids in the area are not distributed according to Poissonian
statistics.

We then tried a different approach by hypothesizing that fluctua-
tions in a uniform distribution could explain the observed asteroid
distribution. Using a random number generator (Press et al. 2001),
we created a set of 253 fictitious bodies with a uniform distribution
in the [a, e, sin (i)] domain. (If a body was in one of the unstable re-
gions discussed in Section 3, we eliminated it from our database and
created a new one. The procedure was repeated until we obtained
a set of 253 bodies.) Using different seeds6 (from a value of −1

6 Random number generators are algorithms that can automatically create
long runs of numbers with good random properties. The string of values
generated by such algorithms is generally determined by a fixed number,
called a seed.

up to −1000, the seed must be a negative number) for the random
number generator, we produced 1000 sets of 253 fictitious aster-
oids with a uniform distribution in the given range of [a, e, sin (i)]
values. The fifth column of Table 2 displays the number of objects
per cell produced with a uniform distribution in the [a, e, sin (i)]
domain that best fitted the results. To check if the observed distri-
bution is compatible with the uniform distribution, we performed
a one-dimensional Kolmogorov–Smirnov (KS) test. This is a test
used to verify the validity of the hypothesis that a random deviate
follows a continuous probability distribution. Given two cumulative
distributions, SN1 (x) and SN2 (x), the KS statistics D is given by the
maximum value of their absolute difference. For given values of D,
we can obtain the probability of the null hypothesis that the two
data sets are drawn from the same distribution (Press et al. 2001).
For the uniform distribution, not even considering the cells with
a94 < a < a21, whose large asteroidal population could be caused
by the larger inclination values of the ν6 resonance in the region
(more asteroids could come to these areas from low-inclinations re-
gions compared to other cells), the KS test yields for the best-fitting
case the low probability of 4.3 × 10−4 that the two data sets are
drawn from the same distribution. This is enough to reject the null
hypothesis.

We have also checked if the uniform distribution can approxi-
mate the number distribution of asteroids in the same range of a, e,
i, as done for the Poissonian distribution. We generated 1000 sets
of uniform deviates in inclination for the observed range of values
of real asteroids in each of the regions in a used for the Poissonian
analysis, and for low- and high-eccentricity and inclined objects,
and we performed a KS test for each of the samples compared with
the observed population. The highest probability results of the KS
test were once again obtained for the Phocaea region, with a prob-
ability that the distribution follows a uniform deviate of 4.26 ×
10−3. The low value of this probability seems to infer that it is un-
likely that asteroids in the Phocaea region follow a uniform deviate.
Surprisingly, we obtained a probability of 1.5 × 10−3 that highly
inclined asteroids are distributed according to uniform deviates. All
other KS tests yielded lower probability results.

Finally, we also created 1000 sets of 253 synthetic objects that
have a Gaussian distribution in (a, e, i), using the same criteria used
for the bodies created with the uniform distribution. Mean values
and variances were computed from the observed sample of large
asteroids. We did not account for possible correlations among the
variables in this preliminary test. We discuss this point in more
detail in Section 4.2. The sixth column in Table 2 displays the
number of objects per cell produced with such a distribution. We
then implemented our KS test for cells up to a = a94 and for all
cells, and we obtained best probability values of the two data set
drawn from the same distribution of 4.3 × 10−4 and 4.1 × 10−4, low
enough to be able to safely reject the null hypothesis. In either case,
it is unlikely that the current asteroid distribution is compatible with
the simulated distribution.

Once again, we checked if the Gaussian distribution can approx-
imate the number distribution of asteroids in the same range of a,
e and i. We repeated the analysis as for the uniform distribution
and we also used the MATLAB built-in routine kstest.m to check for
normality. The results are in agreement with those found previously
for the uniform distribution – there is a low probability (4.5 × 10−3)
that asteroids in the Phocaea region might follow a Gaussian dis-
tribution and a relatively high (when compared with other results)
likelihood (1.7 × 10−3) that highly inclined asteroids might be dis-
tributed normally. All other KS tests gave lower probabilities that
the two distributions, real and fictitious, could both be Gaussian.

C© 2011 The Authors, MNRAS 418, 1102–1114
Monthly Notices of the Royal Astronomical Society C© 2011 RAS

 at FundaÃ
§Ã

£o C
oordenaÃ

§Ã
£o de A

perfeiÃ
§oam

ento de Pessoal de N
Ã

­vel Superior on July 15, 2013
http://m

nras.oxfordjournals.org/
D

ow
nloaded from

 

http://mnras.oxfordjournals.org/


Distribution of highly inclined asteroids 1109

Overall, it seems that current non-uniformities in the number dis-
tributions of asteroids cannot be explained as fluctuations of one-
dimensional statistical distributions, such as the Poissonian, uniform
and Gaussian distributions, except possibly for the Phocaea region,
which is marginally compatible with a Poissonian distribution of the
inclinations. In the following subsection, we investigate if a multi-
variate Gaussian distribution could better approximate our sample
of large asteroids.

4.2 Tri-variate Gaussian probability distributions

In Section 4.1, we checked if one-dimensional probability distri-
butions can account for the observed fluctuations in the number of
large (H < 12) asteroids. Here, we investigate if a multidimensional
distribution in the [a, e, sin (i)] domain, which accounts for possible
correlations among the different dimensions, could better explain
the observed distribution.7

We assume that the distribution in the [a, e, sin (i)] domain has a
tri-variate Gaussian probability density function, given by

f (x) = 1

(2π)3/2
√|�| exp

[
−1

2
(x − μ)′�−1(x − μ)

]
. (6)

Here, x is the vector [a, e, sin (i)] and μ is the vector with com-
ponents equal to the mean values of a, e and sin i of the observed
population (see Section 4.1). � is the covariance matrix, which can
be estimated numerically having n xj vectors using

� = 1

n

n∑
j=1

(xj − μ) · (xj − μ)T. (7)

The numerical elements of � for the 253 highly inclined asteroids
in the region are equal to

� =

∣∣∣∣∣∣∣∣
0.0653 −0.0066 −0.0031

−0.0066 0.0057 0.0007

−0.0031 0.0007 0.0030

∣∣∣∣∣∣∣∣
. (8)

|�| = 9.2913 × 10−7 au is the determinant of the covariance matrix.
We first draw 253 random vectors from this tri-variate Gaussian
distribution using the following procedure. Because � is positive
defined, we use the Choleski decomposition to obtain a matrix A
such that A · AT = �. The A matrix is a triangular matrix whose
elements in the triangle below the main diagonal are all zero. We
then generated 253 tri-dimensional vectors Z, whose components
are N independent standard normal variates, using standard MATLAB

routines, such as randn.m, based on the Box–Muller transform.
Finally, we obtained 1000 sets of 253 vectors in the [a, e, sin (i)]
domain, following the given tri-variate Gaussian distribution, using

X = μ + AZ. (9)

These vectors have the desired distribution because of the affine
transformation property; see Press et al. (2001) for information on
the Choleski decomposition and see the MATLAB support website for
information on randn.m. To check if the current observed asteroidal

7 We limit our analysis to the [a, e, sin (i)] proper-element domain be-
cause the g and, to a lesser extent, s frequencies are strongly perturbed in
the proximity of mean-motion resonances, such as 2J:-1A and 3J:-1A (see
Fig. 3). While there are techniques to normalize the dependence of the g fre-
quency as a function of the mean-motion resonance or of the semimajor axis
(Carruba & Michtchenko 2009), such a normalization requires an in-depth
study of local family membership, which we feel is beyond the purposes of
this study.

population follows a tri-variate normal distribution, we performed
two tests. First, we computed the number of generated fictitious
asteroids per cell, using the cells defined in Section 3. Secondly, we
performed a Pearson χ 2 test by computing a χ 2-like variable defined
by equation (5), where nint is the number of intervals of our variables
(20 in our case), qi is the number of real objects in the ith interval
and pi is the number of simulated objects in the same interval. Once
again, we remind the reader that values of χ 2 of the order of the
number of intervals are consistent with a good fit. Even without
considering the cells with a94 < a < a21, whose large asteroidal
population might be caused by the larger inclination values of the
ν6 resonance in the region (more asteroids could come to these areas
from low-inclinations regions compared to other cells), we find that
the lowest value of χ 2 for regions with a < a94 is 73.22 for 16
intervals for a given set of the simulated 253 asteroids. Using the
standard χ 2 probability distribution (Press et al. 2001), we find that
the null hypothesis has a probability of just 5.3 × 10−10, which is
low enough for us to be able to reject it safely.

In order not to incur possible binning problems associated with
our choice of cells, we also performed the Mardia test (Mardia
1970) on the multivariate normality for the whole 253 asteroid
population.8 The Mardia test is based on multivariate extensions
of skewness and kurtosis measures. For a sample (x1, . . . , xn) of
p-dimensional vectors, we can compute the A and B parameters
given by

A = 1

6n

n∑
i=1

n∑
j=1

[
(xi − μ)T�−1(xj − μ)

]3
(10)

and

B =
√

n√
8p(p + 2)

×
{

1

n

n∑
i=1

[
(xi − μ)T�−1(xi − μ)

]2 − p(p + 2)

}
. (11)

Here, �−1 is the inverse of the covariance matrix given by equa-
tion (7). Under the null hypothesis of multivariate normality, the
statistic A will have approximately a χ 2 distribution with (1/6)p(p +
1)(p + 2) degrees of freedom (10 for p = 3) and B will be approxi-
mately standard normal with a zero mean and a standard deviation
equal to 1. We computed A and B for our sample of 253 asteroids
and we obtained A = 166.98 and B = 7156.20, which have prob-
abilities lower than 10−6 of being associated with χ 2 and normal
distributions. Once again, we can safely assume that, on the whole,
large highly inclined asteroids do not follow a tri-variate Gaussian
distribution.

What about smaller samples of asteroids? To check this hypothe-
sis, we again looked at asteroids in the same range of a, e and sin (i)
as used in Section 4.1. We found that all asteroid samples have
values of the Mardia test B parameter that are too high to satisfy the
null hypothesis. However, the value of the A parameter for asteroids
in the regions of Phocaea, Hansa and Pallas was below 10. This
yielded a non-negligible probability (33 per cent for the lowest val-
ues of A of 7.7 in the Hansa region) that asteroids in the regions can
be fitted by a tri-variate normal distribution. We also found a rela-
tively small value of A for asteroids at very large inclinations (A =
12.8). While the large values of the B parameters should favour

8 Other tests on multivariate normality, such as the tri-dimensional KS test
(Fasano & Franceschini 1987), the Cox–Small test (Cox & Small 1978) and
others, are also adequate choices for testing multivariate normality.

C© 2011 The Authors, MNRAS 418, 1102–1114
Monthly Notices of the Royal Astronomical Society C© 2011 RAS

 at FundaÃ
§Ã

£o C
oordenaÃ

§Ã
£o de A

perfeiÃ
§oam

ento de Pessoal de N
Ã

­vel Superior on July 15, 2013
http://m

nras.oxfordjournals.org/
D

ow
nloaded from

 

http://mnras.oxfordjournals.org/


1110 V. Carruba and J. F. Machuca

the exclusion of the null hypothesis, the fact that large asteroids
in the same range of the semimajor axis are more likely to belong
to tri-variate normal distributions could indicate that the mecha-
nism that populated the highly inclined regions might have acted
‘vertically’ from lower values of inclinations to higher values, rather
than ‘horizontally’, moving asteroids originally at the same incli-
nations to a different semimajor axis – with the possible exception
of very highly inclined objects with sin (i) > 0.45. Further stud-
ies on such a mechanism are needed before this hypothesis can be
confirmed.

5 D ENSITY MAPS

In Section 4, we analysed the statistical likelihood that fluctua-
tions of the orbital distributions of asteroids at high inclination
might be caused by fluctuations in the one- and three-dimensional
statistical distributions. Here, we would like to address the dynam-
ical causes of underpopulated and overpopulated areas by asking
whether the observed underpopulated regions are caused by any
dynamical mechanism.

Information on the dynamics of asteroids can be obtained by an
analysis of the number densities of asteroids. Fig. 5 displays contour
plots of such densities. See Carruba (2009b) for a description of the
procedure used to generate such plots; here, we have used 70 steps of
0.02 au in a and 30 steps of 0.04 in sin (i). Higher number densities
of objects are shown in whiter tones. Once again, we divided the
population of objects into small-eccentricity bodies (e < 0.175;
Fig. 5, panel A) and large-eccentricity bodies (e > 0.175; Fig. 5,
panel B).

As can be seen in Fig. 5 (panels A and B), underpopulated ar-
eas (in black) in the Phocaea region are clearly associated with the
border of the ν6 secular resonance. Therefore, they can be under-
stood within the framework of the interaction of the ν6 resonance
topology, which causes the asteroid eccentricity to increase, and of
planetary close encounters, which remove such particles in short
time-scales (see Carruba 2010a for an in-depth description of these
mechanisms). In the Pallas area, we can identify an underpopulated
area, especially visible at high eccentricities, between the 8J:-3A
and 5J:-2A mean-motion resonances and the ν5 and ν16 secular reso-
nances (roughly, this corresponds to values of sin i between 0.45 and

0.50). This area is shown by white boxes in Fig. 5, and corresponds
to cell 10 in Section 4, which has a low probability that the large
asteroidal population follows a Poissonian distribution. The under-
populated area in the region between 2.5 and 2.82 au in a and 0.36
and 0.42 in sin (i) is caused by the presence of the ν6 resonance.
The destabilizing effect of this resonance on local asteroids has
been studied in Carruba et al. (2011b) and in Carruba & Morbidelli
(2011a). Finally, in the Euphrosyne area, the whole region above
the ν6 resonance and between the 5J:-2A and 9J:-4A mean-motion
resonances is characterized by a remarkable low number density of
objects. We have identified the underpopulated stable areas in Fig. 5
(panels A and B), corresponding to cells 14, 15 and 16 in Section 4,
with white boxes. For comparison purposes, we also selected a re-
gion in a high asteroid population area in the Euphrosyne region.
This region is identified by a green box in Fig. 5, and corresponds
to the last four cells (17, 18, 19 and 20) in Section 4.

To investigate if these regions are dynamically stable or not, we
performed a numerical analysis, the results of which are discussed
in Section 6.

6 DYNAMI CAL MAPS

We start by investigating the region between the 8J:-3A and 5J:-
2A mean-motion resonances, previously identified in Section 3 as
a stable region with a low asteroid population. For this purpose,
we obtained the synthetic proper elements of 2160 test particles
integrated with swift_mvsf, the integrator from the SWIFT package
(Levison & Duncan 1994) modified by Brož (1999), so as to include
on-line digital filtering in order to remove all frequencies with pe-
riods less than 600 yr. We integrated the objects over 10 Myr. The
synthetic proper elements were obtained with the procedure de-
scribed in Knežević & Milani (2003), except for the proper fre-
quencies, which were derived with the frequency modified Fourier
transform (FMFT) of Šidlichovský & Nesvorný (1997) using the
following procedure. We eliminated from the spectra of the Fourier
transform of the equinoctial elements all planetary frequencies and
we then assigned as the proper frequency the largest value in the
spectra that was still observable, rather than fitting the time series
of � f and 	f of the oscillations in the (k, h) and (p, q) planes. The
test particles were integrated under the influence of all planets from

Figure 5. An [a, sin i(i)] density map of low-eccentricity (e < 0.175, panel A) and high-eccentricity (e > 0.175, panel B) asteroids. The white boxes identify
the underpopulated stable areas and the green box displays the stable sample area. The other lines have the same meanings as in Fig. 1.

C© 2011 The Authors, MNRAS 418, 1102–1114
Monthly Notices of the Royal Astronomical Society C© 2011 RAS

 at FundaÃ
§Ã

£o C
oordenaÃ

§Ã
£o de A

perfeiÃ
§oam

ento de Pessoal de N
Ã

­vel Superior on July 15, 2013
http://m

nras.oxfordjournals.org/
D

ow
nloaded from

 

http://mnras.oxfordjournals.org/


Distribution of highly inclined asteroids 1111

Figure 6. Synthetic proper-element map of the region of (183) Istria. The
vertical lines display the mean-motion resonances, the blue asterisks identify
real objects in the region, the red circles are associated with asteroids in linear
secular resonances and the green circles display bodies in non-linear secular
resonances.

Venus to Neptune (Mercury was accounted for as a barycentric cor-
rection of the initial conditions) and we referred the inclination and
the arguments of pericentre and node to the invariable plane of the
Solar system. The initial conditions of planets and particles were
computed on 2000 January 1. Our test particles covered the interval
between 2.7 and 2.83 au, and between 23◦ and 31◦, and we set
the values of initial eccentricity and other angles equal to those of
(183) Istria, the real object in the area with the lowest identification.
Finally, we used 27 intervals in a and 80 in i.

Fig. 6 shows the [a, sin (i)] projection of the synthetic proper
elements of our test particles (black dots). Regions unaffected by
resonances appear as regular vertical alignments of objects in this
map, while areas affected by mean-motion resonances appear as
void of asteroids (see the area near the 5J:-2A resonance). Objects
in secular resonances appear to form inclined alignments. Vertical
red lines show the locations of the main mean-motion resonances
in the region, as identified in Carruba (2010b). Red circles show
the objects to within 0.3 arcsec yr−1 from the centre of the linear
secular resonance in the region, ν5 = g − g5. Green circles show the
locations of the non-linear secular resonances in the area, ν5 − ν16 =
g − g5 − (s − s6) and ν16 + ν5 − ν6 = s − s6 + (g − g5) − (g −
g6). For simplicity, we do not show secular resonances involving
the Uranus frequency g7, such as the ν7 resonance, because these
are close in proper-element space to the resonances involving the
g5 frequency. Finally, real asteroids are identified by blue asterisks.

As can be seen, the region below the ν5 resonance appears to
be stable over the length of the integration (the role of non-linear
secular resonances is investigated in Section 7, yet the number of
real objects observed is very limited. Between sin (i) equal to 0.45
and 0.50, we only observe 21 objects, three of them with an absolute
magnitude lower than 12.

To check if our results are independent of the choice of initial
conditions and of the length of the integration, we also performed
two numerical experiments. We integrated 1530 particles in the
region with different initial conditions (the orbital elements of the
planets and test particles were computed on 2010 January 1) and
with the same initial conditions previously used, but for a longer
time-scale. We find no significant difference concerning lost and
unstable objects between the three simulations; however, there are

Figure 7. Synthetic proper-element map of the region of (705) Erminia.
The magenta circles are associated with bodies in the ν6 secular resonance.
The other symbols are the same as in Fig. 6.

small differences concerning �30 particles in secular resonances,
which have slightly different values of proper frequencies g and,
in some cases, of proper sin(i). Because our goal is to study the
long-term stability of test particles and not to focus on the details
of secular dynamics, we believe that our choice of initial conditions
and integration length should be reasonable.

We now turn our attention to the underpopulated area in the
Euphrosyne region. To investigate the dynamical stability of this
region, we integrated 3000 particles between 2.83 and 3.03 au
(roughly speaking, between the 5J:-2A and 9J:-4A mean-motion
resonances) and between 20◦ and 35◦ in inclination (corresponding
to values of sin (i) between 0.20 and 0.55), under the same condi-
tions as the previous simulation. The initial eccentricities and the
other angles of the test particles were those of (705) Erminia, the
lowest numbered object in the region.

Fig. 7 shows the synthetic proper-element map for this region.
We can see the line of asteroids in ν6 antialigned librating states
(in magenta) that are stable near the centre of the unstable region
as a result of the ν6 resonance (see Carruba & Morbidelli 2011a
for a description of the dynamics of asteroids in ν6 antialigned
librating states). As observed for the region around (183) Istria, the
region between the ν5 resonance (asteroids in this resonance are
shown as red circles) and the ν6 resonance appears to be stable, yet
underpopulated. For this region, we observe the possible role that the
ν5 − ν16, ν6 − ν16 and ν16 + ν5 − ν6 non-linear secular resonances
might play in the dynamical evolution of asteroids in the region.
The long-term effect of these resonances when non-gravitational
forces are present is discussed in Section 7.

Finally, we turn our attention to the highly populated area (the
‘test area’) in the Euphrosyne region discussed in Section 3. Using
the same set-up as for the previous runs, we integrated 4500 particles
between 3.03 and 3.325 au, and between 20◦ and 35◦. The initial
eccentricities and the other angles of the test particles were those
of (31) Euphrosyne. Fig. 8 shows the synthetic proper-element map
for this region. Particles in the ν6 resonance are shown as magenta
circles, those in ν5 are identified by red circles9 and particles in the

9 The large number of red circles near the 2J:-1A mean-motion resonance
is an effect of the perturbation caused by the resonance on the asteroid
frequency g. This does not necessarily indicate that the particles are in a
secular resonance.

C© 2011 The Authors, MNRAS 418, 1102–1114
Monthly Notices of the Royal Astronomical Society C© 2011 RAS

 at FundaÃ
§Ã

£o C
oordenaÃ

§Ã
£o de A

perfeiÃ
§oam

ento de Pessoal de N
Ã

­vel Superior on July 15, 2013
http://m

nras.oxfordjournals.org/
D

ow
nloaded from

 

http://mnras.oxfordjournals.org/


1112 V. Carruba and J. F. Machuca

3.05 3.1 3.15 3.2 3.25

0.36

0.38

0.4

0.42

0.44

0.46

0.48

0.5

0.52

0.54

0.56

(31) Euphrosyne

ν
5

ν
6

ν
16

ν
16

ν
16

+ν
5

ν
6

z
1

z
2

z
3

S
in

e 
of

 p
ro

pe
r 

in
cl

in
at

io
n

Figure 8. Synthetic proper-element map of the region of (31) Euphrosyne.
The magenta circles show bodies inside the z1 resonance, the yellow circles
show asteroids inside the z2 resonance and the cyan circles are associated
with objects inside the z3 resonance. The other symbols are the same as in
Fig. 6.

ν16 resonance are indicated as green circles. The two main unstable
areas in the Euphrosyne regions caused by secular resonances are
associated with the ν6 and ν16 resonances. We also show particles
to within 0.3 arcsec yr−1 of the z1 = ν6 + ν16 resonance (magenta
circles), of the z2 = 2ν6 + ν16 resonance (yellow circles), of the
z3 = 3ν6 + ν16 resonance (cyan circles) and of the ν16 + ν5 − ν6

resonance (green circles). Other secular resonances that are known
to affect the orbits of real asteroids are the ν5 + ν6 and 2ν5 + ν6

resonances (not shown for simplicity). We also observe the effects of
the 1J:3S:-1A and 2J:5S:-2A three-body mean-motion resonances.
Finally, real asteroids are identified by blue asterisks.

Two things can be noticed about the test area. First, the region
is significantly affected by the presence of the 2J:-1A mean-motion
resonance, which, with the exception of the population of asteroids
near the resonance centre (the stable Zhongguos and the unstable
Griquas populations; Roig, Nesvorný & Ferraz-Mello 2002), signif-
icantly destabilized objects near the resonance separatrix. Overall,
the region is less ‘stable’ than the two regions previously studied,

yet it hosts a much larger asteroid population than that found in
previous areas. Secondly, many of the objects found in the region
are at an inclination lower than that of the ν6 resonance centre. This
leads us to wonder if the vast majority of the sin (i) > 0.3 population
in the Euphrosyne region might actually be considered as ‘highly’
inclined. In Section 7, we investigate further the long-term effect
of the local web of resonances when non-gravitational forces are
considered.

7 LONG-TERM STA BI LI TY

To investigate the orbital evolution over a long time-span of particles
in the areas discussed in the previous sections, we numerically
integrated test particles with SWIFT-RMVSY, the symplectic integrator
of Brož (1999), which simulates the diurnal and seasonal versions
of the Yarkovsky effect. Because objects in the hole in the Pallas
region are mostly of S-type (Carruba 2010b), we used values of the
Yarkovsky parameters appropriate for such bodies (Carruba et al.
2003): a thermal conductivity K = 0.001 W m−1 K−1, a thermal
capacity C = 680 J kg−1 K−1, a surface density 1500 kg m−2,
a Bond albedo of 0.1, a thermal emissivity of 0.95 and a bulk
density of 2500 kg m−3. We used two sets of spin-axis orientations,
one with an obliquity of +90◦ and one with −90◦ with respect
to the orbital plane, and periods obtained by assuming that the
rotation frequency is inversely proportional to the object’s radius,
and that a 1-km asteroid had a rotation period of 5 h (Farinella,
Vokrouhlický & Hartmann 1998). We used particles with a radius
of 2 km, and no reorientations were considered, so that the drift
caused by the Yarkovsky effect was the maximum possible. We
integrated the 2160 test particles of the first simulation over 100 Myr,
the time-scale used to study the effect of proximity to the 3J:-
1A mean-motion resonance by Guillens et al. (2002), under the
gravitational influence of seven planets (Mercury was accounted
for as a barycentric correction of the initial conditions).

Fig. 9 shows the maximum permanence times for our particles
that remained in the hole 1 region as a function of the synthetic
proper elements computed in Section 6. The black dots show objects
with permanence times of less than 10 Myr while the red full circles
are associated with particles with permanence times longer than
100 Myr. Blue asterisks show the orbital location of real asteroids

Figure 9. Maximum permanence times for particles in the hole 1 region for asteroids with an initial obliquity of +90◦ (panel A) and −90◦ (panel B). The
black dots show objects with permanence times less than 10 Myr, while the red full circles are associated with particles with permanence times longer than
100 Myr. The other symbols are the same as in Fig. 1.

C© 2011 The Authors, MNRAS 418, 1102–1114
Monthly Notices of the Royal Astronomical Society C© 2011 RAS

 at FundaÃ
§Ã

£o C
oordenaÃ

§Ã
£o de A

perfeiÃ
§oam

ento de Pessoal de N
Ã

­vel Superior on July 15, 2013
http://m

nras.oxfordjournals.org/
D

ow
nloaded from

 

http://mnras.oxfordjournals.org/


Distribution of highly inclined asteroids 1113

in the region. Panel A is associated with asteroids with an initial
obliquity of +90◦ and panel B shows the fate of objects with initial
obliquities of −90◦. We note that most of the km-sized asteroid were
not lost during the integration. Objects lost at times between 10 and
100 Myr were mostly particles close to the boundaries of the 5J:-2A
and 8J:-3A mean-motion resonances. Overall, the region between
the 5J:-2A and 8J:-3A mean-motion resonances seems to be quite
stable for km-sized objects on a 100-Myr time-scale.

We then integrated the same 3000 test particles used for the dy-
namical map of the hole in the Euphrosyne region (see Section 6)
over 100 Myr with SWIFT-RMVSY. The set-up of the simulation was the
same used for the particles in the Pallas hole, but here we used pa-
rameters for the Yarkovsky force typical of C-type asteroids, which
are predominant in the Euphrosyne region (the same parameters as
for an S-type asteroid, but with a bulk density of 1500 kg m−3;
Carruba et al. 2003). See Gil-Hutton (2006) for a discussion of the
spectral types of asteroids in the Euphrosyne region.

Fig. 10 shows the maximum permanence times for our particles
that remained in the hole 2 region as a function of the synthetic
proper elements computed in Section 6. The symbols used are the

same as for Fig. 9. Just as for the simulation for the Pallas hole, we
observe that the layer of particles with permanence times of less
than 10 Myr is limited to the boundaries of the main resonances
in the region (5J:-2A, 7J:-3A, 9J:-4A and ν6), but we observe that
the rest of the region is stable for km-sized objects on a 100-Myr
time-scale.

Finally, we turn our attention to the test area defined in Section 3.
Fig. 11 shows maximum permanence times for the same 4500 par-
ticles in the area as a function of the synthetic proper elements
computed in Section 6, assuming Yarkovsky parameters typical of
C-type asteroids. The symbols used are the same as for Fig. 9. We
observe that the region is overall stable over a 100-Myr time-scale,
with the exceptions of asteroids in the 2J:-1A central island of sta-
bility and a few objects near the borders of the ν6, ν16 and 1J:3S:-1A
resonances.

Because our dynamical analysis has shown that the studied
underpopulated areas are dynamically stable, even when non-
gravitational effects are considered, we have to conclude that the
low population of objects in these regions is real, and is not caused
by fluctuations in statistical distributions, as seen in Section 4. The

Figure 10. Maximum permanence times for particles in the hole 2 region for asteroids with an initial obliquity of +90◦ (panel A) and −90◦ (panel B). All
symbols are the same as in Fig. 9.

3.05 3.1 3.15 3.2 3.25

0.36

0.38

0.4

0.42

0.44

0.46

0.48

0.5

0.52

0.54

0.56

(31) Euphrosyne

S
in

e 
of

 p
ro

pe
r 

in
cl

in
at

io
n

B

3.05 3.1 3.15 3.2 3.25

0.36

0.38

0.4

0.42

0.44

0.46

0.48

0.5

0.52

0.54

0.56

(31) Euphrosyne

S
in

e 
of

 p
ro

pe
r 

in
cl

in
at

io
n

A

Figure 11. Maximum permanence times for particles in the test region for asteroids with an initial obliquity of +90◦ (panel A) and −90◦ (panel B). All
symbols are the same as in Fig. 9.

C© 2011 The Authors, MNRAS 418, 1102–1114
Monthly Notices of the Royal Astronomical Society C© 2011 RAS

 at FundaÃ
§Ã

£o C
oordenaÃ

§Ã
£o de A

perfeiÃ
§oam

ento de Pessoal de N
Ã

­vel Superior on July 15, 2013
http://m

nras.oxfordjournals.org/
D

ow
nloaded from

 

http://mnras.oxfordjournals.org/


1114 V. Carruba and J. F. Machuca

low number of objects observed in these areas might therefore be a
fossil indication of phenomena that occurred in the early phase of
the formation of the Solar system.

8 C O N C L U S I O N S

(i) We have revised the knowledge of the synthetic proper ele-
ments of asteroids in the highly inclined region sin (i) > 0.3. We
have found 10 073 numbered objects from the AstDyS site between
2.0 and 3.4 au.

(ii) We have obtained synthetic proper-element maps for the two
candidate regions and a test region with a high number density of
asteroids. We have identified all the main two-, three- and four-body
mean-motion and non-linear secular resonances in these areas.

(iii) We have verified the possibilities that dynamically stable
regions can be explained as fluctuations of simple statistical one-
dimensional distributions, such as Poissonian, uniform and Gaus-
sian distributions. Our results show that, while asteroids in the
Phocaea region might possibly follow a Poissonian distribution in
inclinations, dynamically stable underpopulated regions cannot be
explained as simple fluctuations of these statistical distributions.
These must be caused either by disuniformities in the mechanism
that populated the highly inclined regions or by some mechanism
that depopulated the currently stable areas.

(iv) We verified the possibility that highly inclined asteroids
might follow a tri-variate Gaussian probability distribution. Once
again, we found that this is not likely. Multidimensional normality
tests show that it is highly unlikely that highly inclined asteroids, as
a whole, might follow a tri-variate normal distribution. Asteroids in
the same range of semimajor axis in the Phocaea, Hansa and Pallas
regions show a small but non-negligible probability of belonging to
a tri-variate normal distribution. This could indicate that the mech-
anisms that populated the highly inclined regions might have acted
‘vertically’ from lower values of inclinations to higher, rather than
‘horizontally’ (i.e. moving asteroids originally at the same inclina-
tions to different semimajor axes). The possible exception is very
highly inclined objects with sin (i) > 0.45. Further studies on such
mechanisms are needed before this hypothesis can be confirmed.

(v) We have obtained density maps of the region and identified
two areas that are characterized by low probabilities that the ob-
served distributions of large asteroids might follow a Poissonian
statistics. The first candidate area is in the central main belt (the
Pallas region), while the second is in the outer main belt (the Eu-
phrosyne region).

(vi) We have simulated the long-term dynamical evolution of test
particles in these regions when the Yarkovsky force is considered.
Despite the low number density of objects observed, the two areas
with low number density were both mostly characterized by large
permanence times, of the order of 100 Myr or more. These showed
no significant qualitative differences from the highly populated test
area.

Overall, our results show that the low number densities of as-
teroids in dynamically stable areas cannot be explained by sim-
ple fluctuations of Poissonian, uniform or Gaussian distributions.
Some other mechanism has to be responsible for the observed
non-uniformities. There are two hypotheses: (i) the current high-
inclination population was created by a non-uniform mechanism;
(ii) some as yet unknown dynamical effect depopulated the currently
dynamically stable regions.

The first hypothesis could be coherent with what was observed
by Ribeiro & Roig (2010) for the external main belt, where most of
background objects seem to be associated with a family, suggesting
that primordial bodies in the main belt are limited and scattered.
The second hypothesis could justify the paucity of objects observed
between the 5J:-2A and 9J:-4A mean-motion resonances, also no-
ticeable in the main belt with a low inclination (sin (i) < 0.3). How-
ever, there could be some difficulty explaining why the region of
low asteroid density in the Pallas area is adjacent to a region of rel-
atively high asteroidal populations at higher and lower inclinations.
In any case, much more work is needed in order to understand the
causes of the current ‘Emmenthal structure’ of the highly inclined
asteroids.

AC K N OW L E D G M E N T S

We thank the referee, Hagai Perets, for hints and suggestions that
significantly improved the quality of this work. We are grateful to
A. Morbidelli for discussions in which this work was conceived.
We would like to thank the São Paulo State Science Foundation
(FAPESP), which supported this work via the grant 06/50005-
5, and the Brazilian National Research Council (CNPq, grants
302183/2008-6 and 473345/2009-9). We are also grateful to the
Department of Mathematics of FEG-UNESP for allowing us to use
their facilities.

REFERENCES

Brož M., 1999, Thesis, Charles Univ., Prague, Czech Republic
Carruba V., 2009a, MNRAS, 395, 358
Carruba V., 2009b, MNRAS, 398, 1512
Carruba V., 2010a, MNRAS, 403, 1834
Carruba V., 2010b, MNRAS, 408, 580
Carruba V., Michtchenko T., 2007, A&A, 475, 1145
Carruba V., Michtchenko T., 2009, A&A, 493, 267
Carruba V., Morbidelli A., 2011, MNRAS, 412, 2040
Carruba V., Burns J. A., Bottke W., Nesvorný D., 2003, Icarus, 162, 308
Carruba V., Machuca F., Gasparino H., 2011, MNRAS, 412, 2052
Cox D. R., Small N. J. H., 1978, Biometrika, 65, 263
Farinella P., Vokrouhlický D., Hartmann W. K., 1998, Icarus, 132, 378
Fasano G., Franceschini A., 1987, MNRAS, 225, 155
Gil-Hutton R., 2006, Icarus, 183, 93
Guillens S. A., Vieira Martins R., Gomes R. S., 2002, AJ, 124, 2322
Knežević Z., Milani A., 2003, A&A, 403, 1165
Levison H. F., Duncan M. J., 1994, Icarus, 108, 18
Mardia K. V., 1970, Biometrika, 57, 519
Michtchenko T. A., Lazzaro D., Carvano J. M., Ferraz-Mello S., 2010,

MNRAS, 401, 2499
Milani A., Knežević Z., 1992, Icarus, 98, 211
Milani A., Knežević Z., 1994, Icarus, 107, 219
Morbidelli A., 2002, Modern Celestial Mechanics. Taylor and Francis, Lon-

don
Press V. H., Teukolsky S. A., Vetterlink W. T., Flannery B. P., 2001, Numer-

ical Recipes in Fortran 77. Cambridge Univ. Press, Cambridge
Ribeiro A. O., Roig F., 2010, in Fernandez J. A., Lazzaro D., Prialnik D.,

Schulz R., eds, Proc. IAU Symp. 263, Icy Bodies of the Solar System.
Kluwer, Dordrecht, p. 237

Roig F., Nesvorný D., Ferraz-Mello S., 2002, MNRAS, 335, 417
Šidlichovský M., Nesvorný D., 1997, Celest. Mech. Dyn. Astron., 65,

137
Yoshikawa M., 1987, Celest. Mech., 40, 233

This paper has been typeset from a TEX/LATEX file prepared by the author.

C© 2011 The Authors, MNRAS 418, 1102–1114
Monthly Notices of the Royal Astronomical Society C© 2011 RAS

 at FundaÃ
§Ã

£o C
oordenaÃ

§Ã
£o de A

perfeiÃ
§oam

ento de Pessoal de N
Ã

­vel Superior on July 15, 2013
http://m

nras.oxfordjournals.org/
D

ow
nloaded from

 

http://mnras.oxfordjournals.org/