1 Introduction

Maybe the secret of the huge success of Kolmogorov–Arnold–Moser (kam) theory relies in the spectacular application, found out by Vladimir Igorevich Arnold, to the planetary problem. Indeed, one decade after Kolmogorov’s announcement, at the 1954’s International Congress of Mathematician, of the “theorem of the conservation of the invariant torus” [14], the brilliant student of Kolmogorov—aged 27—formulated a version of Kolmogorov’s theorem (which he called the “Fundamental Theorem”) suited to the planetary problem [1]. He then used it to prove the “metric stability” of the simplest, albeit non-trivial, planetary system: two planets and a sun constrained on a plane. Strong degeneracies prevented a straightforward application of the Fundamental Theorem to the most general planetary system, which indeed was obtained in the subsequent 50 years, after those degeneracies were completely understood [4, 6, 15, 23, 24]; see [5] for a review. The success of kam theory in classical mechanics boosted other investigations, like instability or finite time stability [2, 21].

Despite the quality and the quantity of results of this kind, some intriguing questions remain open. In this paper we are particularly interested in the problems that have a kind of “global” structure in planetary systems, like, e.g., the mass distribution and/or the stability of the belts of many lightweight objects. In our Solar System this could be relevant to understand the global features of the various asteroids belts, the rings around the planets, or the space debris .

From this respect the basic idea of statistical mechanics, that is the possibility to substitute the exact knowledge of the dynamics of an N-particles system with a probability distribution on its dynamical status in a fixed instant, seems to be promising. Losing the detailed knowledge of the trajectory of the system in the phase space, one gains the possibility to describe global quantities like, e.g., pressure, temperature or density. With this attitude, the study of the analyticity of such global functions may give a microscopical justification of many very interesting phenomena, e.g., the phase transitions.

To fulfill this program, however, one has to assume various technical conditions, and one of the first constraints, understood from the very beginning of the discipline, is the so-called stability condition on the interacting potential. Namely, one has to impose to the potential V(x, y) of the interaction between particles the following condition: it has to exist a positive constant B such that

$$\begin{aligned} \sum _{1\le i< j\le N}V(x_i,x_j)\ge -BN \end{aligned}$$
(1)

for all the possible choices of the positions \(x_1,...,x_N\) of the N particles (see, for instance, [29] and [9]). Such stability condition, obviously, does not hold for the Newtonian potential. In the quantum case it has been possible to show that a non-stable potential, like the Coulomb one, leads to a stable behavior of the matter (see [16]). Many other results (see, e.g., [10, 17, 18]), have been proved exploiting explicitly the features of the Fermi and Bose statistics. In the study of classical systems with Newtonian interaction, a possible way to perform the thermodynamic limit is to define interactions depending on the number of particles, the so-called mean field systems, or, almost equivalently, to rescale suitably the energy with the number of particle in a fixed volume, the so-called Vlasov limit [12, 13, 19]. This approach proves useful to model, for instance, plasma physics and astronomical globular systems. In a different context, namely that of growing planets, some analytical results have been obtained recently by [30] in the case of the planetary systems in which the instability is part of the desired result. In [26, 27] the statistical mechanics of dense stellar clusters have been studied, outlining its similarity with liquid crystals. Moreover, numerical simulations of gravitating systems represent a very active research field. Finally, other recent attempts to perform, numerically, a statistical analysis of the future planet orbits in the solar system include, e.g., [20].

In this paper we propose a different approach, not yet previously investigated, at least to our knowledge. We study a model of interacting particles in a planetary system. Our estimates are mathematically rigorous, although some of the assumptions underlying the definition of the model are motivated by heuristic considerations based on the observed features of the solar system. Some of these features are suitably simplified. The first simplification is the fact that our planetary system is assumed to be planar, and the central star is fixed. It will be clear that this assumption simplifies the details but it does not modify the structure of the bounds we will present. Secondly, we deliberately decide to lose details in the description of the interaction between the center, called hereafter the star, and each light particle, call hereafter asteroid. In particular, we describe each orbit as a probability distribution around a fixed circular orbit. Such probability distribution does not fix the energy of each asteroid. In other words, we try a description of the system in terms close to the familiar idea of canonical distribution. From a physical point of view, this assumption can be justified thinking, for instance, to the main belt of asteroids: the 2-body elliptical trajectory is actually an approximation due to the fact that each asteroid is perturbed by the planets. Hence, the energy of the single asteroid is not conserved. We substitute the computation of the actual trajectory, perturbed by the planets, with the probability distribution mentioned above. The distribution described so far play the role of the reference, or free, measure, in the sense that each asteroid has its independent (i.e., factorized) free measure. Then we introduce an interaction in the probability distribution, adding a gravitational potential among asteroids. As far as the short distance configurations are concerned, a regularization in terms of hard core interaction among asteroids is introduced: if \(x_i\) represent the vector position of the i-th asteroid and \(a_i\) is its radius, then \(|x_i-x_j|\ge a_i+a_j\), or, in other way, the potential \(V_{ij}\) is infinite if \(|x_i-x_j|\le a_i+a_j\). Note that in this context the planets are much bigger than the interacting asteroids, but are very far. The interacting asteroids, on the other side, are light but in principle they can have very small mutual distances, exactly of the order of the sum of their radii, and these colliding configurations will give a huge contribution to the interacting probability measure. We are not assuming in this model neither agglomerations nor disintegrations due to these collisions, we simply compute their static contributions to the interacting probability measure. Even with this strong simplifications, however, it is hopeless to perform the thermodynamical limit. The hard core interaction prevents the possibility of configurations having an infinite probabilistic weight, but the features of the Newtonian interaction, and in particular its very slow decay, do not allow an estimate of the form (1). Nevertheless, assuming that the number of asteroids N is a large but finite parameter, and discussing its value in terms of the masses of the asteroids, we find results that are quite interesting in terms of the description of the real Solar System.

In order to quantify the effect of the gravitational interaction on the trajectories of the asteroids, we define a notion of stability in the following way: each asteroid, with respect to its independent reference measure, has its own variance of the distance from the star. Call \(\sigma ^2_0\) such variance. If the number N of asteroids (and their masses distribution) is appropriately chosen, it is possible, uniformly in the choice of the asteroid, to give for the interacting measure an estimate of the variance \(\sigma ^2\) of the form:

$$\begin{aligned} \sigma ^2=\sigma ^2_0(1+\varepsilon ) \end{aligned}$$
(2)

we then say that the system is \(\varepsilon -N\)-stable. If \(\varepsilon \) is sufficiently small, we can argue that the interaction among asteroids implies small modifications of the asteroid’s orbit.

To give an initial idea of the smallness required on \(\varepsilon \) in order to have an astronomical interpretation of this notion of stability, we are assuming that the energy of each asteroid confines it in the vicinity of the minimum of its effective potential. This means that the eccentricity is small, together with the effect of the interaction with the distant planets, and hence the reference (free) probability of each asteroid has a standard deviation around its average radius very small with respect to the radius itself. If the effect of the interaction among asteroids is such that this standard deviation stays small, i.e., if \(\varepsilon \) is of order 1 or less, we can assume that the probability of a large deviation of an asteroid from its average radius is very small. A quite important point is the relation between this large deviation probability and the astronomical stability of the orbit. The inverse of the probability of a large deviation multiplied by the time scale of the variations of the distance between the asteroid and the star, that is, at least the period of the orbit of the asteroid itself, gives us an idea of the time scale in which we have to assume that such large deviation is not realized, and hence the asteroid will remain close to its present orbit. Such large deviation probabilities estimates could be done with rough but robust tools, like Chebyshev inequality, or with more sophisticated techniques. This delicate point will be discussed in the next section.

We study in detail three different setups:

(1) Similar asteroids

Our first setup is somehow theoretical: the asteroids are very light and their masses are comparable. The average radius of the orbit is similar for each asteroid; we will compute the stability of the system in the worst case, i.e., for equal average radii. Fixing the radius of the asteroids, we find an estimate, depending on \(\varepsilon \), of the maximum value of the number of asteroids N such that the system is \(\varepsilon -N\)-stable. This first computation is important in order to understand that in this context it is not possible to perform a standard thermodynamical limit. However, if the total mass of asteroids depends on their number N, and it goes suitably to zero when N increases, then it is possible to prove the \(\varepsilon -N\)-stability of the system. In this context it is easier to outline the basic problem that one has to face: we want to fix the parameter of the system in such a way that the contributions of the collisions, in which the Hamiltonian is negative and has a large absolute value, are not too relevant for the canonical probability distribution.

(2) Asteroids with a given mass distribution

In this setup the asteroids have always a similar average distance from the star, but they have a well-defined distribution of the masses. We show that a distribution of the form:

$$\begin{aligned} N(>r)=\frac{c}{r^\nu } \end{aligned}$$
(3)

where \(N(>a)\) is the number of asteroids with radius greater than a, \(\nu >1\) and c is a suitable positive constant, guarantees that N can be chosen quite large, and yet the system remains \(\varepsilon -N\)-stable. Remarkably, our assumption about the radii distribution of the asteroids seems to be quite close to the observed one. In particular, the estimates obtained from observed data give a value of \(\nu \) between 1.3 and 3 (see [28] and references therein). Again, the basic problem is to control the contributions of the collisions.

(3) Planets with well-separated orbits.

In the last part of this paper we try to apply the same techniques developed for asteroids to a system of planets, i.e., of object small with respect to the star but larger than asteroids, having orbits with very different average radii: we show that assuming for the average radius of the orbit of the i-th planet the following Titius–Bode law:

$$\begin{aligned} R_i=b+c a^i \end{aligned}$$
(4)

with \(R_i\) the orbit’s average radius of the i-th planet, b and c fixed lengths (0.4 and 0.3 U.A., respectively, for the Solar System) and \(a>1\) a fixed number (\(a=2\) for the Solar System), and assuming N small enough, the system is \(\varepsilon -N\)-stable. In this context we briefly discuss the 1-stability of the Galilean Jupiter’s satellites. Also in this case the main problem is to control the contributions due to the collisions. Here, however, we have to exploit the fact that in order to have a collision the planets have to deviate substantially from their reference distribution, see further comments below, in the beginning of Sect. 5.

In order to obtain these results, we have to interpret the classical meaning of the thermodynamic constants in Gibbs distribution in a different way. The main problem is the interpretation of the temperature in this context. As it will be clear in the next section, it is physically meaningless to define a common temperature for objects with different masses. Since the exact Keplerian orbit and also its correction due to external objects are computed in terms of gravitational interaction, the shape of the orbit and, consequently, the form of the probability distribution that we want to define may not depend on the mass of the orbiting object. On the other side, the contribution to the energy of such object depends linearly on the mass. Hence the contribution of an object of mass m to the probability distribution has to be rescaled by a factor 1/m. This means that each object has its own “temperature”, proportional to 1/m. After this rescaling, it is natural to assign the role of the temperature to a number \(\gamma \) that is related to the free measure deviation of the radii of the asteroids. Small temperature, corresponding in our model to large \(\gamma \), means that the asteroid has a free distribution concentrated in the vicinity of its reference circular orbit. Hence we are assuming that in the low-temperature regime the eccentricity of the orbits and the interactions with the heavy far planets are small, and the asteroids, as far as their free measure distribution is concerned, have an average energy very close to the minimum of their effective potential.

Moreover our particles are obviously distinguishable, and hence the combinatorial Gibbs factors \(\frac{1}{N!}\) are absent in our treatment. In order to estimate the deviations of the radii of the asteroids (and, eventually, planets) in the presence of the interaction among them, we had to use a procedure quite standard in statistical mechanics, usually known as Peierls argument, see Sect. 2, and then judicious combinatorial estimates, similar to the ones introduced in cluster expansion. Such estimates, due to the absence of the Gibbs factor, have some non-standard features.

The results we obtain, despite the simplicity of the estimate we present, may have some interest. From a quantitative point of view, the estimates of the mass and of the number of the asteroids and of the planets are quite different from the ones observed in the Solar System, but the orders of magnitude are not too distant. In the simpler case of Galilean satellites our notion of stability is guaranteed for masses of the satellites close to the actual ones. Moreover the model explains why in order to have stability the number of very light asteroids may be relatively high, while the planets have to be quite far apart and their number has to be very small, of the order of \(N\le 10\).

The model, then, seems to have a reasonable fit with the observed data.

The work is organized as follows: in Sect. 2 we present our model, we define more precisely the notion of “thermodynamical” stability for planetary systems and we define the relation among this notion of stability and the astronomical one; in Sect. 3 we discuss an application of the model to a belt of asteroids having a very narrow distribution of masses; in Sect. 4 we generalize the same results to a more realistic asteroid belt; Sect. 5 is devoted to the application of our model to a planetary system in which the radii of the planets satisfy a kind of Titius–Bode law. Finally in Sect. 6 we discuss some brief final remarks.

2 The Model

2.1 Planetary System

Consider a system of N bodies with mass \(m_i\), constrained on a plane, with pairwise gravitational interaction and interacting gravitationally with a much larger body, the star, of mass M centered at the origin of a reference frame in the plane. The system is described by the Hamiltonian

$$\begin{aligned} H(\vec {p}, \vec {q}) = \sum _{i = 1}^{N} \frac{|p_{i}|^{2}}{2 m_i} - \sum _{i=1}^{N} \frac{k M m_{i}}{\left| q_i\right| } - \sum _{1\le i<j\le N} \frac{k m_{i} m_{j}}{\left| q_{i}-q_j\right| } \end{aligned}$$
(5)

where \(q_i\) are two-dimensional Euclidean coordinates, \(p_i\) the corresponding moments and k the gravitational constant. We remark that in our model the mass M does not move. This appears in (5) from having neglected centrifugal terms coming from taking the reference frame centered at M; compare, e.g., [6], for the general expression of the N-body Hamiltonian in the star-centered frame.

Calling

$$\begin{aligned} H_{0}(\vec {p}, \vec {q}) = \sum _{i = 1}^{N} \frac{|p_{i}|^{2}}{2 m_i} - \sum _{i=1}^{N} \frac{k M m_{i}}{\left| q_i\right| } \end{aligned}$$
(6)

the Hamiltonian describing N uncoupled central interactions with the star, the original Hamiltonian can be seen as the sum of \(H_{0}\) and a perturbing term.

Rewriting each term \(h_0\) of the sum appearing in \(H_{0}\) using polar coordinates (with \(\rho \) the distance from the star and \(\theta \) the true anomaly) and setting \(p_{\theta } = J\) for the conservation of the angular momentum in the central system, we can write \(h_0\) as

$$\begin{aligned} h_{0} = \frac{p_{\rho }^{2}}{2m}+V_{{{\,\mathrm{eff}\,}}}(\rho ) = E \end{aligned}$$
(7)

where

$$\begin{aligned} V_{{{\,\mathrm{eff}\,}}} (\rho )= \frac{J^2}{2m\rho ^{2}} - \frac{kMm}{\rho } \end{aligned}$$
(8)

can be interpreted as an effective potential. If the total energy of the system is close to the minimum of \(V_{{{\,\mathrm{eff}\,}}} (\rho )\), it makes sense to think that a second-order approximation of this potential (harmonic potential) describes reasonably well the gravitational interaction with the star.

Denote by \(R = \frac{J^2}{k m^2 M}\) the value at which the minimum of the potential is attained. A straightforward computation gives, introducing the dimensionless coordinate \(\xi ={\frac{\rho - R}{R}}\), that \(V_{{{\,\mathrm{eff}\,}}} (\rho )\) can be rewritten in terms of \(\xi \) as

$$\begin{aligned} V(\xi ) = \frac{1}{2} \frac{k M m}{R}\left( -1+ \frac{\xi ^2}{(1+\xi )^2}\right) . \end{aligned}$$
(9)

We will call the expansion of this potential in which we neglect the unessential constant \(- \frac{1}{2} \frac{k M m}{R}\) and we keep only the second-order term:

$$\begin{aligned} V_2(\xi ) = \frac{1}{2} \frac{k M m}{R}\ \xi ^2 \end{aligned}$$
(10)

the Gaussian approximation of the central interaction.

Remark 2.1

The Gaussian approximation is apparently a strong assumption, so we need a pair of comments. On one side, neglecting the first term in (9) reflects the precise choice of regarding the R’s as fixed quantities, rather than as thermodynamical variables (see also the next section). Secondly, for what concerns the approximation of V with its quadratic expansion, it will be clear (see Sect. 5 for a discussion) that in the applications of our model to system of very small bodies (asteroids) such assumption is reasonable, because we will show that the interaction between the asteroids keeps the variance of \(\xi \) of the same order of the unperturbed system, and the main terms in the corrections are related to configurations with small \(\xi \), namely colliding asteroids.

2.2 Free Probability Measure

Denote by \(R_{i}\) the radius of the i-th circular orbit, and consider the Gaussian approximation of its central potential.

$$\begin{aligned} V_{2,i}(\xi _i) -=\frac{1}{2} \frac{k M m_i}{R_i}\ \xi _i^2 \end{aligned}$$
(11)

To avoid heavy notations, we will denote such potential with \(V_i(\xi _i)\), dropping the subscript 2 until further notice. We want to define a reference probability measure on the position of the body in the plane in the absence of perturbations. Recalling that \(\xi _i\) and \(\theta _i\) are, respectively, the dimensionless deviation of the distance from the mean radius and the true anomaly of the i-th body, we consider the probability measure

$$\begin{aligned} {{\,\mathrm{d\mathrm {\mu }}\,}}_0(\xi _i, \theta _i) = \frac{e^{\beta _i V_i(\xi _i)}{\text {d}}\!{\theta _i}{\text {d}}\!{\xi _i}}{\int _{0}^{2\pi } {\text {d}}\!{\theta _i} \int _{-\infty }^{\infty }{\text {d}}\!{\xi _i}\ e^{-\beta _i V_i(\xi _i)}}. \end{aligned}$$
(12)

where \(\beta _i\) is a positive parameter. Note that the kinetic part in the Hamiltonian does not play any role in the probability measure since it appears, as a factor, both at the numerator and at the denominator. A similar fate would hold for the terms \(-\frac{1}{2}\frac{kMm_i}{R_i}\) coming from (9).

In statistical mechanics the parameter \(\beta \) plays the role of the inverse temperature. When the inverse temperature is large, the system tends to remain close to the local minimizers of the Hamiltonian. Here each \(\beta _{i}\) is determined so to have a probability distribution on the unperturbed system having a standard deviation of the distance between the asteroid and the star much smaller than its average value. Since we assume small deviations from the average radius of the orbit, the harmonic approximation of the effective gravitational potential is reasonable. The temperature \(\beta _i\) has the dimension of the inverse of an energy. As already remarked, \(\beta _i\) has to be rescaled by the inverse of the mass \(m_i\) in order to be physically meaningful, see the previous section. Moreover, the fact that our coordinates \(\xi _i\) are dimensionless suggests to rescale the temperature by a factor \(R_i\). Hence we choose

$$\begin{aligned} \beta _i = \frac{R_i}{k m_i M} \gamma _i^2 \end{aligned}$$
(13)

where \(\gamma _i\) is a sufficiently large pure number, in order to have a small variance \(\sigma ^2(\xi _i)\) of the deviation \(\xi _i\) from the average radius. As outlined in the introduction, \(\gamma _i\) takes into account both the eccentricity and the interaction with planets.

Introducing (13) in (12) yields:

$$\begin{aligned} {{\,\mathrm{d\mathrm {\mu }}\,}}_0(\xi _i, \theta _i) = \frac{\gamma _i}{(2\pi )^\frac{3}{2}} e^{-\frac{1}{2} \gamma _i^2 \xi _i^2}{\text {d}}\!{\theta _i}{\text {d}}\!{\xi _i} \end{aligned}$$
(14)

meaning that the measure of \(\xi _i\) (without perturbations) is Gaussian with zero mean and variance \(\sigma _i^2 = \frac{1}{\gamma _i^2}\).

Note that as a consequence of the previous considerations, in this model, each body has its own “temperature” \(\beta _i\) that tunes its interaction with the star.

2.3 Interacting Probability Measure

When the interaction between the asteroids is taken into account, the probability distribution of the system is proportional to \(e^{-H^{\mathrm {ad}}}\) where \(H^{\mathrm {ad}}\) is the dimensionless Hamiltonian

$$\begin{aligned} H^{\mathrm {ad}}(\vec {\xi }, \vec {\theta }) = - \sum _{i=1}^{N} \frac{1}{2} \gamma _i^2 \xi _i^2 - \sum _{1\le i<j\le N} \beta _{ij} V_{ij} \end{aligned}$$
(15)

with

$$\begin{aligned} V_{ij} = \frac{k m_i m_j}{\left| \vec {x}_i - \vec {x}_j\right| }; \quad \vec {x}_i = (R_i(1 + \xi _i)\cos \theta _i; R_i(1 + \xi _i)\sin \theta _i) \end{aligned}$$
(16)

and each \(\beta _{ij}\) is a parameter tuning the interaction between the i-th and the j-th body. Here the rescaling of \(\beta _{ij}\) is not obvious, as in the case of the free measure, because \(V_{ij}\) appears in the dynamics of both bodies. In order to have a measure related to the actual interaction among asteroids, it seems very reasonable that the dimensionless expression appearing eventually in our measure should respect the following two conditions

(1) The strength of the interaction has to be of the order of m/M, where m is some kind of average between the masses of the bodies i and j, as suggested by the gravitational nature of the interactions.

(2) The expression of the potential in terms of dimensionless units \(\xi _i\) and \(\xi _j\) should be rescaled by a factor proportional to R, where again R is some kind of average between the radius of the bodies i and j.

To fulfill both request, we argue as follows. Other choices, fulfilling (1) and (2), would affect only the constants appearing in the subsequent estimates.

First of all, order the indices of the asteroids according to their average distance from the star, i.e., say that \(R_i\le R_j\) if \(i<j\). Consider the asteroids with indices i and j (with \(i < j\)) and consider the case where \(\xi _i = \xi _j = 0\). In other words it means that the two planets have both distance from the star equal to the radii \(R_i\) and \(R_j\), respectively. Consider then the scenario where \(\xi _i = \frac{1}{\gamma _i},\ \xi _j = \frac{1}{\gamma _j}\): each of the two asteroids has been moved away from the star by an amount equal to a free standard deviation. Call \(\Delta V_i\) and \(\Delta V_j\) the variation of the gravitational potential describing the interaction of the two bodies with the star associated with this change of scenario, and let \(\Delta V_{ij}\) the corresponding change in the potential describing the gravitational interaction among the two planets.

In order to have a probabilistic weight due to the interaction among asteroids that is comparable with the one due to interaction with the star, we want that, when considering this change of scenario, the ratio \(\frac{\Delta V_{ij}}{\Delta V_i + \Delta V_j}\) is the same as the ratio of the corresponding variation in the exponent of \(e^{-H^{\mathrm {ad}}}\), that is we want that

$$\begin{aligned} \frac{\Delta V_{ij}}{\Delta V_i + \Delta V_j} = \frac{\beta _{ij} \Delta V_{ij} }{\frac{1}{2}\gamma _i^2 (\frac{1}{\gamma _i})^2 + \frac{1}{2}\gamma _j^2 (\frac{1}{\gamma _j})^2 } = \beta _{ij} \Delta V_{ij}. \end{aligned}$$
(17)

Hence

$$\begin{aligned} \beta _{ij}= \frac{1}{\Delta V_i + \Delta V_j}=\frac{1}{kM}\frac{R_i(1+\gamma _i)R_j(1+\gamma _j)}{m_iR_j(1+\gamma _j)+m_jR_i(1+\gamma _i)} \end{aligned}$$
(18)

This means that

$$\begin{aligned} \beta _{ij} V_{ij} = \gamma _{ij} \frac{\sqrt{R_i R_j}}{\left| \vec {x}_i - \vec {x}_j\right| }:= \gamma _{ij} \frac{r_{ij}}{\left| \vec {x}_i - \vec {x}_j\right| } \end{aligned}$$
(19)

with

$$\begin{aligned} \gamma _{ij} = \frac{m_i m_j}{M}\frac{\sqrt{R_iR_j}(1+\gamma _i)(1+\gamma _j)}{m_iR_j(1+\gamma _j)+m_jR_i(1+\gamma _i)}\end{aligned}$$
(20)

The statistical mechanics model that we investigate is, therefore, defined through the following (dimensionless) Hamiltonian (calling again the dimensionless Hamiltonian and the dimensionless potential H and V, respectively, with an abuse of notation):

$$\begin{aligned} H(\vec {\xi }, \vec {\theta }) = \sum _{i=1}^{N} \frac{1}{2} \gamma _{i}^2 \xi _i^2 - \sum _{1\le i<j\le N} \gamma _{ij} \frac{r_{ij}}{\left| \vec {x}_i - \vec {x}_j\right| }=\sum _{i=1}^{N} \frac{1}{2} \gamma _{i}^2 \xi _i^2 - \sum _{1\le i<j\le N} V_{ij} \end{aligned}$$
(21)

Recall that the \(\vec {x}_i\) are constrained by the hard core compatibility condition \(\left| \vec {x}_i - \vec {x}_j\right| \ge a_i + a_j\) where \(a_i\) is the radius of the i-th body. Further note that assuming the asteroids to have constant density \(\delta \), we have \(m_i = \frac{4}{3}\pi \delta a_i^3\).

The probability measure induced by the Hamiltonian (21) that we want to take into account to describe the planetary system is, therefore,

$$\begin{aligned} \mu ({{\,\mathrm{\cdot }\,}}) = \frac{\int {\text {d}}\!{\vec {\xi }} \int {\text {d}}\!{\vec {\theta }} \, ({{\,\mathrm{\cdot }\,}}) \, e^{-H(\vec {\xi }, \vec {\theta })} }{\int {\text {d}}\!{\vec {\xi }} \int {\text {d}}\!{\vec {\theta }} \, e^{-H(\vec {\xi }, \vec {\theta })} }. \end{aligned}$$
(22)

Interpreting the variance of each \(\xi _{i}\) as a quantity linked to the eccentricity of the i-th orbit, assessing the stability of the system amounts to control the variance of the \(\xi _{i}\)’s.

In particular we want to determine the conditions on N and on the physical parameters (mass, radius of the orbits) for which the system is stable in the sense of the following:

Definition 2.2

The system (5) is called \(\varepsilon -N\)-stable if, for a fixed \(\varepsilon \) and for all \(i = 1, \ldots , N\)

$$\begin{aligned} \langle \xi _i^2\rangle \le (1 + \varepsilon )\langle \xi _i^2\rangle _{0} \end{aligned}$$
(23)

where \(\langle \xi _i^2\rangle _{0}\) is the variance of \(\xi _{i}\) with respect to \({{\,\mathrm{d\mathrm {\mu }}\,}}_{0}\).

Indeed, if the previous condition is satisfied for an \(\varepsilon \) small enough, the deviations of the radii of the orbits of the asteroids, with respect to the orbits they would have if the other asteroids were not there, stay small.

As outlined in the introduction, here there are a couple of delicate points deserving a discussion. First, what is the relation among this definition of stability and the evaluation of the stability of the orbit in an astronomical sense? Assume that an asteroid has a stable orbit for a time T if its deviation from its average radius \(R_i\) stays smaller than \(A R_i\) for all \(t<T\), with A a suitable constant. Then we can deduce T by the inverse of the probability \(P(\xi _i>A)\) times the period of revolution \(\tau \) of the asteroid around the star. In probability theory there are many ways to estimate \(P(\xi _i>A)\). One of the rougher, needing just the control of the variance of \(\xi \), is the Chebyshev inequality. A direct application of such inequality gives, for \(A=1\), \(T\approx \gamma _i^2 \tau \). This kind of estimates would hold for objects with masses comparable with the real masses of the asteroids in the main belt. However, this time (order of ten thousand years) is quite short in an astronomical sense. Our control of the variance, as it will turn out in the following sections, allows in principle to use different estimates for \(P(\xi _i>A)\): for instance using the fact that the reference measure is Gaussian, one could use Chernoff inequality, or other methods involving the detailed control of higher order moments of the distribution. This could be easily done in principle, but it would strictly rely on the details of the reference measure. Recall that the standard statistical mechanics is based much more on the geometrical properties of the N-dimensional space, with N of the order of the Avogadro’s number, that on the details of the Gibbs probability measure. While the Chebyshev estimate mentioned above is quite robust, depending only on the fact that a reasonable reference probability should be strongly concentrated around its circular orbit, a more refined estimate will require some additional arguments regarding the faithfulness of the free reference measure. This will be the subject of further investigations.

The second point deserving a discussion is the role of the collisions in the evaluation of the canonical measure defined in (22). As it will be clear by the computations of the following subsection, the main problem of this approach is the control of the probability of the configuration in which n asteroids are very close one with the other (n body collisions), because the corresponding energy turns out to be negative and proportional to \(n^2\). This gives rise to a probabilistic weight proportional to \(C^{n^2}\), with \(C>1\), and it is not clear how to control it from a combinatorial point of view. This seems to suggest that the collision terms are the leading one, and then only a very small amount of asteroids may be considered in order to have \(\varepsilon -N\)-stability with \(\varepsilon \) reasonably small. A relatively standard but judicious control of the structure of the measure (22), however, seems to indicate that the number of asteroids that can be taken into account in this model is not too different from the actual one (see Sect. 4 below), and that this interpretation may suggest the fact that the number of asteroids that we see today is the relic of a much bigger initial set, in which a large part of asteroids has been lost due to the intrinsic instability of the system.

2.4 Estimation of \(\langle \xi _{m}^{2}\rangle \)

The value of \(\langle \xi _{m}^{2}\rangle \) is given by

$$\begin{aligned} \langle \xi _{m}^{2}\rangle&= \frac{\int {\text {d}}\!{\vec {\xi }} {\text {d}}\!{\vec {\theta }}\ \xi _m^2 e^{-H(\vec {\xi }, \vec {\theta })} }{\int {\text {d}}\!{\vec {\xi }} {\text {d}}\!{\vec {\theta }} e^{-H(\vec {\xi }, \vec {\theta })} } = \frac{\int {{\,\mathrm{d\mathrm {\mu }}\,}}_0(\vec {\xi }, \vec {\theta })\ \xi _m^2 e^{-V(\vec {\xi }, \vec {\theta })}}{\int {{\,\mathrm{d\mathrm {\mu }}\,}}_0(\vec {\xi }, \vec {\theta }) e^{-V(\vec {\xi }, \vec {\theta })}} \end{aligned}$$
(24)

where H is defined in (21), \({{\,\mathrm{d\mathrm {\mu }}\,}}_0(\vec {\xi }, \vec {\theta })\) is the product of the measures \({{\,\mathrm{d\mathrm {\mu }}\,}}_0(\xi _i, \theta _i)\) defined in (14) and

(25)

Note that the integral with respect to \({{\,\mathrm{d\mathrm {\mu }}\,}}_0(\vec {\xi }, \vec {\theta })\) must be restricted to “compatible configurations”, that is, those configuration satisfying the hard core compatibility condition.

By writing

$$\begin{aligned} {\prod _{1\le i< j \le N}} e^{-V_{ij}} = {\prod _{1\le i < j \le N}} \left[ (e^{-V_{ij}} - 1) + 1\right] , \end{aligned}$$
(26)

it is possible to rewrite (24) exploiting the following combinatorial identity:

$$\begin{aligned} \prod _{1 \le i < j \le N} (b_{ij} + 1) = \sum _{g \in \mathcal {G}_{N}^{}} \prod _{\{ij\} \in E(g)} b_{ij}, \end{aligned}$$
(27)

where with \(\mathcal {G}_{N}^{}\) we denote the set of all graphs with N vertices and with E(g) the set of all edges of the graph g. Thus we can write

$$\begin{aligned} \langle \xi _m^2\rangle&= \frac{\int {{\,\mathrm{d\mathrm {\mu }}\,}}_0(\vec {\xi }, \vec {\theta })\ \xi _m^2 \prod \nolimits _{i<j} e^{-V_{ij}}}{\int {{\,\mathrm{d\mathrm {\mu }}\,}}_0(\vec {\xi }, \vec {\theta }) \prod \limits _{i<j} e^{-V_{ij}}}\nonumber \\&= \frac{\sum \nolimits _{g \in \mathcal {G}_{N}^{}} \int {{\,\mathrm{d\mathrm {\mu }}\,}}_0(\vec {\xi }, \vec {\theta }) \, \xi _m^2 \prod \nolimits _{\{ij\} \in E(g)} (e^{-V_{ij}} - 1)}{\sum \nolimits _{g \in \mathcal {G}_{N}^{}} \int {{\,\mathrm{d\mathrm {\mu }}\,}}_0(\vec {\xi }, \vec {\theta }) \prod \nolimits _{\{ij\} \in E(g)} (e^{-V_{ij}} - 1)} \end{aligned}$$
(28)

Note that

$$\begin{aligned} \sum _{g \in \mathcal {G}_{N}^{}} \prod _{\{ij\} \in E(g)} b_{ij} = \sum _{k=1}^{N} \sum _{X_1, \ldots , X_k} \prod _{l=1}^{k} \sum _{g \in G_{X_l}} \prod _{\{ij\} \in E(g)} b_{ij} \end{aligned}$$
(29)

where \(X_1, \ldots , X_k\) is a partition of the set \(\{1, \ldots , N\}\) and \(G_{X_l}\) is the set of all connected graphs with vertices in the set \(X_l\).

Using this approach, and denoting by \(X_{0}\) the component of the graph containing the vertex associated to the m-th body, (28) can be rewritten in terms of connected components in the following way:

(30)

that allows to bound \(\langle \xi _m^2\rangle \) as follows (this is what is commonly known as Peierls argument, see for instance [8, 22] for its definition in the case of low-temperature spin systems)

(31)

since

(32)

Indeed, in the previous expression, since the pair potential \(V_{ij}\) are negative in the integration region, the sums are over positive terms, both numerator and denominator are of the same type, but the denominator contains more terms.

We want to rewrite (31) in terms of a sum over trees instead of sum over connected graphs using the so-called Penrose Tree Graph identity introduced by Penrose in [25], see also [7]. To this purpose, denoting \(G_n\) the set of connected graphs on n vertices and \(T_n\) the set of trees, we first give the following:

Definition 2.3

A map \(\mathfrak {M}: T_n \rightarrow G_n\) is called a partition scheme in \(G_n\) if, for all \(\tau \in T_n\), \(\tau \in \mathfrak {M}(\tau )\) and \(G_n = \biguplus _{\tau \in T_n} [\tau , \mathfrak {M}(\tau )]\)

where \(\biguplus \) denotes a disjoint union and \([\tau , \mathfrak {M}(\tau )] = \{g \in G_n : \tau \subset g \subset \mathfrak {M}(\tau )\}\) is a Boolean interval with respect to the set-inclusion. Further, given a partition scheme \(\mathfrak {M}\) and a tree \(\tau \in G_{n}\) write \(m(\tau ) = E(\mathfrak {M}(\tau )) \setminus E(\tau )\) so that, in words, \(m(\tau )\) represents the set of all edges that can be added to \(\tau \) to obtain a connected graph in the Boolean interval \([\tau , \mathfrak {M}(\tau )]\).

With this notation, we have the following:

Lemma 2.4

(General Penrose identity). Let \(n > 2\) and let \(\mathfrak {M}:T_{n} \rightarrow G_{n}\) be a partition scheme in \(G_{n}\). Then

$$\begin{aligned} \sum _{g \in G_{n}} \prod _{\{ij\} \in E(g)} (e^{-V_{ij}} - 1)&= \sum _{\tau \in T_n} \prod _{\{ij\} \in E(\tau )} (e^{-V_{ij}} - 1) \prod _{\{uv\} \in m(\tau )} (e^{-V_{uv}}) \end{aligned}$$
(33)

whose proof is straightforward. Indeed:

Proof

Arguing as in (27), we have

$$\begin{aligned} \sum _{g \in G_{n}} \prod _{\{ij\} \in E(g)} (e^{-V_{ij}} - 1)&= \sum _{\tau \in T_n} \prod _{\{ij\} \in E(\tau )} (e^{-V_{ij}} - 1) \sum _{S \subset m(\tau )} \prod _{\{uv\} \in S} (e^{-V_{uv}} - 1) \end{aligned}$$
(34)
$$\begin{aligned}&= \sum _{\tau \in T_n} \prod _{\{ij\} \in E(\tau )} (e^{-V_{ij}} - 1) \prod _{\{uv\} \in m(\tau )} (e^{-V_{uv}} - 1) + 1 \end{aligned}$$
(35)
$$\begin{aligned}&= \sum _{\tau \in T_n} \prod _{\{ij\} \in E(\tau )} (e^{-V_{ij}} - 1) \prod _{\{uv\} \in m(\tau )} (e^{-V_{uv}}) \end{aligned}$$
(36)

\(\square \)

Using again that \(V_{ij} < 0\) in the integration region, we want to bound suitably the terms of the form \(e^{-V_{ij}} - 1\). To do this we assume a priori that \( |V_{ij}|\le \frac{1}{2}\), we use the fact that \(e^x-1\le \frac{5}{4}x\) for \(0\le x\le \frac{1}{2}\), and then we verify a posteriori that the values of \( |V_{ij}|\) that guarantee the control of the \(\varepsilon -N\)-stability are smaller than \(\frac{1}{2}\). We have

$$\begin{aligned} \langle \xi _{m}^{2}\rangle&\le \sum _{n \ge 1} \sum _{\begin{array}{c} \left| X\right| =n X \ni m \end{array}} \int {{\,\mathrm{d\mathrm {\mu }}\,}}_0(X) \xi _m^2 \,\sum _{g \in G_{n}} \prod _{\{ij\} \in E(g)} (e^{-V_{ij}} - 1) \end{aligned}$$
(37)
$$\begin{aligned}&= \sum _{n \ge 1} \sum _{\begin{array}{c} \left| X\right| =n \\ X \ni m \end{array}} \int {{\,\mathrm{d\mathrm {\mu }}\,}}_0(X) \xi _m^2 \sum _{\tau \in T_n} \prod _{\{ij\} \in E(\tau )} (e^{-V_{ij}} - 1) \prod _{\{uv\}\in m(\tau )} e^{-V_{uv}} \end{aligned}$$
(38)
$$\begin{aligned}&\le \sum _{n \ge 1} \sum _{\begin{array}{c} \left| X\right| =n \\ X \ni m \end{array}} \int {{\,\mathrm{d\mathrm {\mu }}\,}}_0(X) \xi _m^2 \sum _{\tau \in T_n} \prod _{\{ij\} \in E(\tau )} \frac{5}{4}\left| V_{ij}\right| \prod _{\{uv\} \in m(\tau )} e^{-V_{uv}} \end{aligned}$$
(39)

3 Similar Asteroids

Here we consider the case of “asteroids” orbiting with similar radii and similar eccentricities around the star (the asteroids are in the same “belt”). Let N be the total number of asteroids and let \(R_{i} = R,\ \gamma _i=\gamma \) for all \(i \in \{1, \ldots , N\}\). Further let \(a_{i}\) be the diameter of the i-th asteroid and let \(a \le a_{i} \le 2a\).

We want to determine conditions ensuring the \(\varepsilon -N\)-stability of the system. Assume \(V_{ij} \le \frac{1}{2}\). It follows from (39)

$$\begin{aligned} \langle \xi _m^2\rangle \le \langle \xi _m^2\rangle _{0} \left( 1 + \bigg |\bigg | \sum _{n \ge 2} \sum _{\begin{array}{c} \left| X\right| =n\\ X \ni m \end{array}} \sum _{\tau \in T_{n}} \prod _{\{ij\} \in \tau } \frac{5}{4} V_{ij}\prod _{\{ij\}\in m(\tau )} e^{-V_{ij} } \bigg |\bigg |_\infty \right) \end{aligned}$$
(40)

where \(||{\cdot }||_\infty \) is the supremum with respect to feasible configurations and the addend 1 represents the case \(X_{0} = \{m\}\). Denote by \(\delta \) be the (common) density of the asteroids and by \(\delta _{s}\) the density of the star. The dimensionless potential can be written in this case in the form:

$$\begin{aligned} V_{ij} = - (\gamma +1) \frac{m_i m_j}{M(m_i + m_j)} \frac{R}{\left| \vec {x}_i - \vec {x}_j\right| } \end{aligned}$$
(41)

and it is straightforward to verify that

$$\begin{aligned} \left| V_{ij}\right| \le (\gamma +1) \frac{\delta }{\delta _s} \frac{a_{<}^{3}}{a_{>}}\frac{R}{R_s^3} \end{aligned}$$
(42)

with \(a_{<} = \min \{a_{i}, a_{j}\}\) and \(a_{>} = \max \{a_{i}, a_{j}\}\).

By (40), (42) and the previous observations, we have

$$\begin{aligned} \varepsilon \le \sum _{n \ge 2} \sum _{\begin{array}{c} \left| X\right| =n\\ X \ni m \end{array}} \sum _{\tau \in T_n} \left( (\gamma +1) \frac{\delta }{\delta _s} 5 a^{2}\frac{R}{R_s^3}\right) ^{n-1} e^{{n^2 (\gamma +1)\frac{\delta }{\delta _s} 4 a^{2}\frac{R}{R_s^3}}}. \end{aligned}$$
(43)

Note that the fact that in this bound we have used (42) means that we are considering the worst case scenario, i.e., an n-body collision, in order to evaluate the correction to the free measure due to the interaction among asteroids.

Introducing a dependence on N and calling \(A = N ((\gamma +1) \frac{\delta }{\delta _s} 5 a^{2}\frac{R}{R_s^3})\), \({\bar{A}}=\frac{4}{5}A\), (3) can be rewritten as

$$\begin{aligned} \varepsilon&\le \sum _{n\ge 2} \left( {\begin{array}{c}N-1\\ n-1\end{array}}\right) n^{n-2} (\frac{A}{N})^{n-1} e^{{\bar{A}}n} \end{aligned}$$
(44)
$$\begin{aligned}&= \sum _{n \ge 2} \frac{(N-1)(N-2)\cdots (N-n)}{(n-1)!} n^{n-2} \frac{A^{n-1}}{N^{n-1}}e^{{\bar{A}}n} \end{aligned}$$
(45)
$$\begin{aligned}&\le \sum _{n \ge 2} {A}^{n-1} e^{{\bar{A}}n} \end{aligned}$$
(46)
$$\begin{aligned}&\le \sum _{n\ge 2} (Ae^{{\bar{A}}})^{n-1} e^{{\bar{A}}} \end{aligned}$$
(47)
$$\begin{aligned}&\le e^{{\bar{A}}} \frac{Ae^{{\bar{A}}}}{1 - Ae^{{\bar{A}}}} \end{aligned}$$
(48)

since, for \(n\ge 2\)

$$\begin{aligned}&\frac{n^{n-2}}{(n-1)^{n-1}} \le 1, \end{aligned}$$
(49)
$$\begin{aligned}&\quad k! \ge \left( \frac{k}{e} \right) ^k \end{aligned}$$
(50)

and

$$\begin{aligned} \frac{(N-1)(N-2)\cdots (N-n)}{N^{n-1}} < 1. \end{aligned}$$
(51)

Therefore, \(\varepsilon -N\)-stability of the system is guaranteed if \(\frac{A e^{2\bar{A} }}{1 - Ae^{{\bar{A}}}} \le \varepsilon \). Rough numerical estimates show that if \(A<1/5\), then \(\varepsilon \le 2A\).

Note that in order to have \(A=1/5\), \(Na^2\) has to be bounded by a suitable constant. This means that as outlined in the introduction, if N increases the total mass of the asteroids, obviously proportional to \(Na^3\), goes to zero.

4 Asteroids with Power-Law Mass Distribution

Now we consider a more realistic case (see, for instance, [3, 11, 28]): the N asteroids have different masses/diameters (still under the assumption that they have common densities \(\delta \)) and may have different eccentricity. We will keep the assumption \(R_i=R\ \forall i\) because the final estimate on \(\varepsilon \) will represent an upper bound also for a different average radius \(R_i\). Indeed, we are assuming that the deviations around the average radius have always a probabilistic weight of order 1, while a collision among asteroids, that gives the leading contribution in \(V_{ij}\), has a free probability much smaller than 1 if the two average radii \(R_i\) and \(R_j\) are very different.

The distribution of the parameters \(\gamma _i\) appearing in the free measure is also supposed to be not too spread: the eccentricity of the orbits may vary, but the perturbations due to the planets are similar for all the asteroids.

As far as the masses are concerned, we let the diameters of the asteroids satisfy \(a_{\min } \le a_{i} \le a_{\max }\) and we assume the following power-law distribution for their diameters:

$$\begin{aligned} N(> a) = \frac{c}{a^{\nu }} \end{aligned}$$
(52)

where \(N(> a)\) is the number of asteroids with diameter larger than a and c is a suitable constant. This law is assumed for the known asteroids belts in the solar system. To simplify our discussion, in the remainder we will set \(\nu =2\). It will be clear that a different value of \(\nu \) will affect only the constants, provided \(\nu >1\).

We will define the unit of length in order to take \(a_{\min } = 1\) and \(a_{\max } = 2^{L}\) (for some natural number L). Note that \(a_{\min } = 1\) implies \(c = N\). In the applications describing the main asteroid belt in the Solar System the unit will be \(1\,\mathrm {Km}\) (see below).

We partition the asteroids into L classes \(A_{1}, \ldots , A_{L}\). The i-th asteroid belongs to the l-th class if \(2^{l-1} \le a_{i} < 2^{l}\). In this case we write \(i \in A_{l}\). Denoting by \(N_{l}\) the number of asteroids in the l-th class we have \(N_{l} = N(> 2^{l-1}) - N(> 2^{l}) = \frac{3}{4^{l}}N\).

Let \(i \in A_{l}\) and \(j \in A_{m}\) with \(l > m\). From (19), setting \(\gamma =\max _{1\le i\le N}\gamma _i\), it follows

$$\begin{aligned} \left| V_{ij}\right| \le w_{lm} := (\gamma +1) \frac{\delta }{\delta _{s}} \frac{a_{m}^{3}}{a_{l}} \frac{R}{R_{s}^{3}} = (\gamma +1) \frac{\delta }{\delta _{s}} 4^{m} 2^{-(l-m-1)} \frac{R}{R_{s}^{3}} \end{aligned}$$
(53)

For \(l=m\) we have

$$\begin{aligned} w_{ll}= (\gamma +1) \frac{\delta }{\delta _{s}} 4^{l} \frac{R}{R_{s}^{3}} \end{aligned}$$
(54)

Hence

$$\begin{aligned} \varepsilon \le \sum _{n \ge 2} \sum _{\begin{array}{c} n_1, \ldots n_l\\ \sum _i n_i = n \end{array}} \prod _{l=1}^{L} \left( {\begin{array}{c}N_l\\ n_l\end{array}}\right) n^{n-2} \max _{\tau } \prod _{\{ij\} \in \tau } (e^{w_{lm}} - 1) \prod _{\{ij\} \in m(\tau )} e^{w_{lm}}. \end{aligned}$$
(55)

Since the estimates of the interactions \(w_{lm}\) decay exponentially in \(|l-m|\), the worst case is the tree \(\tau \) having \(n_L - 1\) connections among asteroids in class L and \(n_l\) connections among asteroids of class l. Then

$$\begin{aligned} \max _{\tau } \prod _{\{ij\} \in \tau } (e^{w_{lm}} - 1)&\le \left( (\gamma +1) \frac{\delta }{\delta _s} \frac{R}{R_s^3}4^L \right) ^{n_L - 1}\prod _{l=1}^{L-1} \left( (\gamma +1) \frac{\delta }{\delta _s} \frac{R}{R_s^3}4^l \right) ^{n_l} \end{aligned}$$
(56)
$$\begin{aligned}&\le A^{n-1}\left( \prod _{l=1}^{L-1} \left( \frac{1}{N_l}\right) ^{n_l}\right) \left( \frac{1}{N_L}\right) ^{n_L - 1} \end{aligned}$$
(57)

where in the last step we have set \(A = \gamma \frac{\delta }{\delta _s} \frac{3R}{R_s^3}N\).

In addition, we have

(58)

Finally we obtain (assuming the asteroid m is in class 1):

$$\begin{aligned} \varepsilon&\le \sum _{n \ge 2} \sum _{\begin{array}{c} n_1 \ldots n_l\\ \sum _i n_i = n\\ n_i \ge 1 \end{array}} \prod _{l=2}^{L} \left( {\begin{array}{c}N_l\\ n_l\end{array}}\right) \left( {\begin{array}{c}N_1 - 1\\ n_1 - 1\end{array}}\right) n^{n-2} A^{n-1} \prod _{l=1}^{L-1} \left( \frac{1}{N_l} \right) ^{n_l}\left( \frac{1}{N_L} \right) ^{n_L -1} e^{nA} \end{aligned}$$
(59)
$$\begin{aligned}&\le e^{A+1}\sum {(ALe^{A})}^n = e^{A+1}\frac{ALe^{A}}{1 - ALe^{A}} \end{aligned}$$
(60)

Remark 4.1

It is interesting, in this slightly more realistic framework, to compare this result with the actual main belt of asteroids of Solar System. The parameter A, setting \(\gamma =50\), \(\frac{\delta }{\delta _s}=2\) and the real values for R, \(R_s\), has a value \(A\approx \frac{N}{5\times 10^5}\). Setting \(L=10\), and considering only the asteroids with a diameter \(a\ge 1\, \mathrm {Km}\), one finds that to obtain \(\varepsilon \le 1\) the condition on A is \(A\le 1/4\). This means that with our (rough) approximations \(N\approx 10^5\). The actual number of asteroids in the main belt having diameter larger that \(1\, \mathrm {Km}\) is \(N=10^6\)

Remark 4.2

The computation above assumes a minimal size of the asteroids. Here we present an indication of the fact that the power law mass distribution for the very light asteroids has to have an exponent \(\nu <1\). Calling dN(a) the number of asteroids having the diameter between a and \(a+da\), we clearly have that if \(N(> a) = \frac{N_1}{a^{\nu }}\), then

$$\begin{aligned} dN(a)=N_1\nu \frac{da}{a^{\nu +1}} \end{aligned}$$
(61)

Considering that in the estimate of \(\varepsilon \) we have to give a bound of the quantity \(\sum _{ij}|V_{ij}|\) and using (42) and (61) we get

$$\begin{aligned} \sum _{ij}|V_{ij}|\le (\gamma +1) \frac{\delta }{\delta _s} \frac{R}{R_s^3}\int _{a_\mathrm{min}}^{a_\mathrm{max}}da \int _{a}^{a_\mathrm{max}}db N_1^2 \nu ^2\frac{1}{a^{\nu +1}}\frac{1}{b^{\nu +1}}\frac{a^3}{b} \end{aligned}$$
(62)

It is now clear that \(a_\mathrm{max}\) has to be such that \(N(>a_\mathrm{max})=1\), and hence \(a_\mathrm{max}=N_1^{\frac{1}{\nu }}\). Performing the elementary integrals in (62), we get

$$\begin{aligned} \sum _{ij}|V_{ij}|\le (\gamma +1) \frac{\delta }{\delta _s} \frac{R}{R_s^3}\frac{\nu ^2N_1^2}{\nu +1} \left[ \frac{N_1^{\frac{2-2\nu }{\nu }}-a_\mathrm{min}^{2-2\nu }}{ 2-2\nu }- \frac{N_1^{\frac{2-2\nu }{\nu }}-a_\mathrm{min}^{3-\nu }N_1^{-\frac{\nu +1}{\nu }}}{3-\nu }\right] \nonumber \\ \end{aligned}$$
(63)

This expression shows that if \(a_\mathrm{min}\) is a finite value, say \(a_\mathrm{min}=1\), then the conditions ensuring the control of \(\varepsilon \) are \(\nu >1\) and

$$\begin{aligned} (\gamma +1) \frac{\delta }{\delta _s} \frac{R}{R_s^3}\frac{\nu ^2N_1^2}{\nu +1}<K \end{aligned}$$
(64)

with K suitably chosen. If we want to consider small \(a_\mathrm{min}\), we have to assume that the distribution \(N(> a) = \frac{N_1}{a^{\nu }}\), with \(\nu >1\), is valid for \(a>1\), while defining \(N_<(>a)\) as the number of asteroids having a diameter between a and 1, it has to be of the form \(N_<(> a) = \frac{N_1}{a^{\nu '}}\) with \(\nu '<1\). To our knowledge we do not have many observations on the mass distribution of the very small asteroids. However experimental data seem to show (see, for instance, [28]) that the exponent in the distribution tends to decrease for smaller asteroids.

5 Planets

The basic idea developed in the previous sections is to describe the effect of the perturbation given by other distant objects, say planets, to the orbits of a large number of asteroids living in a single belt, i.e., with similar radii, by a probability distribution centered around a circular orbit. In this section we try to apply the same idea to a system of relatively few planets having well separated orbits. In this case the free measure, i.e., the system obtained neglecting the interactions with the other planets, can be completely determined in terms of an elementary two-body problem. However we shall see that a toy model in which N, the number of planets, is small (\(\sim 10\)), the masses of the planets may be quite different and the eccentricity of the orbits is very small (large \(\gamma _i\)) for all planets, keeps some interesting forecast performance, even when we substitute the well-known Keplerian orbit with a probability distribution.

The computations involved in this case, however, are quite different. Indeed, in the asteroids case the quantity to be controlled is the probability of collisions, and such collisions do not imply large deviations, in terms of the free measure, from the average value R of the distance from the star, that is the same for all asteroids. In other words, the detailed structure of the free measure does not play any role, and the Gaussian approximation of the free measure is simply a way to compute very easily the free variance of the distribution of the distance from the star. The estimates, therefore, can be done always in the sense of an \(L_\infty \) norm, and the fact that with a reasonable choice of the parameters we can keep \(\varepsilon \) small means that the collisions give a negligible contribution to the interacting probability.

In the case of planets we will show that the system is \(\varepsilon -N\)-stable if the radii \(R_i\) of the planets are very different, namely if the condition \(R_{i} - R_{j} = c (a^{i} - a^{j})\) holds. This assumption amounts to saying that the radii of the orbits satisfy the Titius–Bode law, that is, \(R_{i} = b + c a^{i}\). Note that the Titius–Bode law is fulfilled quite accurately in the Solar System.

Since our main task is to verify that even in this case, with larger masses, the collisions give a negligible contribution to the interacting measure, we have to modify the previous computations: collisions are events with a very small probability with respect to the free measure, and hence we cannot use \(L_\infty \) estimates in order to evaluate the collisions. On the other side, since a collision is possible only when at least one planet has a very large fluctuation around his free orbit, the Gaussian approximation loses its meaning, and we need some initial estimates about the free complete measures.

The first important observation is that the probability density

$$\begin{aligned} dw_0(\xi ,\theta )=\exp \left( -\frac{\gamma ^2}{2}\frac{\xi ^2}{(1+\xi )^2} \right) d\theta d\xi \end{aligned}$$
(65)

cannot be normalized on the whole space. Indeed

$$\begin{aligned} \int _0^{2\pi }d\theta \int _{-1}^\infty d\xi \exp \left( -\frac{\gamma ^2}{2}\frac{\xi ^2}{(1+\xi )^2} \right) =\infty \end{aligned}$$
(66)

The simplest way out is to define the free measure on a finite space, say on a sphere of radius \(2R_N\). This means that for the i-th planet \(-1<\xi _i\le A_i=\frac{2R_N}{R_i}\). Since our task is to show that the collisions among planets have a negligible probability in the interacting measure, we will show that for large \(\gamma \) the main contribution to the interacting measure will be given by the configurations in which each planet i will have a distance from the star quite close to \(R_i\), i.e., a \(\xi _i\) of the order of \(1/\gamma \). In other words, we are saying that a planet is inside the planetary system if it is not too far from the star. Note that in the Solar System \(A_\mathrm{Mercury}=200\). Hence we will call

$$\begin{aligned} Z_i=\int _0^{2\pi }d\theta _i\int _{-1}^{A_i}d\xi _i \exp \left( -\frac{\gamma _i^2}{2}\frac{\xi _i^2}{(1+\xi _i)^2} \right) \end{aligned}$$
(67)

Note that the main contribution in the integral comes from the interval \(-1/2\le \xi _i\le 1/2\) since the obvious \(L_\infty \) estimate

$$\begin{aligned} \int _{|\xi |>1/2}d\xi _i \exp \left( -\frac{\gamma _i^2}{2}\frac{\xi _i^2}{(1+\xi _i)^2} \right) \le A e^{-\frac{\gamma _i^2}{18}} \end{aligned}$$
(68)

holds. It is a standard algebraic task, then, to show that for large values of \(\gamma \) the variance of \(\xi _i\) of the free measure for all planets i is proportional to \(\sigma _-^2=\frac{1}{4\gamma _i}\), as in the Gaussian approximation. Note that

$$\begin{aligned} \langle \xi _{i}^{2}\rangle _{0}:=\int {{\,\mathrm{d\mathrm {\mu }}\,}}_0(\xi _i) \xi ^{2} :=\frac{1}{Z_i} \int _0^{2\pi }d\theta _i\int _{-1}^{A_i}d\xi _i\ \xi _i^2\exp \left( -\frac{\gamma _i^2}{2}\frac{\xi _i^2}{(1+\xi _i)^2} \right) \end{aligned}$$
(69)

and hence

$$\begin{aligned} \langle \xi _{i}^{2}\rangle _{0}=\frac{1}{Z_i} \int _0^{2\pi }d\theta _i\int _{|\xi _i|\le 1/2}d\xi _i\ \xi _i^2\exp \left( -\frac{\gamma _i^2}{2}\frac{\xi _i^2}{(1+\xi _i)^2} \right) + O \left( A_i^3e^{-\frac{\gamma _i^2}{18}} \right) \nonumber \\ \end{aligned}$$
(70)

The interesting values for application to the Solar system are \(A_i\le 200\) and \(\gamma _i\ge 50\), and therefore the correction of order \(A^3e^{-\frac{\gamma _i^2}{18}}\) is completely negligible. Hence the leading part of the integral is dominated from above and from below by two Gaussian measures with variance \(\sigma _-^2=\frac{1}{4\gamma _i}\) and \(\sigma _+^2=\frac{9}{4\gamma _i}\), respectively, and these are two bounds, both proportional to \(\gamma _i^{-1}\), for the variance.

In what follows we will assume for simplicity that \(\gamma _i=\gamma \) for all planets. For all \(m=1,...,N\) we want to give an estimate of the quantity

$$\begin{aligned} \varepsilon \langle \xi _{m}^{2}\rangle _{0} = \sum _{n = 2}^{N} \sum _{\begin{array}{c} \left| X\right| =n\\ X \ni x_{m} \end{array}} \int {{\,\mathrm{d\mathrm {\mu }}\,}}_0(X) \xi ^{2} \sum _{\tau \in T_{n}} \prod _{\{ij\} \in \tau } (e^{-V_{ij}}-1)\prod _{\{ij\}\in m(\tau )} e^{-V_{ij}} \end{aligned}$$
(71)

where \(V_{ij} = -2\gamma \frac{m_im_j}{M} \frac{R_iR_j}{R_im_j+R_jm_i}\frac{1}{\left| \vec {x_i} - \vec {x_j}\right| }\)

We call \(\xi \) typical when \(|{\xi _i}| < k\frac{1}{\gamma }\). For a fixed X, we write \(X = T \bigcup T^{c}\) with \(T = \{i \in X | \xi _i \text { is typical}\}\).

Let us consider first the case in which i and j are both typical. Standard algebra shows that for \(j>i\)

$$\begin{aligned} \left| \vec {x_i} - \vec {x_j}\right|&\ge R_{j}\left( 1 - \frac{k}{\gamma } \right) - R_{i}\left( 1 + \frac{k}{\gamma } \right) \end{aligned}$$
(72)
$$\begin{aligned}&\ge c_1 (a^j - a^i) \end{aligned}$$
(73)

with \(c_1 = c - \frac{2k}{\gamma } (c + b)\)

On the other side, if i and/or j are not typical

$$\begin{aligned} \left| \vec {x_i} - \vec {x_j}\right| \ge r_j+r_i \end{aligned}$$
(74)

where \(r_i\) is the radius of planet i. To obtain an upper bound of \(V_{ij} \) we observe that, recalling \(j>i\)

$$\begin{aligned} \frac{R_iR_j}{R_im_j+R_jm_i}\le \frac{b+ca^i}{m_{ij}^\mathrm{min}} \end{aligned}$$
(75)

where obviously \(m_{ij}^\mathrm{min}\) is the smallest mass between the planet i and j. If \(i_\mathrm{min}\) is the smallest i in the planetary system (recall for instance that in the Solar System Mercury corresponds to \(i=-1\)), calling \(c_2=c+\frac{b}{a^{i_\mathrm{min}}}\) we obtain

$$\begin{aligned} \frac{R_iR_j}{R_im_j+R_jm_i}\le \frac{c_2a^i}{m_{ij}^\mathrm{min}} \end{aligned}$$
(76)

Then if i and j are both typical, \(j>i\), we get

$$\begin{aligned} |V_{ij}|\le 2\gamma \frac{{m_{ij}^\mathrm{max}}}{M}\frac{c_2}{c_1}\frac{1}{a^{j-i}-1} \le \frac{2a\gamma }{a-1}\frac{m_{ij}^\mathrm{max}}{M}\frac{c_2}{c_1}a^{-(j-i)}:={{\bar{V}}}_{ij} \end{aligned}$$
(77)

Otherwise

$$\begin{aligned} |V_{ij}|\le {{\bar{V}}}_{ij}\frac{c_1(a^j-a^i)}{r_j+r_i}= {{\widetilde{V}}}_{ij} \end{aligned}$$
(78)

The strategy will be the following: first we evaluate in (71) the case \(X=T\), using for \(V_{ij}\) the estimate \({{\bar{V}}}_{ij}\). We will call \(\bar{\varepsilon }\) the estimate obtained in this way. In this case we will proceed as in the previous cases, with an \(L_\infty \) estimate.

$$\begin{aligned} {\bar{\varepsilon }} \le&\sum _{n = 2}^{N} \sum _{\begin{array}{c} \left| X\right| =n\\ X \ni x_{m} \end{array}} \sum _{\tau \in T_{n}} \prod _{\{ij\} \in \tau } (e^{{{\bar{V}}}_{ij}}-1)\prod _{i<j}e^{{{\bar{V}}}_{ij}} \end{aligned}$$
(79)
$$\begin{aligned} \le&\left[ \prod _{i\ne m}\left( 1+\sum _{j\ne i}\left( e^{{{\bar{V}}}_{ij}}-1\right) \right) -1\right] \prod _{i<j}e^{{{\bar{V}}}_{ij}} \end{aligned}$$
(80)

Remark 5.1

Note that in (79) we gave a quite rough estimate of the combinatorics on trees. In particular we used that for trees rooted in m, since every vertex but m has a unique “predecessor”, \(\sum _{\tau _X}\prod _{ij\in \tau }e_{ij}\le \prod _{i\ne m}\sum _{j\ne i}e_{ij}\). The addend 1 takes into account the sum on X. The last addend \(-1\) takes into account the fact that the sum in n starts from 2. Since N is small and the orbits are well separated, this estimate is reasonable.

Calling now \(c_3=\frac{2a\gamma }{a-1}\frac{m_{ij}^\mathrm{max}}{M}\frac{c_2}{c_1}\) we have \({{\bar{V}}}_{ij}\le c_3a^{-(j-i)}\). Assume \(c_3<1/2\). Since \(\sqrt{e}<5/3\) we have that \(e^{{{\bar{V}}}_{ij}}-1<\frac{5}{3}{{\bar{V}}}_{ij}\) and hence

$$\begin{aligned} {\bar{\varepsilon }} \le \left[ \prod _{i\ne m}\left( 1+\frac{5}{3}c_3\sum _{j\ne i}a^{-|j-i|}\right) -1\right] e^{c_3\sum _{i< j}a^{-(j-i)}} \end{aligned}$$
(81)

Using now the elementary inequalities \(\sum _{i< j}a^{-(j-i)}\le \frac{1}{a-1}\) and \(1+x\le e^x\) we finally get

$$\begin{aligned} {\bar{\varepsilon }} \le \left( e^{\frac{5}{3}c_3N\frac{2}{a-1}}-1\right) e^{N\frac{c_3}{a-1}} \end{aligned}$$
(82)

This concludes the estimate for \(T=X\). The crucial relation to control the general case is the following. Call \({\widetilde{E}}\) the set of pairs i, j of planets such that for their estimate we cannot use (77) (collisions). For a fixed \(T^c\) the contribution to \(\varepsilon \), that we will denote \(\varepsilon (T^c)\) can be bounded by

$$\begin{aligned} \varepsilon (T_c)\le \int d\mu _0(T^c)\ \prod _{ij\in {{\widetilde{E}}}}e^{{{\widetilde{V}}}_{ij}}(e^{{{\bar{V}}}_{ij}}-1)^{-1} \end{aligned}$$
(83)

To prove (82) it is enough to observe that for all \(\tau \)

$$\begin{aligned} \prod _{\{ij\} \in E(\tau )} (e^{-V_{ij}}-1) \prod _{\{ij\}\in m(\tau )} e^{-V_{ij} }\le \prod _{\{ij\} \in E(\tau )} (e^{{{\bar{V}}}_{ij}}-1) \prod _{\{ij\}\in m(\tau )} e^{{{\bar{V}}}_{ij} }\prod _{ij\in {{\widetilde{E}}}}\frac{e^{{{\widetilde{V}}}_{ij}}}{(e^{{{\bar{V}}}_{ij}}-1)} \end{aligned}$$
(84)

and then bound with 1 the contribution of the integral \(\int d\mu _0(T)\)

The idea is then to bound the very large contribution due to \(\prod _{ij\in {{\widetilde{E}}}}e^{{{\widetilde{V}}}_{ij}}(e^{{{\bar{V}}}_{ij}}-1)^{-1}\) with the smallness of \(\int d\mu _0(T^c)\).

Let us start with the simplest case in which \(T^c=\{i\}\) and the collision is with planet \(i+1\). In this case the only estimate we can do for the probability of collision with respect to the free measure is the probability of i to be non-typical, \(\mu _0(T^c)\le e^{-\frac{2k^2}{9}}\). On the other side the weight in the interacting measure of the collision is proportional to \(e^{{{\bar{V}}}_{i,i+1}\frac{c_1a^i(a-1)}{r_i+r_{i+1}}}\approx e^{c_3\frac{c_1a^i(a-1)}{a(r_i+r_{i+1})}}\). Hence our condition in order to control the single collisions will be

$$\begin{aligned} \frac{2k^2}{9}>c_3\frac{c_1a^i(a-1)}{a(r_i+r_{i+1})}>\frac{2a\gamma }{a-1}\frac{m_{ij}^\mathrm{max}}{M}c_2\frac{a^i}{r_i+r_{i+1}} \end{aligned}$$
(85)

We outline that for N not too large, say \(N\ge 10\), the case \(T^c=\{i\}\) is the leading one: in order to evaluate the l-body collisions the contribution \(c_3\frac{c_1a^i(a-1)}{a(r_i+r_{i+1})}\) has to be multiplied by \({l\atopwithdelims ()2}\), while the contribution \(\frac{2k^2}{9}\) becomes much larger, because at least \(l-2\) planets have to undergo a deviation in \(\xi \) of order 1, and hence the factor becomes of the order of \(\gamma ^2\) instead of \(k^2\).

We end this section outlining that (85) and (83) can be specified in the case of the planets of the Solar System and in the case of the Galilean satellites. Note, however, that the numerical estimates we stated in the generic case may be specified better once we know the actual value of the parameter. In the case of the planets we solve simply both conditions, in the sense of the 1-stability, using as free parameter \(m^\mathrm{max}\). k can be chosen in order to have the largest possible value of \(m^\mathrm{max}\). Reasonable values for the parameters are:

  • \(\gamma =150\), since the eccentricity of the orbits are very small.

  • \(k=30\),

  • \(c_2=1 UA\)

  • \(a=2\)

  • \(a^{i^\mathrm{max}}= 128\)

Then it is possible to satisfy (85) and (82) in order to have \(\varepsilon <1\) with a value of \(m^\mathrm{max}\) similar to the Earth’s one. In the case of Galilean satellites, in which \(N=4\) and, most of all, \(b=0\), we can control the various steps of the estimates much better. The combinatorics on trees and the sums on \(V_{ij}\) can be written more explicitly, obtaining eventually that a ratio \(m^\mathrm{max}/m_J\approx 10^{-4}\), which is the actual value, ensures 1-stability.

6 Conclusions and Open Problems

The aim of this work is to outline the fact that with a judicious but quite standard use of results typical of equilibrium statistical mechanics one can evaluate some global features of the systems of particles rotating around a much bigger body. The estimates presented here are quite rough, and they can be surely improved by a careful numerical evaluation of the constants appearing in the theory. Nevertheless, the results we got, namely an evaluation of the “thermodynamical” stability of the main asteroid belt, of the planets in the solar system and of the Galilean satellites, give quantitative estimates not too distant from the real data, and seem therefore to indicate that this approach to the planetary system gives a reasonable possibility to understand the global structure of the Solar System. More precisely, our model seems to indicate that in order to have a thermodynamically stable system the masses of the particles orbiting around the fixed large body have to be very small if the orbit’s parameters of the particle are similar, but they can increase if the objects are far apart. It would be nice to have some data about the very small objects in the belts of the Solar system (main belt of asteroids, trans-Neptunian belts, rings around the planets) because our model seems to indicate that the distribution of the very light objects in a belt has to have a different scaling law with respect to the one of the heavier ones.